PARI's commands

[science_man_88]

Originally Posted by axn

? n=0
%12 = 0
? ispower(1024,,&n)
%13 = 10
? n
%14 = 2

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

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]

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]

Originally Posted by CRGreathouse

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.

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

[axn]

Originally Posted by science_man_88

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

use

b^valuation(n,b)==n

[science_man_88]

Originally Posted by axn

use

b^valuation(n,b)==n

or just:

(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

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]

Originally Posted by science_man_88

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

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

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));
}

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

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;
}

which is nearly as fast for small inputs, maybe 25% slower.

[CRGreathouse]

Originally Posted by science_man_88

or just:

(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

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%

Your code isn't testing if the numbers are powers of two, just if they're odd. Shorter version:

sum(n=1,#mersenne1,((2^mersenne1[n]-1)%n)%2)

[science_man_88]

Originally Posted by CRGreathouse

Your code isn't testing if the numbers are powers of two, just if they're odd. Shorter version:

sum(n=1,#mersenne1,((2^mersenne1[n]-1)%n)%2)

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

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)

using this on the primes I get:

(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

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

[CRGreathouse]

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

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.");

(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]

Originally Posted by CRGreathouse

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

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.");

(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.

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]

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