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)

 JeppeSN 2018-03-15 10:13

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

I was creating a new OEIS entry as a kind of procrastination...

I wonder what the first (probable) prime of the form [$]a^{16384}+b^{16384}[/$] is. Since this type of (extended) generalized Fermat numbers is mentioned in many sources (in the web), I would think someone had determined the answer? I had no lucky googling.

For each [$]m<14[/$], brute force will relatively early find a probable prime [$]a^{2^m}+b^{2^m}[/$]. The last of these is [$]72^{8192} + 43^{8192}[/$] which can be found i Caldwell's database: [URL="https://primes.utm.edu/primes/page.php?id=109871"]72^8192 + 43^8192[/URL]

So what about [$]m \ge 14[/$]?

/JeppeSN

 science_man_88 2018-03-15 10:37

[QUOTE=JeppeSN;482389]I was creating a new OEIS entry as a kind of procrastination...

I wonder what the first (probable) prime of the form [$]a^{16384}+b^{16384}[/$] is. Since this type of (extended) generalized Fermat numbers is mentioned in many sources (in the web), I would think someone had determined the answer? I had no lucky googling.

For each [$]m<14[/$], brute force will relatively early find a probable prime [$]a^{2^m}+b^{2^m}[/$]. The last of these is [$]72^{8192} + 43^{8192}[/$] which can be found i Caldwell's database: [URL="https://primes.utm.edu/primes/page.php?id=109871"]72^8192 + 43^8192[/URL]

So what about [$]m \ge 14[/$]?

/JeppeSN[/QUOTE]

a and b must be coprime. Their powers can't be addiive inverse remainders of each other.Etc.

 axn 2018-03-15 11:22

these are of the form a^16384+(a+1)^16384

 axn 2018-03-15 11:33

what is the form of factors for these numbers? It will be k*2^14+1, but will k itself has other restrictions (like k is a multiple of 2 or 4 or something? other modular restrictions?)

 Dr Sardonicus 2018-03-15 11:35

[QUOTE=JeppeSN;482389]I was creating a new OEIS entry as a kind of procrastination...

I wonder what the first (probable) prime of the form [$]a^{16384}+b^{16384}[/$] is. [/QUOTE]

I would say 2 (a = b = 1).

 JeppeSN 2018-03-15 12:12

[QUOTE=axn;482393]what is the form of factors for these numbers? It will be k*2^14+1, but will k itself has other restrictions (like k is a multiple of 2 or 4 or something? other modular restrictions?)[/QUOTE]

I would say the multiplier must be even, so it becomes (under a new convention where k is odd) [$]k\cdot 2^n + 1[/$] where [$]n \ge 15[/$]. See also [URL="http://www.prothsearch.com/GFNfacs.html"]New GFN factors (Wilfrid Keller)[/URL].

/JeppeSN

 JeppeSN 2018-03-15 12:19

[QUOTE=Dr Sardonicus;482395]I would say 2 (a = b = 1).[/QUOTE]

I should have said I was interested in odd primes only (and the number 1 is not prime).

Therefore a and b cannot both be odd. And as said by science_man_88 above, additionally a and b must be coprime.

These comments imply that a and b are distinct. Without loss of generality, a > b > 0.

/JeppeSN

 Dr Sardonicus 2018-03-15 13:23

[QUOTE=JeppeSN;482398]I should have said I was interested in odd primes only (and the number 1 is not prime).[/QUOTE]
There are some large "generalized Fermat" primes listed, e.g. [url=https://primes.utm.edu/top20/page.php?id=12]here[/url]; these have the forms

a^(2^20) + 1, a^(2^19) + 1, and a^(2^18) + 1.

Of course, they might not be the smallest prime a^(2^k) + b^(2^k) for their exponents.

In looking at smaller exponents, it did occur to me to look at cases

(a^(2^k) + b^(2^k))/2 with a*b odd. I noticed that (3^2 + 1)/2 and (3^4 + 1)/2 were primes, and, since I knew 2^32 + 1 isn't prime, I tried n = (3^32 + 1)/2. Pari's ispseudoprime(n) returned 1...

 VBCurtis 2018-03-15 15:33

[QUOTE=JeppeSN;482398]I should have said I was interested in odd primes only (and the number 1 is not prime).
/JeppeSN[/QUOTE]

Your OP never mentioned a and b should be prime. I don't see why "1 is not prime" is relevant. Do you mean to now require a>b>1?

 JeppeSN 2018-03-15 16:10

[QUOTE=VBCurtis;482416]Your OP never mentioned a and b should be prime. I don't see why "1 is not prime" is relevant. Do you mean to now require a>b>1?[/QUOTE]

No. I wanted to forbid [$]1 = 1^{16384} + 0^{16384}[/$] as another trivial "solution".

I do not want [I]any[/I] trivial solutions.

I know [$](a, b) = (67234, 1)[/$] gives an acceptable solution ([$]b=1[/$] is allowed), but it is certainly not minimal (although we do not have the proof until someone gives an example (EDIT: axn's first post in this thread already gave examples)).

Nitpicking is fine, but hopefully everyone sees I am just asking for the equivalent of [$]72^{8192} + 43^{8192}[/$] with exponents [$]16384[/$].

/JeppeSN

 JeppeSN 2018-03-15 16:46

Me: [QUOTE=JeppeSN;482389]For each [$]m<14[/$], brute force will relatively early find a probable prime [$]a^{2^m}+b^{2^m}[/$]. The last of these is [$]72^{8192} + 43^{8192}[/$] which can be found i Caldwell's database: [URL="https://primes.utm.edu/primes/page.php?id=109871"]72^8192 + 43^8192[/URL][/QUOTE]

To be explicit, this is what brute force finds:

[CODE]m=0, 2^1 + 1^1
m=1, 2^2 + 1^2
m=2, 2^4 + 1^4
m=3, 2^8 + 1^8
m=4, 2^16 + 1^16
m=5, 9^32 + 8^32
m=6, 11^64 + 8^64
m=7, 27^128 + 20^128
m=8, 14^256 + 5^256
m=9, 13^512 + 2^512
m=10, 47^1024 + 26^1024
m=11, 22^2048 + 3^2048
m=12, 53^4096 + 2^4096
m=13, 72^8192 + 43^8192[/CODE]

The next line takes more time. The question is: Hasn't this been considered before?!

/JeppeSN

 paulunderwood 2018-03-15 17:26

[QUOTE=JeppeSN;482422]

The next line takes more time. The question is: Hasn't this been considered before?!

/JeppeSN[/QUOTE]

Are you computing with PFGW :question:

Edit:

I ran PFGW for bases less than or equal to 50 and found no PRP:

[CODE]cat Jeppe.abc2
ABC2 $a^16384+$b^16384
a: from 2 to 50 step 2
b: from 3 to 50 step 2
[/CODE]
[CODE]./pfgw64 -f -N Jeppe.abc2[/CODE]

 science_man_88 2018-03-15 19:00

[QUOTE=JeppeSN;482422]Me:

To be explicit, this is what brute force finds:

[CODE]m=0, 2^1 + 1^1
m=1, 2^2 + 1^2
m=2, 2^4 + 1^4
m=3, 2^8 + 1^8
m=4, 2^16 + 1^16
m=5, 9^32 + 8^32
m=6, 11^64 + 8^64
m=7, 27^128 + 20^128
m=8, 14^256 + 5^256
m=9, 13^512 + 2^512
m=10, 47^1024 + 26^1024
m=11, 22^2048 + 3^2048
m=12, 53^4096 + 2^4096
m=13, 72^8192 + 43^8192[/CODE]

The next line takes more time. The question is: Hasn't this been considered before?!

/JeppeSN[/QUOTE]
Probably has, think modular tricks like fermats little theorem and extentions. Mod 8 it becomes a^0+b^0 for example, mod 9 it's a^4+b^4 these work for all bases coprime to 8 or 9 respectively.

 JeppeSN 2018-03-16 06:18

[QUOTE=paulunderwood;482426]Are you computing with PFGW :question:
[/QUOTE]
See [URL="https://oeis.org/A291944"]A291944 in OEIS[/URL]; it is not public yet, so see its history.

I used PARI/GP [FONT="Fixedsys"]ispseudoprime[/FONT] in a loop, like the code shown there, and I suspect Robert G. Wilson v used Mathematica. Maybe PFGW is faster?

There is no point in all of us running the same tests, except whoever uses the best tools will "win" the competition. I just thought maybe this had been established already.

/JeppeSN

 axn 2018-03-16 08:30

[QUOTE=JeppeSN;482473]There is no point in all of us running the same tests, except whoever uses the best tools will "win" the competition. I just thought maybe this had been established already.[/QUOTE]

I will try to write a custom sieve during the weekend. After that I can post the sieve output here, so that interested people can test ranges.

PFGW should indeed be faster than Pari or Mathematica.

EDIT:- Testing 71^16384+46^16384, PFGW took about 20s, while Pari took 2mins and change. So PFGW is about 6x faster.

 paulunderwood 2018-03-16 10:33

[QUOTE=axn;482479]I will try to write a custom sieve during the weekend. After that I can post the sieve output here, so that interested people can test ranges.

PFGW should indeed be faster than Pari or Mathematica.

EDIT:- Testing 71^16384+46^16384, PFGW took about 20s, while Pari took 2mins and change. So PFGW is about 6x faster.[/QUOTE]

PFGW will be more than 6x faster with much bigger numbers :smile:

 JeppeSN 2018-03-16 10:36

Something to note:

Use here the convention [$]a > b > 0[/$]. There is a slight chance that the smallest odd prime [$]a^{16384}+b^{16384}[/$] does not minimize [$]a[/$].

As an example, [$]677 < 678[/$], but still [$]677^{128}+670^{128} > 678^{128}+97^{128}[/$] (both of these sums of like powers are prime).

However, for the smallest one with that exponent, [$]27^{128}+20^{128}[/$], the value [$]a=27[/$] is also minimal. And I think this will be the case generally, because the bases [$]a[/$] and [$]b[/$] will be relatively small (I conjecture). But we will check for that with 16384 once axn's excellent initiative has come to fruition.

/JeppeSN

 ATH 2018-03-16 12:19

[QUOTE=axn;482479]I will try to write a custom sieve during the weekend. After that I can post the sieve output here, so that interested people can test ranges.

PFGW should indeed be faster than Pari or Mathematica.

EDIT:- Testing 71^16384+46^16384, PFGW took about 20s, while Pari took 2mins and change. So PFGW is about 6x faster.[/QUOTE]

I'm already working on a<b<=1000, or b<a<=1000 which ever convention you use

I used fbncsieve to sieve the factors k*2^14+1. It took only ~2min up to k=10^9.

Then I used these prime factors in a quickly written GMP program to sieve an array 1000x1000 of a,b. First I removed all values where b>=a, a<2, b<2, a%2=b%2 (both odd or both even), and gcd(a,b)>1. Down to 61K candidates at k=462M.

I'm running pfgw while continuing to trial factor. So far no PRP in 2<=b<=16 and b<a<=1000.

 axn 2018-03-16 13:15

[QUOTE=ATH;482490]I'm already working on a<b<=1000, or b<a<=1000 which ever convention you use

I used fbncsieve to sieve the factors k*2^14+1. It took only ~2min up to k=10^9.

Then I used these prime factors in a quickly written GMP program to sieve an array 1000x1000 of a,b. First I removed all values where b>=a, a<2, b<2, a%2=b%2 (both odd or both even), and gcd(a,b)>1. Down to 61K candidates at k=462M.[/quote]
Cool. But a custom sieve will be much more efficient. Hopefully that will be useful for >= 15.

[QUOTE=ATH;482490]I'm running pfgw while continuing to trial factor. So far no PRP in 2<=b<=16 and b<a<=1000.[/QUOTE]
Since the objective is to find the smallest, you should test in a different order.

2<=a<=1000, 1<=b<a

 a1call 2018-03-16 13:51

Stating the obvious for the sake of having it stated

a+b | a^q + b^q for all odd q

And

a+bi | a^q + b^q for all even q

So the result will be definitely not prime over the imaginary field.

Corrections are welcome.:smile:

 JeppeSN 2018-03-16 13:53

[QUOTE=axn;482493]Since the objective is to find the smallest, you should test in a different order.[/QUOTE]

I was about to say just that. If the (a,b) space is visualized as this triangle:
[CODE]
(2,1)
(3,2)
(4,1) (4,3)
(5,2) (5,4)
(6,1) (6,3) (6,5)
(7,2) (7,4) (7,6)
(8,1) (8,3) (8,5) (8,7)
. .
. .
. .
[/CODE]
it is best to search from the top and down, not from left to right.

/JeppeSN

 axn 2018-03-16 14:21

[QUOTE=axn;482493]But a custom sieve will be much more efficient.[/QUOTE]

That was somewhat presumptuous of me. What is the rate at which your sieve is progressing thru the factor candidates?

If you can post your sieve source, I can use that as a starting point.

 Batalov 2018-03-16 14:39

Were these the droids you were looking for?
[SPOILER]216^16384+109^16384[/SPOILER]

 JeppeSN 2018-03-16 15:08

[QUOTE=Batalov;482504]Were these the droids you were looking for?
[SPOILER]216^16384+109^16384[/SPOILER][/QUOTE]

Yes! And is it certain that 216^16384+109^16384 is minimal?

Did you search and find this now, or was it something that was known (to you) in advance?

Of course, now I want then minimal odd prime a^(2^m)+b^(2^m) for all the next m values :smile:

/JeppeSN

 paulunderwood 2018-03-16 16:30

[QUOTE=Batalov;482504]Were these the droids you were looking for?
[SPOILER]216^16384+109^16384[/SPOILER][/QUOTE]

For added confidence in Serge's result:

[CODE]time ./pfgw64 -k -f0 -od -q"[SPOILER]216^16384+109^16384[/SPOILER]" | ../../coding/gwnum/hybrid -

Testing (x + 1)^(n + 1) == 2 + 3 (mod n, x^2 - 3*x + 1)...
Likely prime!

real 1m45.793s
user 2m36.120s
sys 0m55.244s
[/CODE]

 Dr Sardonicus 2018-03-16 16:30

If one wanted to check exponents 2^m for several values of m, one might want to compile a table of pairs (a,b) with a < b, b <= N (given upper bound), gcd(a, b) = 1, [i]and[/i] a + b odd. This isn't entirely trivial. And I'm quite poor at this sort of thing.

Without the condition that a + b be odd, it would be equivalent to calculating the Farey sequence of order N. Perhaps the way to go is to calculate the Farey sequence, and then skip pairs where a and b are both odd.

One idea I considered by rejected was, for [i]even[/i] b, find the positive integers r < b with gcd(r, b) = 1, then fill in the entries (b, r + k*b) with r + k*b not exceeding N. Alas, you would have to make sure, for each r, that the multiplier k was relatively prime to r...

Of course, viewing table computation as a preliminary calculation, it probably makes little difference [i]how[/i] clumsily it is compiled, if N isn't very big. For large exponents, the value of b will be a good indicator of the size of a^(2^m) + b^(2^m)

The idea of sieving out multiples of "small" primes p congruent to 1 (mod 2^(m+1)) makes sense. Each pair (a, b) sieved out as giving a multiple of p means you can skip a PRP test. I looked at the primes p congruent to 1 (mod 2^15) up to the limit 50 million. There are 174 such primes. Grasping at straws(*), I computed the product (1 - 2^14/p) over these primes, and got .507+ which would be great if it meant eliminating almost half the pairs (a, b) actually under consideration.

(*) I reckoned that the relevant thing here was the 2^14 solutions to Mod(r,p)^(2^14) + 1 = Mod(0,p). Each pair (Mod(a,p), Mod(r, p)*Mod(a,p)) gives a multiple of p. This ignores the question of which such pairs might reduce to pairs of small coprime integers (a, b) of different parity. But hey -- any PRP test you can skip cheaply is to the good, right?

I did a quick check of the idea with the exponent 16. I looked at the smallest prime (p = 97) congruent to 1 (mod 32), and easily found from the root Mod(19,97) to

x^16 + 1 = Mod(0,97)

that 5^16 + 2^16 is a multiple of 97.

 Batalov 2018-03-16 16:52

[QUOTE=JeppeSN;482508]Yes! And is it certain that 216^16384+109^16384 is minimal?

Did you search and find this now, or was it something that was known (to you) in advance?

Of course, now I want then minimal odd prime a^(2^m)+b^(2^m) for all the next m values :smile:

/JeppeSN[/QUOTE]
It is minimal.
What puzzles me is why is this taking others so long to find the same? It only took overnight on 8 cores. The workflow is very simple: generate (a,b)s in 1 second with a trivial Pari script (opposite parity, gcd > 1, say, up to a<=400), then split into a few workfiles (hawever many cores you have) and run [c]pfgw -f"{65536}" -N -k -l input[/c] on each chunk, periodically check, ..., done!

For m=15, it is perhaps best to sieve, but for m=14, this -f parameter takes care of the small factors just well enough.

 ATH 2018-03-16 17:15

It is just a question of resources. For this random question I only used 1 core for sieving and 1 core for pfgw, I had no need to use more, rest is running Prime95 on this machine, it is not like we had to find the answer within a certain time limit.

I only tried briefly using pfgw for trial factoring long ago, and my impression was it was horribly slow at it, so I have not tried again, I might try next time I need it.

Edit: and yes I should have sorted the list so it ran pfgw on candidates by size, but I made this quick program last night before I went to bed, and I only got 4-5 hours of sleep before work.

 axn 2018-03-16 17:41

[QUOTE=Batalov;482516]run [c]pfgw -f"{65536}" -N -k -l input[/c] on each chunk, periodically check, ..., done![/QUOTE]
The should be f{32768}. You will miss a prolific factor 10*16384+1

 Batalov 2018-03-16 18:28

A typo, this command-line is already prepared for m=15. :-)

 science_man_88 2018-03-16 18:49

Technically if a and b are squares you find one for the next level up as well.

 CRGreathouse 2018-03-16 19:36

[QUOTE=science_man_88;482533]Technically if a and b are squares you find one for the next level up as well.[/QUOTE]

Yes, but 216 and 109 aren't (though the former is a cube).

 Batalov 2018-03-16 20:00

Another droid crossed my path. [SPOILER]260^32768+179^32768[/SPOILER]
Not proven minimal yet, though.

 paulunderwood 2018-03-16 20:18

[QUOTE=Batalov;482538]Another droid crossed my path. [SPOILER]260^32768+179^32768[/SPOILER]
Not proven minimal yet, though.[/QUOTE]

Firming up this result:

[CODE]time ./pfgw64 -k -f0 -od -q"[SPOILER]260^32768+179^32768[/SPOILER]" | ../../coding/gwnum/hybrid -

Testing (x + 2)^(n + 1) == 7 (mod n, x^2 - x + 1)...
Likely prime!

real 3m34.502s
user 5m52.580s
sys 1m13.776s
[/CODE]

 science_man_88 2018-03-16 20:28

[QUOTE=CRGreathouse;482537]Yes, but 216 and 109 aren't (though the former is a cube).[/QUOTE]

m=0, 2^1 + 1^1
m=1, 2^2 + 1^2
m=2, 2^4 + 1^4
m=3, 2^8 + 1^8
m=4, 2^16 + 1^16
m=5, 9^32 + 8^32
m=6, 11^64 + 8^64
m=7, 27^128 + 20^128
m=8, 14^256 + 5^256
m=9, 13^512 + 2^512
m=10, 47^1024 + 26^1024
m=11, 22^2048 + 3^2048
m=12, 53^4096 + 2^4096
m=13, 72^8192 + 43^8192

Aka

(72^2)^4096+(43^2)^4096
(53^2)^2048+(2^2)^2048
(22^2)^1024+(3^2)^1024
(47^2)^512+(26^2)^512
(13^2)^256+(2^2)^256
(14^2)^128+(5^2)^128
(27^2)^64+(20^2)^64
(11^2)^32+(8^2)^32
(9^2)^16+(8^2)^16
(2^2)^8+(1^2)^8
(2^2)^4+(1^2)^4
(2^2)^2+(1^2)^2
(2^2)^1+(1^2)^1

Minimal next level= minimal square based current level. Aka Serges previous result has a>14 or b>10 or both for the next level.

 Batalov 2018-03-16 23:00

[QUOTE=Batalov;482538]Another droid crossed my path. [SPOILER]260^32768+179^32768[/SPOILER]
[STRIKE]Not proven minimal yet, though[/STRIKE].[/QUOTE]
Now proven minimal.

m=16, anyone?

 science_man_88 2018-03-17 00:20

[QUOTE=Batalov;482562]Now proven minimal.

m=16, anyone?[/QUOTE]

Not unless we find or can easily apply more because with opposite parity and gcd(x,y)=1, my calculations had under 9000 pairs up to (209,208) but more than 18200 up to (300,299)

 Batalov 2018-03-17 01:06

[QUOTE=science_man_88;482567]my calculations had under 9000 pairs up to (209,208) but more than 18200 up to (300,299)[/QUOTE]
Yeah, my calculations also show that 209/300 ~= 0.7 is a damn good approximation of the 1/ sqrt2. Consequently, but of course, there are almost exactly 2 times more pairs for a<=300 than for a<=209. So what? The number of pairs with a<=n will be O(n[SUP]2[/SUP]) (in reality, slightly less: the gcd criterion will help).

And under a<=600 there are ~4 times more pairs than under a<=300. So what, just keep sieving and testing! What seems to be the problem?

It is excellent that the number of tries goes quadratic -- you go almost flat in size (and in test runtime); and one likely need to test γ*2^m pairs, i.e. ~30,000 to 40,000. Just keep sieving and testing!

 science_man_88 2018-03-17 01:59

[QUOTE=Batalov;482571]Yeah, my calculations also show that 209/300 ~= 0.7 is a damn good approximation of the 1/ sqrt2. Consequently, but of course, there are almost exactly 2 times more pairs for a<=300 than for a<=209. So what? The number of pairs with a<=n will be O(n[SUP]2[/SUP]) (in reality, slightly less: the gcd criterion will help).

And under a<=600 there are ~4 times more pairs than under a<=300. So what, just keep sieving and testing! What seems to be the problem?

It is excellent that the number of tries goes quadratic -- you go almost flat in size (and in test runtime); and one likely need to test γ*2^m pairs, i.e. ~30,000 to 40,000. Just keep sieving and testing![/QUOTE]

I might see if I can figure out my additive inverse property to help. It may even be relatable to the arithmetic mean of the pair.

 Batalov 2018-03-17 02:39

Aw well, m=16 was too easy. Even a pair for the good measure:
124^65536+57^65536
143^65536+106^65536

 axn 2018-03-17 02:58

[QUOTE=Batalov;482574]Aw well, m=16 was too easy. Even a pair for the good measure:
124^65536+57^65536
143^65536+106^65536[/QUOTE]

Are all of these being reported to TOP PRP site? I know you're not reporting them to factordb.

 paulunderwood 2018-03-17 03:08

[QUOTE=Batalov;482574]Aw well, m=16 was too easy. Even a pair for the good measure:
124^65536+57^65536
143^65536+106^65536[/QUOTE]

Tests for Serge's latest finds:
[CODE]time ./pfgw64 -k -f0 -od -q"124^65536+57^65536" | ../../coding/gwnum/hybrid -

Testing (x + 1)^(n + 1) == 2 + 3 (mod n, x^2 - 3*x + 1)...
Likely prime!

real 8m40.139s
user 15m28.568s
sys 1m48.060s[/CODE]

[CODE]time ./pfgw64 -k -f0 -od -q"143^65536+106^65536" | ../../coding/gwnum/hybrid -

Testing (x + 2)^(n + 1) == 7 (mod n, x^2 - x + 1)...
Likely prime!

real 8m29.679s
user 15m1.912s
sys 1m47.804s[/CODE]

 Dr Sardonicus 2018-03-17 03:45

[QUOTE=Batalov;482571]Yeah, my calculations also show that 209/300 ~= 0.7 is a damn good approximation of the 1/ sqrt2. Consequently, but of course, there are almost exactly 2 times more pairs for a<=300 than for a<=209. So what? The number of pairs with a<=n will be O(n[SUP]2[/SUP]) (in reality, slightly less: the gcd criterion will help).[/QUOTE]
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
;-)

 Batalov 2018-03-17 04:28

[QUOTE=axn;482577]Are all of these being reported to TOP PRP site? I know you're not reporting them to factordb.[/QUOTE]
Sure did.

 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.

 Batalov 2018-03-18 20:51

[QUOTE=a1call;482727]There is not much saving if any, in doing more than one such test.[/QUOTE]
That's right.
You will spend time better doing trial factoring mod k*2[SUP]m+1[/SUP]+1.
You do realize that gcd() is nothing other than a bunch of divisions?

 a1call 2018-03-18 20:57

I meant per candidate. Trial by division is very efficient for filtering out candidates with small factors. GCD of the form I suggest can probably knock out a handful of candidates with much larger factors. Any little bit should help when you have very expensive PRP testing to do.

 science_man_88 2018-03-18 21:40

[QUOTE=a1call;482730]I meant per candidate. Trial by division is very efficient for filtering out candidates with small factors. GCD of the form I suggest can probably knock out a handful of candidates with much larger factors. Any little bit should help when you have very expensive PRP testing to do.[/QUOTE]

But k and k' can't be arbitrary, if a and b are opposite parity then a+k has same parity as b+k' if both pairs are opposite parity pairs. If they are both same parity pairs then the whole expression value is even ( as is the previous case as well). So really with the parity check already in place we have done this for mod 2. And the gcd knocks out other factors. It's because additive inverses can't easily be done that we can't really do that much with it.

 a1call 2018-03-18 21:51

You need to generate testing pairs which are Mod (1,q) including their odd factors.
Correct me if I wrong but the ks are arbitrary.

For the same of GCD the generated partner does not need to be odd, just appropriately large. Not t0o large to contradict the time saved.
Now I am out to attend the Noruz party.
Will check in later
ETA I think if your for steps generate say odd a and even b then the ks will have to be both odd. The point is they could be any positive or negative odd addends.

 Dr Sardonicus 2018-03-19 14:25

For the sieving part, you need, presumably, the primes p == 1 (mod 2*e), e = 2^m, up to some limit. My stick-pokes at the problem indicate that the limit will be much, [i]much[/i] larger than that for precomputed tables of primes. So, you have to roll your own. This would be another initial computation. The question is, what's a reasonable limit for a given exponent e = 2^m?

I note that, if you proceed from exponent 2^m to exponent 2^(m+1), you can cull the list of primes p == 1 (mod 2^(m+1)) for at least the first part of your list of primes p == 1 (mod 2^(m+2)); but I'm guessing that the larger the exponent e = 2^m, the higher your limit on primes will be.

Another thing occurs to me: The numbers a^e + b^e being sieved here are, as noted in the OP to this thread, divisible [i]only[/i] by primes congruent to 1 (mod 2*e). The number of such numbers less than or equal to X is, asymptotically[sup]*[/sup], c*X*log(X)^(1/e - 1); I'm too lazy to work out the value of c. Anyhow, this might affect how well sieving is likely to work.

[sup]*[/sup]Bateman, Paul T. and Diamond, Harold G. [u]Analytic Number Theory[/u] [i]An Introductory Course[/i] Theorem 10.1.

The intrepid web surfer might try [url=https://www.worldscientific.com/worldscibooks/10.1142/5605]here[/url].

 science_man_88 2018-03-19 14:52

[QUOTE=a1call;482736]
ETA I think if your for steps generate say odd a and even b then the ks will have to be both odd. The point is they could be any positive or negative odd addends.[/QUOTE]

No, because of mod 3 etc. if you have a and b, both 0 mod 3, then both k's can't be without leading to a number the GCD we already do would eliminate. If a and b are additive inverses mod 3 adding k that are additive inverse pairs causing additive inverse pairs also cause a number gcd would eliminate already.

 Batalov 2018-03-19 14:55

Set n = limit of search for a; e = 2^m.
[CODE]Initialize an array of all eligible pairs (a,b), 0 < b < a <= n, gcd(a,b)=1.
Allocate array uint64_t s[n].

For each prime p == 1 (mod 2*e) up to 63 bits {
For each small prime q < n { s[q] = q^e by powmod }
then all composite s[b] as a prod of s[q]
for each surviving pair of (a,b): remove pair from list if s[a]+s[b] == p
}[/CODE]Up to 63 bits, the above code will work - even without gmp.
There are multiple existing cuda sieves that could be adapted to the above code.

 Dr Sardonicus 2018-03-20 02:23

[QUOTE=Batalov;482789]Set n = limit of search for a; e = 2^m.
[snip]
For each prime p == 1 (mod 2*e) up to 63 bits
[snip]
[/QUOTE]
Hmm, use primes up to 63 bits. I guess "sieve as much as possible" is indeed the way to go. Alas, my computing resources are limited.

Meanwhile, out of curiosity, I looked at pairs (a, b) with a < b, gcd(a, b) = 1, and a and b both [i]odd[/i]. Of course, a^e + b^e is even, but with e = 2^m, (a^e + b^e)/2 is odd. The following are the smallest a and b I found for which (a^e +b^e)/2 is a pseudoprime.

With e = 2^3 and e = 2^7, I got an unexpected bonus! At e = 2^12, my Pari script seems to have hit the wall. It will either grind through my set of pairs (a,b) without result, or it will find a pseudoprime. Either way, I'm done with this variant.

e = 2^1: a = 1, b = 3

e = 2^2: a = 1, b = 3

e = 2^3: a = 1, b = 9

e = 2^4: a = 1 , b = 3

e = 2^5: a = 1, b = 3

e = 2^6: a = 1, b = 3

e = 2^7: a = 9, b = 49

e = 2^8: a = 3, b = 7

e = 2^9: a = 9, b = 35

e = 2^10: a = 57, b = 67

e = 2^11: a = 49, b = 75

 Batalov 2018-03-20 02:43

The odd b GFNs were kept on a close backburner for many years already:
[URL]http://www.prothsearch.com/GFN03.html[/URL]
[URL]http://www.prothsearch.com/GFN05.html[/URL]
[URL]http://www.prothsearch.com/GFN07.html[/URL]
[URL]http://www.prothsearch.com/GFN11.html[/URL]
All of them are (b^2^n+1)/2. There are some primes there (for b<12) and [URL="http://www.primenumbers.net/prptop/searchform.php?form=%28b%5En%2B1%29%2F2&action=Search"]tons of PRPs[/URL], of course.

If you replace b by b/a (with gcd(a.b)=1), the modular divisibility is not changing and intrinsic factor 2 is sticking around. ...There is [URL="http://www.prothsearch.com/GFNsmall.html"]a 217-digit prime[/URL] (7^2^8+3^2^8)/2 between small xGFNs.

 Batalov 2018-03-20 02:52

1 Attachment(s)
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. )

 Dr Sardonicus 2018-03-20 13:03

[QUOTE=Batalov;482844][snip] If you replace b by b/a (with gcd(a.b)=1), the modular divisibility is not changing and intrinsic factor 2 is sticking around. ...There is [URL="http://www.prothsearch.com/GFNsmall.html"]a 217-digit prime[/URL] (7^2^8+3^2^8)/2 between small xGFNs.[/QUOTE]

Thanks for the URLs. That b/a trick is nice for the "homogeneous" GFN's, saves a powmod every time while sieving.

Yesirree, (3^256 + 7^256)/2 is indeed prime. I'd checked it using isprime(), since it wasn't so awfully big. Took many times longer than ispseudoprime(), but only seconds, not minutes.

EDIT: I have now checked (9^512 + 35^512)/2 with isprime(), which took minutes. Also prime.

If someone wants to check (57^1024 + 67^1024)/2 and (49^2048 + 75^2048)/2, you're welcome to it. The manual has convinced me that isprime() would take a [i]very[/i] long time trying to deal with these.

 Dr Sardonicus 2018-03-20 16:24

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 [b]Batalov[/b]'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 2018-03-21 14:26

Once upon a time, long long ago -- [url=http://www.mersenneforum.org/showpost.php?p=482722&postcount=52]here[/url] in fact -- I posted the following ineptly-written script, along with its running time:

[quote]? 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.[/quote]

It suddenly occurred to me -- [i][b]D'OH![/b][/i] -- 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 [i]said[/i] 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 2018-03-22 16:52

The entry [URL="https://oeis.org/A291944"]A291944[/URL] is public now. /JeppeSN

 Batalov 2018-03-22 18:15

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

 CRGreathouse 2018-03-22 20:50

[QUOTE=JeppeSN;483072]The entry [URL="https://oeis.org/A291944"]A291944[/URL] is public now. /JeppeSN[/QUOTE]

Very nice!

 Dr Sardonicus 2018-03-23 15:34

[QUOTE=JeppeSN;483072]The entry [URL="https://oeis.org/A291944"]A291944[/URL] is public now. /JeppeSN[/QUOTE]
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 [i]definitely[/i] 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 2018-03-23 17:42

[QUOTE=Dr Sardonicus;483175] (If someone knows any decent heuristic estimates, though, please post them, or a link).[/QUOTE]
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 [FONT=&quot]~γ *[/FONT] 2[SUP]m[/SUP] pairs, and you have ~0.2 * a^2 pairs up to a, hence:
a(m) ~ 1.15 * 2[SUP]m/2[/SUP]
Plot the values and you will find great concordance with this estimate.

 axn 2018-03-24 04:15

1 Attachment(s)
[QUOTE=Batalov;482845]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. )[/QUOTE]

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 2018-03-24 08:22

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

 JeppeSN 2018-03-24 13:48

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

[CODE]
(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[/CODE]

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 2018-03-24 13:53

[QUOTE=JeppeSN;483252]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}[/$$][/QUOTE]. Not really, it happens any time both bases are squares.

 JeppeSN 2018-03-24 14:19

[QUOTE=science_man_88;483253]. Not really, it happens any time both bases are squares.[/QUOTE]

Of course, but I conjecture this is the only instance where, for a given [$]n[/$], the [I]smallest[/I] prime of form [$]F_n'(a,b)=\frac{a^{2^n}+b^{2^n}}{2}[/$] has [$]a[/$] and [$]b[/$] both perfect squares. /JeppeSN

 Dr Sardonicus 2018-03-24 14:49

[QUOTE=JeppeSN;483252]For both a and b odd, the first primes are (as others already found but did not post):

[CODE]
[snip]
(5^8+3^8)/2
[snip]
(67^1024+57^1024)/2
[/CODE]
[/QUOTE]

Sheesh, I don't know HOW I messed that one up. I had posted [url=http://www.mersenneforum.org/showpost.php?p=482841&postcount=63]here[/url] (1^8 + 9^8)/2. But (3^8 + 5^8)/2 is correct.

I also posted, for e = 2^11, n = (49^e + 75^e)/2.

 Batalov 2018-03-24 17:04

[QUOTE=JeppeSN;483260]Of course, but I conjecture this is the only instance where, for a given [$]n[/$], the [I]smallest[/I] prime of form [$]F_n'(a,b)=\frac{a^{2^n}+b^{2^n}}{2}[/$] has [$]a[/$] and [$]b[/$] both perfect squares. /JeppeSN[/QUOTE]
Not necessarily. It is the other way around - if the hit for the F[SUB]n+1[/SUB](a,b) happens very early, then it will decide the fate of the minimum in F[SUB]n[/SUB](a[SUP]2[/SUP],b[SUP]2[/SUP]).

Here, in sequence [URL]https://oeis.org/A246119[/URL], the same happened three time already.

 Dr Sardonicus 2018-03-24 17:37

[QUOTE=Batalov;483280]Not necessarily. It is the other way around - if the hit for the F[SUB]n+1[/SUB](a,b) happens very early, then it will decide the fate of the minimum in F[SUB]n[/SUB](a[SUP]2[/SUP],b[SUP]2[/SUP]).

Here, in sequence [URL]https://oeis.org/A246119[/URL], the same happened three time already.[/QUOTE]

Squares are rare, so if either a or b is a square it is worthy of note. And, lo and behold, for m = 11 (though my results should be checked), one of the numbers (49, 75) is a square. So I wouldn't dismiss the possibility of another instance of [i]both[/i] being squares out of hand.

Alas, the number of cases any of us will ever know of is, likely, quite small.

 science_man_88 2018-03-24 18:08

[QUOTE=Dr Sardonicus;483283]Squares are rare, so if either a or b is a square it is worthy of note. And, lo and behold, for m = 11 (though my results should be checked), one of the numbers (49, 75) is a square. So I wouldn't dismiss the possibility of another instance of [i]both[/i] being squares out of hand.

Alas, the number of cases any of us will ever know of is, likely, quite small.[/QUOTE]

They are more common, than the solutions to the next stage up... I think you mean minimal ones are of note.

 JeppeSN 2018-03-24 21:31

[QUOTE=Dr Sardonicus;483269]Sheesh, I don't know HOW I messed that one up. I had posted [url=http://www.mersenneforum.org/showpost.php?p=482841&postcount=63]here[/url] (1^8 + 9^8)/2. But (3^8 + 5^8)/2 is correct.

I also posted, for e = 2^11, n = (49^e + 75^e)/2.[/QUOTE]

Apparently I forgot your post. Sorry. I agree on exponent 2^11.

In the search for the next one, I find:

[CODE](157^(2^12)+83^(2^12))/2 is 3-PRP! (4.1038s+1.9638s)
(157^(2^12)+111^(2^12))/2 is 3-PRP! (4.0938s+3.3192s)
(161^(2^12)+157^(2^12))/2 is 3-PRP! (4.1053s+1.9778s)[/CODE]

/JeppeSN

 JeppeSN 2018-03-24 21:42

[QUOTE=Batalov;483280]Not necessarily. It is the other way around - if the hit for the F[SUB]n+1[/SUB](a,b) happens very early, then it will decide the fate of the minimum in F[SUB]n[/SUB](a[SUP]2[/SUP],b[SUP]2[/SUP]).

Here, in sequence [URL]https://oeis.org/A246119[/URL], the same happened three time already.[/QUOTE]

Yes, also interesting. I could be wrong, of course. With [I]two[/I] "free parameters" [$]a,b[/$] in [$]F_n'(a,b)[/$] is it likely to occur again for [$]n>11[/$]? I do not think so.

The form in A246119, [$]k^{2^n} (k^{2^n}-1)+1[/$], has only one parameter [$]k[/$]. I still think this square phenomenon should occur only finitely many times?

/JeppeSN

 JeppeSN 2018-03-24 21:56

[QUOTE=JeppeSN;483297]Apparently I forgot your post. Sorry. I agree on exponent 2^11.

In the search for the next one, I find:

[CODE](157^(2^12)+83^(2^12))/2 is 3-PRP! (4.1038s+1.9638s)
(157^(2^12)+111^(2^12))/2 is 3-PRP! (4.0938s+3.3192s)
(161^(2^12)+157^(2^12))/2 is 3-PRP! (4.1053s+1.9778s)[/CODE]

/JeppeSN[/QUOTE]

[B]EDIT[/B] The below numbers are not primes. PFGW picks an "unlucky" base:

And then:

[CODE](3^(2^13)+1^(2^13))/2 is 3-PRP! (0.6454s+0.4075s)[/CODE]

EDIT: And again:

[CODE](3^(2^14)+1^(2^14))/2 is 3-PRP! (2.8285s+1.3930s)[/CODE]

And:

[CODE](9^(2^15)+1^(2^15))/2 is 3-PRP! (57.8865s+21.0036s)[/CODE]

So easy(?):

[CODE](3^(2^16)+1^(2^16))/2 is 3-PRP! (57.7878s+24.4478s)[/CODE]

/JeppeSN

 Batalov 2018-03-24 23:03

Use base 2, obviously

 JeppeSN 2018-03-25 06:38

[QUOTE=JeppeSN;483297]Apparently I forgot your post. Sorry. I agree on exponent 2^11.

In the search for the next one, I find:

[CODE](157^(2^12)+83^(2^12))/2 is 3-PRP! (4.1038s+1.9638s)
(157^(2^12)+111^(2^12))/2 is 3-PRP! (4.0938s+3.3192s)
(161^(2^12)+157^(2^12))/2 is 3-PRP! (4.1053s+1.9778s)[/CODE]

/JeppeSN[/QUOTE]

For exponent 2^13, we have:

[CODE](107^(2^13)+69^(2^13))/2[/CODE]

/JeppeSN

 Dr Sardonicus 2018-03-25 16:00

[QUOTE=JeppeSN;483302][B]EDIT[/B] The below numbers are not primes. PFGW picks an "unlucky" base:

[code](3^(2^13)+1^(2^13))/2 is 3-PRP! (0.6454s+0.4075s)[/CODE]

[CODE](3^(2^14)+1^(2^14))/2 is 3-PRP! (2.8285s+1.3930s)[/CODE]

[CODE](9^(2^15)+1^(2^15))/2 is 3-PRP! (57.8865s+21.0036s)[/CODE]

[CODE](3^(2^16)+1^(2^16))/2 is 3-PRP! (57.7878s+24.4478s)[/CODE]
[/QUOTE]
Hardly surprising. We have

3^2 == 1 (mod 2^3). Therefore,

3^(2^k) == 1 (mod 2^(k+2)) for every positive integer k.

Now, let k be a positive integer, and N = (3^(2^k) + 1)/2. Then N - 1 = (3^(2^k) - 1)/2. By the above, we have

2^(k+1) divides N - 1, so

3^(2^(k+1)) - 1 divides 3^(N - 1) - 1. Since

3^(2^(k+1)) - 1 = (3^(2^k) - 1)*(3^(2^k) + 1), and (3^(2^k) + 1)/2 = N, we have

N divides 3^(N-1) - 1. That is, N = 3^(2^k) + 1 is a base-3 pseudoprime for every positive integer k.

I haven't checked the Rabin-Miller criterion in this case, though.

(The present instance is reminiscent of the fact that, if p is prime, 2^p - 1 is a base-2 pseudoprime, regardless of whether it's prime or composite.)

 JeppeSN 2018-03-25 20:46

I found one for 2^14 and now have:

[CODE]
(3^(2^0)+1^(2^0))/2
(3^(2^1)+1^(2^1))/2
(3^(2^2)+1^(2^2))/2
(5^(2^3)+3^(2^3))/2
(3^(2^4)+1^(2^4))/2
(3^(2^5)+1^(2^5))/2
(3^(2^6)+1^(2^6))/2
(49^(2^7)+9^(2^7))/2
(7^(2^8)+3^(2^8))/2
(35^(2^9)+9^(2^9))/2
(67^(2^10)+57^(2^10))/2
(49^(2^11)+75^(2^11))/2
(157^(2^12)+83^(2^12))/2
(107^(2^13)+69^(2^13))/2
(71^(2^14)+1^(2^14))/2[/CODE]

Would it be interesting to add this sequence, ..., 35, 67, 49, 157, 107, 71, ..., to OEIS?

/JeppeSN

 Batalov 2018-03-25 21:17

Yes, why not.
Any sequence is a sequence.
Most of the text can be reused from the previous sequence, and use this sister sequence for the upper boundary [URL]https://oeis.org/A275530[/URL]

 JeppeSN 2018-03-26 09:15

[QUOTE=Batalov;483379]Yes, why not.
Any sequence is a sequence.
Most of the text can be reused from the previous sequence, and use this sister sequence for the upper boundary [URL]https://oeis.org/A275530[/URL][/QUOTE]

To appear as [URL="https://oeis.org/A301738"]A301738[/URL] (until it is approved, see History to see the proposed versions). /JeppeSN

 JeppeSN 2018-03-26 09:44

We can also restrict ourselves to [I]consecutive[/I] odd bases:
[$$]\frac{a^{2^n}+(a-2)^{2^n}}{2}[/$$]
Can also be parametrized in other ways, such as the [$]k[/$] in:
[$$]\frac{(2k+1)^{2^n}+(2k-1)^{2^n}}{2}[/$$]
OEIS does not seem to have it either (searching for a, or for a-2, or for 2k=a-1, or for k). /JeppeSN

 axn 2018-03-26 11:02

[QUOTE=JeppeSN;483412]We can also restrict ourselves to [I]consecutive[/I] odd bases:
[$$]\frac{a^{2^n}+(a-2)^{2^n}}{2}[/$$]
Can also be parametrized in other ways, such as the [$]k[/$] in:
[$$]\frac{(2k+1)^{2^n}+(2k-1)^{2^n}}{2}[/$$]
OEIS does not seem to have it either (searching for a, or for a-2, or for 2k=a-1, or for k). /JeppeSN[/QUOTE]

[CODE](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
(179^128+177^128)/2
(169^256+167^256)/2
(935^512+933^512)/2
(663^1024+661^1024)/2[/CODE]

 axn 2018-03-26 13:23

1 Attachment(s)
So, I sieved n=17 for 0 < b < a <= 2048 till 2^53. Not sure what range or what depth Serge has sieved on this n, but it looks "sieved enough". This could be used to divide up PRP work.

[QUOTE=Batalov;483240]Neat!
...Now time for a quick cuda sieve facelift, based on gfnsieve? ^_^[/QUOTE]

The logic doesn't exactly correspond to gfn/cyclo sieves, so a cuda sieve will have to be built from scratch, more or less; I have to think about it a bit.

 Dr Sardonicus 2018-03-26 16:26

[QUOTE=JeppeSN;483375]Would it be interesting to add this sequence, ..., 35, 67, 49, 157, 107, 71, ..., to OEIS?
[/QUOTE]

No. The sequence itself is merely a curiosity. The real interest of this thread, at least for me, has been in the methods described for determining pairs (a, b) for which

a^(2^m) + b^(2^m), [or (a^(2^m) + b^(2^m))/2, if a and b are both odd]

is a (pseudo)prime in a reasonable length of time.

 Batalov 2018-03-26 21:28

axn's sieve is easily modified for odd-odd pairs. In fact nearly nothing needs to be changed (only the intake parity filters, if they are there; they were in mine, just drop them).

 axn 2018-03-27 03:02

[QUOTE=Batalov;483480]axn's sieve is easily modified for odd-odd pairs. In fact nearly nothing needs to be changed (only the intake parity filters, if they are there; they were in mine, just drop them).[/QUOTE]

No. The hash matching also needs to change (it currently puts even indexed residues in the hash and uses odd-indexed residues to probe the hash).

All times are UTC. The time now is 12:25.