Smallest prime of the form a^2^m + b^2^m, m>=14

[Dr Sardonicus]

I looked at the exponent e = 2^14. Rather than compute the gcd's of pairs of numbers a^e + b^e -- which would have been a godawful amount of computation -- I did the following: I used my list of primes p == 1 (mod 2^15) up to 2^26, and went through all the pairs [a,b] with a < b, a+b odd, gcd(a,b) = 1, and b < 201. For each prime p, I simply counted the number of pairs for which a^e + b^e was divisible by p (employing Batalov's b/a trick). This took seconds, not minutes.

What I found was, the first prime on my list, p = 65537, divided over 2,000 of the numbers a^e + b^e -- nearly a quarter of them. The number of a^e + b^e divisible by p decreased rapidly with p, but IIRC was at least 2 for the first 120 primes on the list, at least 3 for the first 103 primes on the list, and at least 4 for the first 63 primes on the list.

[Dr Sardonicus]

Once upon a time, long long ago -- here in fact -- I posted the following ineptly-written script, along with its running time:

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

It suddenly occurred to me -- D'OH! -- that I could insure that i + j was odd in more or less the same way I had kept track of parity in my script computing the number s of pairs (i,j) in advance. (Hey, I said I'm not very good at writing scripts!) This really helps, because it halves the number of steps in the inner loop.

? k=1;v=vector(s);f=0;for(j=2,200,forstep(i=f+1,j-1,2,if(gcd(i,j)==1,v[k]=[i,j];k++));f=!f)
time = 23 ms.

Half the steps, half the time...

For e = 2^14 I also looked at how many p == 1 (mod 2*e), p < 2^26, divide the numbers a^e + b^e. The champion pair was (116,127); 116^e + 127^e is divisible by 65537, 1146881, 1179649, 17367041, and 46497793.

[JeppeSN]

The entry A291944 is public now. /JeppeSN

[Batalov]

a(17) > 412
and sieved to a<=600 up to 52 bits (and higher in some ranges)

[CRGreathouse]

Originally Posted by JeppeSN

The entry A291944 is public now. /JeppeSN

Very nice!

[Dr Sardonicus]

Originally Posted by JeppeSN

The entry A291944 is public now. /JeppeSN

That's pretty good. Since x^(2^N) + y^(2^N) is a "non-full form" for N > 1, there isn't much theory to guide estimates of the size of the minimal prime of that form, etc. So, numerical computation seems to be the way to go. (If someone knows any decent heuristic estimates, though, please post them, or a link).

It occurred to me that, for a given N, odd (pseudo)primes of the form a^(2^N) + b^(2^N) are of two types: smaller number is odd, and smaller number is even. I computed the smallest of both types up to the limit N = 10, without doing any sieving, because it would have taken me longer to set up the sieving, than it took the clunky test-all-pairs scripts to run. (Pressing on to N = 11 and higher, though, it would definitely be worth the effort of doing a LOT of sieving.)

My default assumption is that the minimal [a, b] are of each type (approximately) equally often as the upper bound for N increases. I have no idea how much larger the minimal pair of the other type may be compared to the minimal pair; whether there may be several pairs of the same type as the minimal pair which smaller than the minimal pair of the other type; etc.

N = 1: [1, 2], [2, 3]
N = 2; [1, 2], [2, 3]
N = 3; [1, 2], [2, 13]
N = 4; [1, 2], [6, 7]
N = 5; [8, 9], [13, 18]
N = 6; [8, 11], [7, 16]
N = 7; [20, 27], [31, 32]
N = 8; [5, 14], [34, 43]
N = 9; [2, 13], [35, 38]
N = 10; [26, 47], [39, 56]
N = 11; [3, 22], ??
N = 12; [2, 53], ??
N = 13; [43, 72], ??
N = 14; [109, 216], ??
N = 15; [179, 260], ??
N = 16; [57,124], ??

Finally, the "two types" issue doesn't occur for the (pseudo)primes (a^(2^N) + b^(2^N))/2 with a and b both odd.

[Batalov]

Originally Posted by Dr Sardonicus

(If someone knows any decent heuristic estimates, though, please post them, or a link).

Clearly there is no precise estimate, just like with any Poisson process. The estimate with 70%confidence interval will be +- 100% of the trend value. (and the 95% CI will be from 0 to several times trend.)

Trend value's estimate is very easy. Nearly all candidates will have an equal chance to be prime which only depends on their size which is almost exactly the same in bits. So, one will have to try on the order of ~γ * 2^m pairs, and you have ~0.2 * a^2 pairs up to a, hence:
a(m) ~ 1.15 * 2^m/2
Plot the values and you will find great concordance with this estimate.

[axn]

Originally Posted by Batalov

Here is an unpolished sieve (and a lazy one - I didin't implement the pseudocode above fully; that is, - I am simply exponentiating all a's). It is really ugly but it works. I am sieving m=17 case to 50-51 bits with this.

( PFGW -f{...} essentially sieves to 44.5-45 bits de facto, so this is a bit of an improvement above simply prefactoring with PFGW. )

So I got bored, and modified your source to add the missing bits. And some bells and whistles.

This one uses explicit k values in cmdline (instead of bits). Also reads/writes files instead of stdin/stdout. Reads sieve.txt for candidates, appends factor.txt for factors, and at the end of the range, writes sieve.txt back with survivors. It also self starts -- if you give starting k as 1, it will generate the initial sieve itself instead of reading from file (but in this case, a third parameter maxa needs to be specified to define the sieve region).

The algorithm is O(maxa), and so it is best to just sieve as large a region as you want from the start (rather than trying to increase it later).

EDIT:- Right now, hard coded to N=17. Change the NN #define, and you're set.

[Batalov]

Neat!
...Now time for a quick cuda sieve facelift, based on gfnsieve? ^_^

[JeppeSN]

For both a and b odd, the first primes are (as others already found but did not post):

(3^1+1^1)/2
(3^2+1^2)/2
(3^4+1^4)/2
(5^8+3^8)/2
(3^16+1^16)/2
(3^32+1^32)/2
(3^64+1^64)/2
(49^128+9^128)/2
(7^256+3^256)/2
(35^512+9^512)/2
(67^1024+57^1024)/2

It is interesting that the exponent 128 gives a number that is identical to that for 256, i.e. \[\frac{49^{128}+9^{128}}{2}=\frac{7^{256}+3^{256}}{2}\]

[science_man_88]

Originally Posted by JeppeSN

It is interesting that the exponent 128 gives a number that is identical to that for 256, i.e. \[\frac{49^{128}+9^{128}}{2}=\frac{7^{256}+3^{256}}{2}\]

. Not really, it happens any time both bases are squares.