Is there utility in implementing an LLT-like sequence for GIMPS?

[kriesel]

Was looking at http://tony.reix.free.fr/Mersenne/index.html this morning, and saw his LLT-like primality test using S₀=5, S_i+1 = 2 S_i²-1 mod Mp proven, from 2005. http://tony.reix.free.fr/Mersenne/Pr...nneNumbers.pdf

It appears to me to be slightly slower to calculate than the standard LLT, S₀=4, S_i+1=s_i²-2 Mod Mp.
Decrement or double decrement are trivial, O(1), but the doubling in front seems to me to be O(n log n log log n) additional work (where n is the exponent) in an fft based implementation. A big multiword bit shift left by one. Smaller contribution to the run time constant in front, but of similar order.

Perhaps it would be useful, as a quick verification method of a new-found Mersenne prime.
Computing a different conclusive sequence on the fastest available software and hardware might be faster than using a different software and different hardware architecture combination.
Consider the cases where a new Mersenne prime is found with the usual LL sequence, with prime95/mprime or CUDALucas, or a PRP is found with gpuOwL or prime95/mprime and determined prime with CUDALucas.
Then what are the verification hardware and software candidates? Presumably Mlucas, and Glucas again, which can take a while. Perhaps Reix's sequence offers a faster verification alternative.

Would it be worthwhile to also include this different sequence in the standard prime95 software, to immediately launch a verification run on a new-found prime on the lucky owner's own system? That check run would inherently have a head start over any verification run by another user.
The different seed, and different sequence, may give different error behavior on the same hardware/OS/software combination, just as different offsets may affect error results. (It might be more palatable for the unlucky user that gets a false positive somehow, to have his own system break the news that it did not duplicate, than to wait days or perhaps weeks for someone else to complete a run and tell him. Such a self check if positive would not constitute sufficient proof of primality for GIMPS or EFF, but may constitute near proof of compositeness if negative.)

Given the current and future exponents in the tens or hundreds of millions, error checking is increasingly important.
Does the Jacobi check apply to this sequence?
Does some other more effective check apply?
Being a different sequence, it should hit the round off errors differently. Because of the factor of two, I'd expect it to require slightly different exponent cutoffs for the various fft lengths. If oriented to verification rather than first time runs, its fft limits could be set rather more conservatively than for the standard LL sequence.

Maybe this ground has all been thoroughly covered and there's nothing new or sufficiently useful to GIMPS here. But I don't recall seeing this sequence discussed. Especially not in the context of the recently added Jacobi and Gerbicz error checks. (If I just missed it, a URL pointer in the right direction would be welcome.)

If it's worthwhile, it would eventually have Primenet consequences; another work type to handle, and another possible form of discovery information leak to be addressed. Or the client could report the self check run locally only, I suppose.

[Batalov]

Originally Posted by kriesel

Was looking at http://tony.reix.free.fr/Mersenne/index.html this morning, and saw his LLT-like primality test using S₀=5, S_i+1 = 2 S_i²-1 mod Mp proven, from 2005.

It appears to me to be slightly slower same to calculate than the standard LLT, S₀=4, S_i+1=s_i²-2 Mod Mp.

S₀=5, S_i+1 = 2 S_i²-1 mod Mp
is equivalent to
S₀=10, S_i+1=S_i²-2 Mod Mp
...you are just tracking half-values with the former of the latter.

S₀=10 has no difference to the S₀=4, strategically. Of course it can be run, at all times, as is; there is a configuration setting for it.

[kriesel]

Originally Posted by Batalov

S₀=5, S_i+1 = 2 S_i²-1 mod Mp
is equivalent to
S₀=10, S_i+1=S_i²-2 Mod Mp
...you are just tracking half-values with the former of the latter.

S₀=10 has no difference to the S₀=4, strategically.

Are you saying Tony Reix's formula is effectively the same as LL with seed 10, offset -1? Following Tony Reix's algebra, the computation is a different form; 2 s_i²-1 mod Mp. And it's not just half value, necessarily. A few simple examples evaluated numerically demonstrate that.

It seems somewhat similar to me, to George's resort sometimes to computing (x-r)(x+r)+r² if x² hits a bad patch; using different algebra to accomplish the same end while avoiding numerical evaluation issues that can turn up.
For Reix, square, bit shift left with carries, decrement, is different algebra and code than for the usual LL, square, decrement by two, and the operands are different values for the same iteration generally.

Of course it can be run, at all times, as is; there is a configuration setting for it.

In what software, and where's that documented?
Result of search by any of various keywords, and =, for configuration settings compatible with the above statement:
prime95: readme.txt no
whatsnew.txt no
undoc.txt no
GUI menus no
prime.txt no
gwnum.txt no
mprime: same doc files as prime95
(there are some interesting controls for PRP, not relevant to LL)

CUDALucas 2.06: cudalucas.ini no
readme.txt no
source code nothing obvious

Glucas 2.9.0 doc no

Mlucas readme, no
source no, appears to exclude it.
Mlucas source snippet:
/* Always use 3 as the p-1 and Fermat-test seed, and 4 for the LL-test seed: */
if(MODULUS_TYPE == MODULUS_TYPE_MERSENNE && TEST_TYPE == TEST_TYPE_PRIMALITY)
iseed = 4;
else
iseed = 3;

gpuowl not applicable (PRP)

Being standardized on one seed value makes interim residues comparable. That's useful for comparisons between runs for where they begin to diverge in case of an error or more than one.

[GP2]

Originally Posted by kriesel

In what software, and where's that documented?

In mprime/Prime95, and it's not.

Use InitialLLValue=10 in prime.txt.

[Batalov]

Originally Posted by kriesel

Are you saying Tony Reix's formula is effectively the same as LL with seed 10?

Of course. Can't you see it?
Look:
First method, in Russian usually called, 'методом тыка' (literally, 'prodding' - if you want to know if something is soft or hard, just prod it with a stick):

$ gp -q
? s=5;for(i=1,5,s=2*s*s-1;print(2*s))
98
9602
92198402
8500545331353602
72259270930397519221389558374402
? s=10;for(i=1,5,s=s*s-2;print(s))
98
9602
92198402
8500545331353602
72259270930397519221389558374402

Second (proper): algebraically, do a substitution S=s*2 and observe that

S=10; repeat { S := S*S-2 }
is equivalent to

s*2=10; repeat { s*2 := s*2*s*2-2 } // now, divide by 2
-->
s=5; repeat { s := 2*s*s-1 }


___________________
Next: Re: your PDF. and your words: "And it's not just half value, necessarily. A few simple examples evaluated numerically demonstrate that."

Well, not if you don't know how to divide! You may want to learn that.

Let's start with page 1. Mp=8191
Your table says (you even put mod Mp in the beginning of each line, so you have to live by the mod Mp rules. Floating division in this context is quite silly.)

row 4: mod Mp (sic!) 1411 / 4801 = 0.293897105
No it's not!
mod Mp : 1411 / 4801 = 2 exactly

and of course it will be 2 for every row, e.g.
row 8: mod Mp 2113 / 5152 = 0.410131988 2

[Batalov]

Originally Posted by GP2

In mprime/Prime95, and it's not.

Use InitialLLValue=10 in prime.txt.

I think this is deliberately undocumented. If you leave an exposed configuration setting in the framework of the large distributed project, sooner or later some fraction of folks will set it. And then all the results from that computer will be useless to the GIMPS server, until this is reset.

The double-checker's computer also needs to be set to S₀ = 10, then with different shifts, these two results would form a valid data pair, in the context of the project (and the S₀ = 4, 10 or 2/3 should be recorded to the database). If the GIMPS server was interested in a randomized control trial of this kind, it could be sending some tasks with the seed of 10, but twice.

By the way, quick, easy homework for all interested:
calculate S₀ = 2/3 (Mod 2^p-1) without searching the forum.
Just take a sheet of paper and find out what it is. It is an "integer" modular value.

[kriesel]

Originally Posted by Batalov

Of course. Can't you see it?
Look:
First method, in Russian usually called, 'методом тыка' (literally, 'prodding' - if you want to know if something is soft or hard, just prod it with a stick):

$ gp -q
? s=5;for(i=1,5,s=2*s*s-1;print(2*s))
98
9602
92198402
8500545331353602
72259270930397519221389558374402
? s=10;for(i=1,5,s=s*s-2;print(s))
98
9602
92198402
8500545331353602
72259270930397519221389558374402

Second (proper): algebraically, do a substitution S=s*2 and observe that

S=10; repeat { S := S*S-2 }
is equivalent to

s*2=10; repeat { s*2 := s*2*s*2-2 } // now, divide by 2
-->
s=5; repeat { s := 2*s*s-1 }


___________________
Next: Re: your PDF. and your words: "And it's not just half value, necessarily. A few simple examples evaluated numerically demonstrate that."

Well, not if you don't know how to divide! You may want to learn that.

Let's start with page 1. Mp=8191
Your table says (you even put mod Mp in the beginning of each line, so you have to live by the mod Mp rules. Floating division in this context is quite silly.)

row 4: mod Mp (sic!) 1411 / 4801 = 0.293897105
No it's not!
mod Mp : 1411 / 4801 = 2 exactly

and of course it will be 2 for every row, e.g.
row 8: mod Mp 2113 / 5152 = 0.410131988 2

Well, sure, if you substitute a different term in place of what Reix wrote, you can make the sequence whatever you want. In this case an additional doubling each time, not present in his paper, which is a different S_i, different sequence, producing different values.

Re row 4 and 8: I'm looking at the actual magnitudes that need to be represented and processed by the algorithm. Perhaps it would have been better if I'd more fully expressed the rows as "Value after the mod Mp operation" instead of merely "mod Mp" but that takes up a lot of space. Or "remainder". You seem to be saying
1411(plus an imputed 8191) / 4801 = 2. That's true. (It seems a bit arbitrary to choose adding a single Mp when it's convenient to your point; magnitude ratios go all over the place for differing multiples of Mp added back, and there's practically speaking a disincentive to add back anything in performing the actual calculations.) And that seems to me irrelevant to the problem at hand, of ensuring correct execution; the values that must be represented and used for the next iteration and handled without error are not padded with multiples of Mp that have already been removed, at the mod Mp operation. The point of the ratio was to show that the relative magnitudes that must be handled per iteration fluctuate differently and do not between the two series have values always corresponding to a 2:1 ratio in magnitude. That seems to me a source for possibly different error behavior to occur. The algorithm does not have to be capable of squaring 9602 error-free, if the modulo reduction was complete. Now, to look at that more closely, I'd need to break the residues into parts for a multiword representation, but that is more than I'm inclined to do right now. Perhaps you'd be happier if I had used deltas rather than ratios. I used ratios of the actual values because you earlier made reference to a constant ratio. Which does not hold in all iterations. Equivalence is not equality.
1411 is a SMALLER magnitude than 4801 to represent in storage and process in the next iteration and survive the roundoff error and other errors. Its magnitude is not larger by a factor of two than 4801.


A separate matter is your impolite tone. George and Ernst are worth emulating in that respect. I don't understand why you seem hostile to someone trying to understand better. And who are you that you believe you can dictate the terms of discussion? Maybe it's a language thing, or a bad day, I don't know.

[ATH]

Originally Posted by Batalov

By the way, quick, easy homework for all interested:
calculate S₀ = 2/3 (Mod 2^p-1) without searching the forum.
Just take a sheet of paper and find out what it is. It is an "integer" modular value.

3 * (2*(2^p+1)/3-1) = 2*(2^p+1)-3 = 2*(2^p-1)+1 = 1 (Mod 2^p-1) =>
3^-1 = 2*(2^p+1)/3-1 = 2*W_p - 1 (Mod 2^p-1)

S₀ = 2/3 = 2*3^-1 = 2*(2*W_p - 1) = 4*W_p - 2 = 2^p+1 + W_p - 2 = W_p (Mod M_p)


Spoiler tag does not work with SUP and SUB tags for some reason.

[Batalov]

Noch eine Minute, bitte.

Originally Posted by kriesel

Well, sure, if you substitute a different term in place of what Reix wrote, you can make the sequence whatever you want.

I haven't changed anything in what you wrote. What Reix wrote is irrelevant. I was replying to what you wrote, wasn't I?
Maybe you quoted from him correctly, maybe you didn't.
You asked if S0=5, Si+1 = 2 Si^2-1 mod Mp was a new method, and you received an answer.

Originally Posted by kriesel

...(It seems a bit arbitrary to choose adding a single Mp when it's convenient to your point; magnitude ratios go all over the place for differing multiples of Mp added back, and there's practically speaking a disincentive to add back anything in performing the actual calculations.)
... I used ratios of the actual values because you earlier made reference to a constant ratio.

Who said anything about "to choose adding a single Mp when it's convenient to your point"?
Working mod Mp means you operate on mods:
Mod(1411,8191) / Mod(4801,8191) = 2 exactly
The exact way you can visualize it? - there are many shortcuts, including what you wrote (with one important change: No, not "when it's convenient to your point". When numerator is odd. That's all. ... Mod an odd number, a half always exists and this how it can be visualized: if numerator is even, divide by 2; else add modulo and divide by 2)

Originally Posted by kriesel

A separate matter is your impolite tone. George and Ernst are worth emulating in that respect. I don't understand why you seem hostile to someone trying to understand better. And who are you that you believe you can dictate the terms of discussion? Maybe it's a language thing, or a bad day, I don't know.

Impolite is one thing. Hostile is another.

A suggestion to learn about modular division was impolite? Ok, it was curt. It certainly wasn't hostile and was not out of place while discussing the modular algorithm that you suggested.
It helps you get to the answer and get to many more answers.
While, in contrast, by refusing to do so you will continue to run in darkness in circles and lose time.
Which is fine if that's what you want - nobody is forcing you. It was just a suggestion, and it was polite: "You may want to learn that".


"I'm here to help - if my help's not appreciated then lotsa luck, gentlemen." (W.Wolf)

[Batalov]

Originally Posted by ATH

Wp

Correct!

[kriesel]

Originally Posted by Batalov

I haven't changed anything in what you wrote. What Reix wrote is irrelevant. I was replying to what you wrote, wasn't I?
Maybe you quoted from him correctly, maybe you didn't.

What Reix wrote in the linked paper is central. It's what I was asking about at the outset; in the thread title and post one.

I provided a link to his paper in the beginning of post one, for the convenience of any reader. Here it is again, and, with red boxes added by me to highlight the most relevant initial portions, on the screen capture attached. http://tony.reix.free.fr/Mersenne/Pr...nneNumbers.pdf
Actually following what Reix wrote occurred to me as different computer code:
LL: square and subtract two mod MP
LL-Reix: square and double and subtract one mod MP
And there's precedent for using nonoptimal-compute-time code sequences as workarounds for numerical issues that arise at times.

You asked if S0=5, Si+1 = 2 Si^2-1 mod Mp was a new method, and you received an answer.

I asked if computing it had utility for GIMPS. Not if computing the usual LL with a different seed had utility for GIMPS, which is an interesting separate question, and very fortunate that that is already implemented somewhere.

Working mod Mp means you operate on mods:
Mod(1411,8191) / Mod(4801,8191) = 2 exactly

But I wasn't working mod Mp at that point, in a side calculation, comparing magnitudes of integers that are the result of operations and destined to be operands for the next computer operations. We switch numbers in or out of modulo math routinely. How many hours is it between noon today and 3pm tomorrow, for example; 27 hours, not 3. How many hours is it between noon and 3am? (same day, -9 hours, not 15 or 3.) What's the ratio between the time between 9am and noon, and 9am to 9pm, 1/4.

A part of the possible point of the Reix series is the numerical differences in interim computation values. Another part is the different iteration formula.

Kriesel: Are you saying Tony Reix's formula is effectively the same as LL with seed 10, offset -1?

(I was taught by a previous employer that the most effective way of ensuring accurate communication was to paraphrase back what was heard. This creates an opportunity for the other to check for accuracy or misunderstanding, iterating if necessary until details or nuance acceptably match.)

Batalov: Of course. Can't you see it?

Maybe it's all moot, of no further interest since the Reix sequence I initiated this thread about equivalent to another possibility that's already implemented in one software. And the talent and time to write, maintain, document, and tune software is the scarcest resource in the project.

The conclusion that base is already well covered is not surprising.

Does the Jacobi check apply to this sequence?
Does some other more effective check apply?

The Reix sequence being equivalent to LL seed 10 at offset -1, as shown by Batalov, the Jacobi check should apply, and Gerbicz not, and any as yet undiscovered error check for the LL sequence equally.

There's more, but



(because I couldn't find the "beating a dead horse" smily I've seen here before)