PARI's commands

[science_man_88]

Originally Posted by CRGreathouse

Does (4^4655 - 1)/3 fail also?

no. edit: but for me it fails at a(4651)

[CRGreathouse]

Originally Posted by science_man_88

no. edit: but for me it fails at a(4651)

So the problem is with your function a(). You'll need to rewrite it to not fill the stack.

[science_man_88]

 = ()->for(k=0,2^199,if((isprime(floor(log(6*k+2)/log(2)))||(issquare(floor(log(6*k+2)/log(2)))&&isprimepower(floor(log(6*k+2)/log(2)))==2
&&isprime(2^floor(sqrt(log(6*k+2)/log(2)))-1)))&&#factor(6*k+1)[,1]==2&&sum(x=1,#factor(6*k+1)[,2],factor(6*k+1)[x,2])==2,print1(floor(log(6*k+2)/log(2))","));k=4*k)

I tried different versions like this before just going back through my pari history but I can't get much improvement if at all in fact some things I thought could be improvements don't appear to be like pre-calculating

floor(log(6*k+2)/log(2))

seems to not really help much and I have been debating if I should use issquarefree(x) or core(x)==x because one seems faster above a certain number once I take the limit away, I even tried k as the exponent and a forstep and nothing seems to improve it even though I know there are checks I could do to improve it more like the squarefree scenario ( of course I can limit mod 120 etc first to lower the amount of squarefree checks when checking the prime part of the OR) any other ideas that might speed it up

[CRGreathouse]

Since your function is an unreadable single line I won't dissect it, but I notice that you're using floor(log(6*k+2)/log(2)) a lot. Let's define t = floor(log(6*k+2)/log(2)), which you'll notice stays the same for long streches of k and then goes up by 1. Why don't you loop over values of t directly, and then re-create the associated values of k if t passes your tests?

[science_man_88]

Originally Posted by CRGreathouse

Since your function is an unreadable single line I won't dissect it, but I notice that you're using floor(log(6*k+2)/log(2)) a lot. Let's define t = floor(log(6*k+2)/log(2)), which you'll notice stays the same for long streches of k and then goes up by 1. Why don't you loop over values of t directly, and then re-create the associated values of k if t passes your tests?

each loop through the for loop k is multiplied by 4 and then adds one more. I'm looping through k values that produce odd exponent mersenne numbers I took care of it in my newest version it still doesn't trim a lot of time like i thought. I replaced it with p=logint(6*k+2,2) for each k used and more but I realize I used %39 to call it so in my new session I can't call it back to the quickest form I got to:

try() ={
n=0;
m=vector(120,n,gcd(n-1,120)==1);
for(k=0,2^199,
     p=logint(6*k+2,2);
     if(
        (
        (
        m[(p%120)+1] && core(p)==p && isprime(p)
        )
        ||
        (
        issquare(p) && isprimepower(p,&n)==2 && isprime(2^n-1)
        )
        )
        && #factor(6*k+1)[,1]==2 && sum(x=1,2,factor(6*k+1)[x,2])==2,
        print1(p",")
        )
;        k=4*k
    )}

is the best I could do at rewriting my best implementation.

[CRGreathouse]

OK, that gives me a lot more to go on.

First, you should always declare scope for your variables. In about 100% of cases this should be done with my(), with some residual cases using local() instead. Otherwise you overwrite global variables which is BAD and will cause unintended side-effects down the road.

try()={
	my(m=vector(120,n,gcd(n-1,120)==1), q);
	for(k=0,2^199,
		my(p=logint(6*k+2,2));
		if(
		(
			(m[(p%120)+1] && core(p)==p && isprime(p))
		||
			(issquare(p) && isprimepower(p,&q)==2 && isprime(2^q-1))
		) && #factor(6*k+1)[,1]==2 && sum(x=1,2,factor(6*k+1)[x,2])==2,
			print1(p",")
		);
		k=4*k
	)
};

Second you should split out your handling of the test of p from the rest of the program. This is just good general form. A first attempt would be a straight split:

m=vector(120,n,gcd(n-1,120)==1); \\ global

is(p,k)={
	my(q);
	(
		(m[(p%120)+1] && core(p)==p && isprime(p))
	||
		(issquare(p) && isprimepower(p,&q)==2 && isprime(2^q-1))
	) && #factor(6*k+1)[,1]==2 && sum(x=1,2,factor(6*k+1)[x,2])==2
};

try()={
	for(k=0,2^199,
		my(p=logint(6*k+2,2));
		if(is(p, k),
			print1(p",")
		);
		k=4*k
	)
};

Let's work on is() first. There's no need for m, it would be better to do the gcd each time. Much faster than core(n)==n is the function issquarefree(n), but you don't even need this at all because primes are always squarefree. That gives:

is(p,k)={
	my(q);
	(
		(gcd(p,120)==1 && isprime(p))
	||
		(issquare(p, &q) && isprime(q) && isprime(2^q-1))
	) && #factor(6*k+1)[,1]==2 && sum(x=1,2,factor(6*k+1)[x,2])==2
};

try()={
	for(k=0,2^199,
		my(p=logint(6*k+2,2));
		if(is(p, k),
			print1(p",")
		);
		k=4*k
	)
};

Next I notice that you're factoring 6k+1 three times, which seems extreme. One way to handle this would be to store the results of factoring, but in this case there's an easier solution: just look at the exponent vector. Numbers relatively prime to 120 are just those relatively prime to 30, but we can simplify further since that path is just looking for primes so we can just test that p > 5. In fact that test can be pulled outside (since the only prime square <= 5 is 4 which doesn't pass). Oh, and we need a definition of logint, so here we go:

logint(N,b)=my(k);while(N>=b,N\=b;k++);k;

is(p,k)={
	my(q);
	p>5 && 
	(isprime(p) || (issquare(p, &q) && isprime(q) && isprime(2^q-1)))
	&& factor(6*k+1)[,2]==[1,1]~
};

try()={
	for(k=0,2^199,
		my(p=logint(6*k+2,2));
		if(is(p, k),
			print1(p",")
		);
		k=4*k
	)
};

Now let's look at try(). I don't like abusing for() loops like this, a while() loop would make the transformation k -> 4k+1 clearer. But I think that it would be even clearer if you just wrote (4^n - 1)/3, and since the calculations inside is() take far longer than this calculation I'd just do that directly. (k <= 2^199 corresponds to n <= 100 here.) That also makes it clear that p = 2n+1 which was not at all clear before! Also we can remove our new function logint() since it isn't used.

is(p,k)={
	my(n);
	p>5 && 
	(isprime(p) || (issquare(p, &n) && isprime(n) && isprime(2^n-1)))
	&& factor(6*k+1)[,2]==[1,1]~
};

try()={
	for(n=0,100,
		my(k=(4^n-1)/3, p=2*n+1);
		if(is(p, k),
			print1(p",")
		)
	)
};

Now the code is much clearer, and you can see that you're just looking for semiprimes of the form 2^m - 1 subject to some conditions on m -- probably necessary conditions for semiprimality. But we can go further! You're ultimately testing 2^p - 1 for primality, and if p is the square of a prime q then the factorization of 2^p - 1 must be 2^q-1 times another prime, so we don't need to factor 2^p-1 in generality in that case. This also lets us clean up the code to remove k, which is no longer needed.

is(p)={
	my(q);
	if(p < 7, return(0));
	if(issquare(p, &q) && isprime(q) && isprime(2^q-1)),
		return(isprime((2^p-1)/(2^q-1)))
	);
	isprime(p) && factor(2^p-1)[,2]==[1,1]~
};

try()={
	for(n=0,100,
		my(p=2*n+1);
		if(is(p),
			print1(p",")
		)
	)
};

Evidently we can clean up the loop in try() to go directly over the odd numbers, no need for n. Also, since we're requiring it to be at least 7, let's just start the loop there and remove that condition from is().

is(p)={
	my(q);
	if(issquare(p, &q) && isprime(q) && isprime(2^q-1),
		return(isprime((2^p-1)/(2^q-1)))
	);
	isprime(p) && factor(2^p-1)[,2]==[1,1]~
};

try()={
	forstep(p=7,201,2,
		if(is(p), print1(p","))
	)
};

At this point I'm done simplifying, but you could improve this further. For example, you could replace the test isprime(2^q-1) with a function LL(q) which uses the Lucas-Lehmer test to see if 2^q-1 is prime, and you could write trial division code that checks for small prime factors of 2^p-1 (using only the ones which could be factors, of course) -- if you find one you don't have to complete the factorization, you can just test the cofactor for primality.

Edit: Another good step would be to replace try() with try(lim=201) and using forstep(p=7,lim,2 so you can test different regions without editing the function.

[science_man_88]

sometimes I think you know too much about PARI in the sense that my all inbuilt function using codes are half as fast as the ones you make. I was doing the other testing to decrease the number of checks but I guess my checks took so much time that the decrease didn't matter. I thought about using k as the exponent but could never get it to work. my code caught 9 but not 4 and yours appears to forget about 9. I'll see what I can implement of your further suggestions.

[CRGreathouse]

I get

9,11,23,37,41,49,59,67,83,97,101,103,109,131,137,139,149,167,197,199,

in 7 seconds. This could be improved, no doubt.

Edit: I did have to remove a superfluous paren from my code, try it now to see if it works.

[science_man_88]

Originally Posted by CRGreathouse

I get

9,11,23,37,41,49,59,67,83,97,101,103,109,131,137,139,149,167,197,199,

in 7 seconds. This could be improved, no doubt.

Edit: I did have to remove a superfluous paren from my code, try it now to see if it works.

(16:51) gp > is(p)=my(q);if(issquare(p,&q)&&isprime(q)&&isprime(2^q-1),(isprime((2^p-1)/(2^q-1))));isprime(p)&&factor(2^p-1)[,2]==[1,1]~;
(17:57) gp > try(lim=201)={forstep(p=7,lim,2,if(is(p), print1(p",")))};
(17:57) gp > try()
11,23,37,41,59,67,83,97,101,103,109,131,137,139,149,167,197,199,
(17:58) gp > ##
  ***   last result computed in 3,970 ms.

oh and I tried to make the LL check but it made it take longer probably because I suck I was thinking of just using your M code ( if I can find it) or a vector of known to be mersenne prime exponents because even with a complete list I could use the code in theory to check up to 924309391636849 or a little higher.

[CRGreathouse]

The ultimate optimization would be something like
2^p-1
2^(p^2)-1

[science_man_88]

Originally Posted by CRGreathouse

The ultimate optimization would be something like
2^p-1
2^(p^2)-1

I'm trying now to optimize to check if the division is prime from the formula I got for when the exponent's a certain multiple plus something and the reduced binomial theorem on wikipedia for when you have a variable plus 1 to a power I get that the result of division is:

my(n=q);sum(k=0,n-1,binomial(n,k+1)*y^k)

but I'm trying to think the best way to use this to test primality.