 mersenneforum.org (https://www.mersenneforum.org/index.php)
-   Math (https://www.mersenneforum.org/forumdisplay.php?f=8)
-   -   Smallest prime of the form a^2^m + b^2^m, m>=14 (https://www.mersenneforum.org/showthread.php?t=23160)

 science_man_88 2018-03-17 11:44

[QUOTE=Dr Sardonicus;482596]The number of pairs (a, b) with 1 <= a <= b <= N, and gcd(a, b) = 1, is SUM, b = 1 to N, eulerphi(b). This is the number of fractions in the Farey sequence of order N.

Asymptotically, this is 6/pi^2 * N^2. Long known.

The additional condition that a + b is odd, knocks out the pairs (a, b) with a and b both odd. In this case, b - a is even and gcd(b - a, b) = 1, so for odd b > 1, the condition knocks out exactly 1/2*eulerphi(b) pairs. For b = 1 it knocks out the pair (1, 1).

I'm too lazy to work out the modified asymptotic. I'm sure it is long known, but alas I'm also too lazy to look it up
;-)[/QUOTE]

Using me and serges calculations it seems to be roughly 0.2*N^2 I get a>11 or b>10 or both for the next one.

 Dr Sardonicus 2018-03-17 15:14

[QUOTE=Dr Sardonicus;482596]The number of pairs (a, b) with 1 <= a <= b <= N, and gcd(a, b) = 1, is SUM, b = 1 to N, eulerphi(b). This is the number of fractions in the Farey sequence of order N.

Asymptotically, this is 6/pi^2 * N^2. Long known.

The additional condition that a + b is odd, knocks out the pairs (a, b) with a and b both odd. In this case, b - a is even and gcd(b - a, b) = 1, so for odd b > 1, the condition knocks out exactly 1/2*eulerphi(b) pairs. For b = 1 it knocks out the pair (1, 1).

I'm too lazy to work out the modified asymptotic. I'm sure it is long known, but alas I'm also too lazy to look it up
;-)[/QUOTE]

I screwed up. Mea culpa.

The coefficient 6/pi^2 is for the number of pairs of positive integer (a, b), both less than or equal to N, with gcd(a, b) = 1 -- but [i]without[/i] the requirement a < b. This coefficient may be viewed as the limiting value, as N increases without bound, of the probability that two positive integers, chosen at random in [1, N], are relatively prime. By considering separately the possible prime common factors 2, 3, 5, etc and a bit of handwaving, you obtain the coefficient as being the product

(1 - 1/2^2)*(1 - 1/3^2)*(1 - 1/5^2)...

whose inverse is obviously the sum, from k = 1 to infinity, of 1/k^2, whose value is well known to be pi^2/6. Thus, the coefficient is 6/pi^2.

Imposing the requirement a < b obviously cuts the number of pairs in half, giving the correct asymptotic for the Farey sequence of order N as 3/pi^2 * N^2. The coefficient 3/pi^2 is

0.30396355092701331433163838962918291671, approximately.

Being too lazy either to work out or look up the variant with the further requirement that a and b be of different parity, I simply had Pari-GP find the actual value numerically up to some ridiculous limit, and print out the value every so often. It was interesting to watch as an increasing number of decimal digits remained fixed, but the convergence was slow, so I Ctrl-C'd it out of my misery.

The asymptotic for the variant clearly is 2/pi^2 * N^2, the coefficient 2/pi^2 being

0.20264236728467554288775892641945527781 approximately.

 science_man_88 2018-03-17 17:00

[QUOTE=Dr Sardonicus;482619]I screwed up. Mea culpa.

The coefficient 6/pi^2 is for the number of pairs of positive integer (a, b), both less than or equal to N, with gcd(a, b) = 1 -- but [i]without[/i] the requirement a < b. This coefficient may be viewed as the limiting value, as N increases without bound, of the probability that two positive integers, chosen at random in [1, N], are relatively prime. By considering separately the possible prime common factors 2, 3, 5, etc and a bit of handwaving, you obtain the coefficient as being the product

(1 - 1/2^2)*(1 - 1/3^2)*(1 - 1/5^2)...

whose inverse is obviously the sum, from k = 1 to infinity, of 1/k^2, whose value is well known to be pi^2/6. Thus, the coefficient is 6/pi^2.

Imposing the requirement a < b obviously cuts the number of pairs in half, giving the correct asymptotic for the Farey sequence of order N as 3/pi^2 * N^2. The coefficient 3/pi^2 is

0.30396355092701331433163838962918291671, approximately.

Being too lazy either to work out or look up the variant with the further requirement that a and b be of different parity, I simply had Pari-GP find the actual value numerically up to some ridiculous limit, and print out the value every so often. It was interesting to watch as an increasing number of decimal digits remained fixed, but the convergence was slow, so I Ctrl-C'd it out of my misery.

The asymptotic for the variant clearly is 2/pi^2 * N^2, the coefficient 2/pi^2 being

0.20264236728467554288775892641945527781 approximately.[/QUOTE]

To be fair, my "calculation" was done with PARI/gp code that counted them for me. I still wish I could figure my additive inverse version out. If it weren't for the 14,15 cases it seemed most often a composite m had a odd and b even. Of course so did 7( and 8 breaks the general composite rule). Just looking at the data.

 Dr Sardonicus 2018-03-18 15:55

Just for practice, I (tried to) write Pari-GP scripts to perform various tasks associated with this problem. I compiled a vector of all pairs [a, b] with b up to 200 for which

a < b, gcd(a, b) = 1, and a + b is odd.

I deliberately left the testing of the conditions appallingly clumsy. Since I didn't bother to compute in advance how many pairs there were going to be, I used concat() to append each new pair to my vector. That took [i]way[/i] more time than computing the pairs [a, b]. No matter, I had my vector of pairs, 8156 of them. If I had done the computation of the number of pairs in advance, this part of the operation would have been much faster.

Next, I decided to see how effective sieving might be. Using the original exponent e = 2^14, I compiled a vector of primes p < 2^26, p congruent to 1 (mod 2^15). There are 228 of them.

I then wrote a script that initialized a count at 0 and went through my vector of pairs [a,b]. If my script is writ right (I invite anyone to check the results, I'm not good at writing scripts), for each pair [a,b] it goes through my vector of primes p, and computes

r = Mod(a,p)^e + Mod(b,p)^e.

Then if r is zero, it increments the count, and immediately moves on to the next pair [a,b].

(I note that the smallest value 2^(2^14) + 1 of a^e + b^e is [i]way[/i] bigger than my prime limit 2^26.)

At the end of the run, the count of pairs sieved out was 4065. Now 4065./8156 = 0.4984, approximately. The "grasping at straws" product of (1. - 2^14/p), taken over the 228 primes on my list, is 0.4995, approximately. That is much closer agreement than I had hoped for!

 science_man_88 2018-03-18 16:26

my code for counting them was just:

[CODE]my(a=0);for(x=1,200,forstep(y=x-1,1,-2,if(gcd(x,y)==1,a++)));a[/CODE]

Done in about 50 ms I guess. Where as the version using concat was 19591 ms I may be able to do it with two vectors with fewer concats to the main vector though.

 Batalov 2018-03-18 17:12

Just in case, if anyone is sieiving or testing for m=17, a(17) > 238.

 paulunderwood 2018-03-18 17:20

[QUOTE=Batalov;482708]Just in case, if anyone is sieiving or testing for m=17, a(17) > 238.[/QUOTE]

That's 311503 decimal digits -- PRP'ing with generic mod reduction is really going to bite :smile:

 Dr Sardonicus 2018-03-18 19:47

[QUOTE=science_man_88;482706]my code for counting them was just:

[CODE]my(a=0);for(x=1,200,forstep(y=x-1,1,-2,if(gcd(x,y)==1,a++)));a[/CODE]

Done in about 50 ms I guess. Where as the version using concat was 19591 ms I may be able to do it with two vectors with fewer concats to the main vector though.[/QUOTE]

I went back and did the counting in advance (up to the limit b = 200) thusly:

[code]? s=0;r=1;for(i=2,200,if(r,t=eulerphi(i),t=eulerphi(i)/2);s+=t;r=!r);print(s)
8156
time = 0 ms.[/code]

For creating the vector of pairs [a, b] I then did this: I'm sure the testing of whether i+j is odd can be handled more expeditiously, but it probably wouldn't make much difference in the time.

[code]? k=1;v=vector(s);for(j=2,200,for(i=1,j-1,if((i+j)%2==1&&gcd(i,j)==1,v[k]=[i,j];k++)))
time = 47 ms.
? print(v[s])
[199, 200][/code]

It looks to me like improving times in looking for pseudoprimes would depend on doing as much sieving as possible. This would mean extending the "factor base" of primes congruent to 1 (mod 2*e) way, [i]way[/i] beyond my tiny limit of 2^26. Doubling the [i]log[/i] of the prime limit would perhaps sieve out half the pairs left by sieving to the old limit...

 a1call 2018-03-18 20:04

My 2¢ worth,

The fastest filtering is size based since the objective is to find the smallest PRP.
You can set a narrow test range based on size, assign each band to a different thread and advance upwards after each completion. The size filtering before any trial by division or PRP test seems most efficient to me.
It should avoid a lot of unnecessary lager PRP tests.

 science_man_88 2018-03-18 20:19

[QUOTE=Dr Sardonicus;482722]I went back and did the counting in advance (up to the limit b = 200) thusly:

[code]? s=0;r=1;for(i=2,200,if(r,t=eulerphi(i),t=eulerphi(i)/2);s+=t;r=!r);print(s)
8156
time = 0 ms.[/code]

For creating the vector of pairs [a, b] I then did this: I'm sure the testing of whether i+j is odd can be handled more expeditiously, but it probably wouldn't make much difference in the time.

[code]? k=1;v=vector(s);for(j=2,200,for(i=1,j-1,if((i+j)%2==1&&gcd(i,j)==1,v[k]=[i,j];k++)))
time = 47 ms.
? print(v[s])
[199, 200][/code]

It looks to me like improving times in looking for pseudoprimes would depend on doing as much sieving as possible. This would mean extending the "factor base" of primes congruent to 1 (mod 2*e) way, [i]way[/i] beyond my tiny limit of 2^26. Doubling the [i]log[/i] of the prime limit would perhaps sieve out half the pairs left by sieving to the old limit...[/QUOTE]
You can also treat composite and prime a values seperately. Odd primes always have all opposite parity less than them, composites don't always except maybe powers of two. But again doubtful this decreases the time. EDIT: you actually have that all odd numbers have the powers of two less than them. So really the ones in need of findng, are only when 4n+2 and odd numbers crash.

 a1call 2018-03-18 20:47

Chances are you can eliminate some 10's of very expensive PRP tests by performing a GCD test between each candidate :

GCD (a^q+b^q,(a+k)^+(b+k')^q)

For arbitrary values of k and k'.

Such a test can speed things up, but gets less and less effective as the exponent/PRP gets larger.
There is not much saving if any, in doing more than one such test.

All times are UTC. The time now is 17:14.