mersenneforum.org

mersenneforum.org (https://www.mersenneforum.org/index.php)
-   PARI/GP (https://www.mersenneforum.org/forumdisplay.php?f=155)
-   -   PARI's commands (https://www.mersenneforum.org/showthread.php?t=13636)

science_man_88 2011-02-10 19:11

[QUOTE=axn;252089][CODE]? n=0
%12 = 0
? ispower(1024,,&n)
%13 = 10
? n
%14 = 2[/CODE]

You must pass a variable as the third parameter so that the function can pass back the "base".[/QUOTE]

thanks for the tip I never quite figured it out I asked because i was seeing something that I sent to prime95 as weird, now I don't as much but if you look (2^mersenne1[n]-1)%n seems to be a mersenne number a lot if mersenne1 is the mersenne prime exponents. in fact it looks like 35/40 have this property in the first 100 prime numbers if you do 2^(prime(n))-1 over there index into the primes it seems about 49% of the first 100 work out.

CRGreathouse 2011-02-10 21:37

ispower(n) checks if n is a nontrivial power. ispower(n,k) checks if n is a k-th power. ispower(n,,&b) checks if n is a nontrivial power and, if so, stores the base in b. ispower(n,k,&b) checks if n is a k-th power and, if so, stores the base in b.

science_man_88 2011-02-10 21:40

[QUOTE=CRGreathouse;252111]ispower(n) checks if n is a nontrivial power. ispower(n,k) checks if n is a k-th power. ispower(n,,&b) checks if n is a nontrivial power and, if so, stores the base in b. ispower(n,k,&b) checks if n is a k-th power and, if so, stores the base in b.[/QUOTE]

it should be able to check if n is a power of b.

axn 2011-02-10 22:19

[QUOTE=science_man_88;252113]it should be able to check if n is a power of b.[/QUOTE]

use [CODE]b^valuation(n,b)==n[/CODE]

science_man_88 2011-02-10 22:56

[QUOTE=axn;252115]use [CODE]b^valuation(n,b)==n[/CODE][/QUOTE]

or just:

[CODE](18:55)>a=0;for(n=1,#mersenne1,if(valuation(((2^mersenne1[n]-1)%n)+1,2)!=0,a=a+1));return(a)
%51 = 35[/CODE] to check my count, which is accurate, for the amount in mersenne1 . my count wasn't as accurate for the primes and the indexes withing the primes. looks like from primes 1 to 100 we have 75 that fit my criteria not 67.5%

CRGreathouse 2011-02-10 23:57

[QUOTE=science_man_88;252113]it should be able to check if n is a power of b.[/QUOTE]

I've needed that too, at least in the case of powers of two. axn's solution is good. For a faster method, I implemented
[code]long
ispow2(GEN n)
{
if (typ(n) != t_INT)
pari_err(arither1, "ispow2");
if (signe(n) < 1)
return 0;
pari_sp ltop = avma;

GEN xp = int_LSW(n);
long lx = lgefint(n);
ulong u = *xp;
long i = 3;
for (; i < lx; ++i)
{
if (u != 0)
{
avma = ltop;
return 0;
}
xp = int_nextW(xp);
u = *xp;
}
avma = ltop;
return !(u & (u-1));
}[/code]
in C and load it into gp. For small numbers it takes ~11.6 ns vs. 205 ns for axn's version or 175 ns for the shifting variant. For large numbers, especially those divisible by large powers of two, the difference is greater.

Of course, this is comparing apples and oranges. The internal form of axn's algorithm (when appropriately optimized) is
[code]long
ispow2(GEN n)
{
if (typ(n) != t_INT)
pari_err(arither1, "ispow2");
if (signe(n) < 1)
return 0;
pari_sp ltop = avma;

long ret = expi(n) == vali(n);
avma = ltop;
return ret;
}[/code]
which is nearly as fast for small inputs, maybe 25% slower.

CRGreathouse 2011-02-11 00:14

[QUOTE=science_man_88;252118]or just:

[CODE](18:55)>a=0;for(n=1,#mersenne1,if(valuation(((2^mersenne1[n]-1)%n)+1,2)!=0,a=a+1));return(a)
%51 = 35[/CODE] to check my count, which is accurate, for the amount in mersenne1 . my count wasn't as accurate for the primes and the indexes withing the primes. looks like from primes 1 to 100 we have 75 that fit my criteria not 67.5%[/QUOTE]

Your code isn't testing if the numbers are powers of two, just if they're odd. Shorter version:
[code]sum(n=1,#mersenne1,((2^mersenne1[n]-1)%n)%2)[/code]

science_man_88 2011-02-11 00:42

[QUOTE=CRGreathouse;252128]Your code isn't testing if the numbers are powers of two, just if they're odd. Shorter version:
[code]sum(n=1,#mersenne1,((2^mersenne1[n]-1)%n)%2)[/code][/QUOTE]

okay it still gave a correct answer, here's the version you seem to want:

[CODE]a=0;for(n=1,#mersenne1,if(2^valuation((((2^mersenne1[n])-1)%n)+1,2)==(((2^mersenne1[n])-1)%n)+1,a=a+1));return(a)[/CODE]

using this on the primes I get:

[CODE](20:38)>a=0;for(n=1,100,if(2^valuation((((2^prime(n))-1)%n)+1,2)==(((2^prime(n))-1)%n)+1,a=a+1));return(a)
%107 = 50
(20:39)>a=0;for(n=1,1000,if(2^valuation((((2^prime(n))-1)%n)+1,2)==(((2^prime(n))-1)%n)+1,a=a+1));return(a)
%108 = 256
(20:39)>a=0;for(n=1,10000,if(2^valuation((((2^prime(n))-1)%n)+1,2)==(((2^prime(n))-1)%n)+1,a=a+1));return(a)
%109 = 971[/CODE]

each time as you get higher fewer and fewer fit what 35/40 of the Mersenne prime exponents appear to fit.

CRGreathouse 2011-02-11 02:15

Incidentally, since you're working with Mersenne primes and exponents a lot, here's an idea for a script for your file:

[code]MeVec=[2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049, 216091, 756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593, 13466917];
Me(n)={
if(n<=#MeVec,
MeVec[n]
,
error("not known")
)
};
addhelp(Me, "Me(n): Returns the n-th Mersenne prime exponent, Sloane's A000043.");

M(n,modulus=0)={
if(n>#MeVec, error("not known"));
if(modulus,
lift(Mod(2,modulus)^MeVec[n]-1)
,
1<<MeVec[n]-1
)
};
addhelp(M, "M(n, {modulus}): Returns the n-th Mersenne prime, Sloane's A000668. If modulus is given, instead return the n-th Mersenne number mod the modulus.");[/code]

(Of course you can expand the vector to include the numbers whose primality is known but whose order is not.)

This should not only save you keystrokes but time on longer calculations since the modular calculations are efficient compared to creating a large number and dividing it.

science_man_88 2011-02-11 19:21

[QUOTE=CRGreathouse;252145]Incidentally, since you're working with Mersenne primes and exponents a lot, here's an idea for a script for your file:

[code]MeVec=[2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049, 216091, 756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593, 13466917];
Me(n)={
if(n<=#MeVec,
MeVec[n]
,
error("not known")
)
};
addhelp(Me, "Me(n): Returns the n-th Mersenne prime exponent, Sloane's A000043.");

M(n,modulus=0)={
if(n>#MeVec, error("not known"));
if(modulus,
lift(Mod(2,modulus)^MeVec[n]-1)
,
1<<MeVec[n]-1
)
};
addhelp(M, "M(n, {modulus}): Returns the n-th Mersenne prime, Sloane's A000668. If modulus is given, instead return the n-th Mersenne number mod the modulus.");[/code]

(Of course you can expand the vector to include the numbers whose primality is known but whose order is not.)

This should not only save you keystrokes but time on longer calculations since the modular calculations are efficient compared to creating a large number and dividing it.[/QUOTE]

I tried something with this can it be confirmed that all Mersenne primes so far found>31 fit (31 or 43) mod 84 I bet I'm missing something trivial. If not can we limit the prime exponents by this ( did test saves 2 prime exponent checks under 1000000 both under 1000)?

CRGreathouse 2011-02-11 21:59

Can you prove that it behaves the way you expect mod 7? Mod 3? Mod 4?


All times are UTC. The time now is 23:02.

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.