Double-Double Arithmetic

[R. Gerbicz]

More solutions:
depth=29;x=5638349967449211460;p=7049705358222163969;
depth=29;x=4669797239304715507;p=7128230911014862849;
depth=29;x=1068131362738116904;p=7301107626761256961;
depth=29;x=5914744651162318648;p=7612253723611365377;
depth=29;x=4114583113667178426;p=7931690348588302337;
depth=29;x=3191062263392789539;p=8585113616301162497;

and there is no more up to p=2^63 (the useful range),
searched for the p=k*2^29+1 form. Note that here the listed "depth" is limited at f, where 2^f||(p-1),
because in some cases this can be arbitrarily large, small example:
2^(2^(2*t))==2 mod 7 for every t integer.

Interestingly for these primes there is a single Google hit, also a mersenneforum(!)
topic: http://www.mersenneforum.org/showthread.php?t=1081
I see no missed solution for my search.

[ewmayer]

Originally Posted by Mysticial

Yes, there will be 2^28 different deep roots. Given any root-of-two of depth k, you can generate another one of the same depth by multiplying it by a root-of-unity of the same depth or less.

I prefer to work things the other way around - first find a primitive root of maximal order (call said ring order 'ord'), which seems to require only a handful of trials with various small candidates that pass the simplest filter (power-mod to raise to power = ord/2, pass if result == -1), to ensure they also yield non-unity when raised to ord divided by each of its odd-prime factors, then power one's primitive root as needed to get any desired root of unity in the tree.

Which brings up a new point that I hadn't realized before. If there exists a root-of-two of depth k, then there must also exist a root-of-unity of depth k.

Perhaps I'm missing something obvious, but how do you figure?

Originally Posted by Mysticial

I believe the answer is yes. But I'm becoming less and less sure the more you ask since I'm in no way an expert in the field.

The reason why I think the answer is yes is because once you find any root-of-two (R) of depth (say 2^28), you will automatically be able to generate the remaining (2^28-1) others of the same depth by multiplying them by roots-of-unity. (assuming those exist at the same depth)

You seem to be mixing your roots up here - you multiply roots of 2 by each other to get other roots of 2, roots of 1 by each other to get other roots of 1.

These 2^28 roots-of-two must also be unique. Since p is prime and a multiplicative inverse always exists, you will be able to divide them by R to re-obtain the roots-of-unity.

The roots of 2 which I demonstrated for FGT(M(p)) are not unique - they are only as many distinct ones of them as there are powers of 2 (mod 2^p-1), i.e. p distinct ones, one of which is unity - and yet they are of full order > p.

[Mysticial]

Originally Posted by R. Gerbicz

More solutions:
depth=29;x=5638349967449211460;p=7049705358222163969;
depth=29;x=4669797239304715507;p=7128230911014862849;
depth=29;x=1068131362738116904;p=7301107626761256961;
depth=29;x=5914744651162318648;p=7612253723611365377;
depth=29;x=4114583113667178426;p=7931690348588302337;
depth=29;x=3191062263392789539;p=8585113616301162497;

and there is no more up to p=2^63 (the useful range),
searched for the p=k*2^29+1 form. Note that here the listed "depth" is limited at f, where 2^f||(p-1),
because in some cases this can be arbitrarily large, small example:
2^(2^(2*t))==2 mod 7 for every t integer.

Interestingly for these primes there is a single Google hit, also a mersenneforum(!)
topic: http://www.mersenneforum.org/showthread.php?t=1081
I see no missed solution for my search.

Interesting. Looks like that search was already done over a decade ago!

Originally Posted by ewmayer

You seem to be mixing your roots up here - you multiply roots of 2 by each other to get other roots of 2, roots of 1 by each other to get other roots of 1.

Wait, how am I wrong?

p = 10358097392821075969

Here's a root-of-two of depth 2^28.
36411718325718249 ^ (2^28) = 2 mod p

Here's a root-of-unity.
3175727014930767360^ (2^27) = -1 mod p
3175727014930767360^ (2^28) = 1 mod p

Multiply them together to get a different root-of-two of the same depth.
36411718325718249 * 3175727014930767360 mod p = 9888747190045455352
9888747190045455352 ^ (2^28) = 2 mod p

Using this method, I can generate 2^28 distinct roots-of-two of the same depth by using each of the different 2^28 different roots-of-unity of depth 2^28 or less. (This is somewhat tangent to the IBDWT since you only need one root. But I was mainly answering ATH's question about finding other roots of the same order.)

The roots of 2 which I demonstrated for FGT(M(p)) are not unique - they are only as many distinct ones of them as there are powers of 2 (mod 2^p-1), i.e. p distinct ones, one of which is unity - and yet they are of full order > p.

Yes. I realized that once Robert gave me an example where 2 itself is part of the cycle generated by squaring.

And you should be able to generate a lot more of them using the method above. But they'll be complex roots. Again, not particularly useful for IBDWT itself, but they do exist.

[ewmayer]

Yes, it came to me over dinner just now that if r and w are nontrivial roots of 2 and 1, respectively, each of order k, then (rw)^k = r^k.w^k = 2*1 = 2.

Which begs the question: Are the real power-of-2 roots of 2 which I found for FGT(M(p)) not in fact of full order, or do we simply have the choice there between real or complex roots of 2?

If not of full order, then that raises the issue of whether length-N IBDWT needs N distinct roots of 2 or not. My long-ago Fortran90 proof-of-principle FGT implementation indicated not, since it was a soup-to-nuts program (albeit restricted to power-of-2 transform lengths), and the resulting LL-test residues precisely matched those of my pure-floating-point-FFT version.

Another indication to me that IBDWT may not need distinct roots of 2 is from my Fermat-mod code path, which (unlike the codes used at time of the F24 primality test) supports non-power-of-2 transform lengths, in which case we combine both the standard negacyclic-transform-effecting weighting using (2N)th complex roots of unity and a Mersenne-mod-style variable-base-plus-fractional-roots-of-2-weights IBDWT. But for an FFT length of form odd.2^k said powers of 2 repeat every (odd)th term, i.e. there are only (odd) distinct ones of them.

[Mysticial]

Originally Posted by ewmayer

Which begs the question: Are the real power-of-2 roots of 2 which I found for FGT(M(p)) not in fact of full order, or do we simply have the choice there between real or complex roots of 2?

Perhaps the word "full order" isn't the right term to use. Unlike for roots-of-unity, each "order" for roots-of-two are not necessarily inclusive of all the lower orders. And it looks like the elements in each order will be mutually exclusive of all other orders iff 2 is not a root-of-two of itself.

This is not the case for Mersenne primes. And for M31: 2, 4, 16, 256, and 65536 together cover all orders to infinity. But the # of elements in each order is the # of distinct roots-of-unity of that order.

So there are 2^32 different roots-of-two of order 2^123456789 for M31 over the complex field.

[R. Gerbicz]

Just to complete the search to 2^64, the solutions for p in (2^63,2^64):

depth=29;x=7836412484265765590;p=9657578508943097857;
depth=29;x=5246133613423626286;p=9857256047176056833;
depth=29;x=8147572229873106544;p=9983467653877989377;
depth=28;x=1349217705146881478;p=10358097392821075969;
depth=28;x=3907128162577576875;p=12548406757803687937;
depth=28;x=12913619447375217277;p=14251526830518435841;
depth=28;x=10397636786069575256;p=14358224048763174913;
depth=30;x=13143085125098926212;p=14380913871361671169;
depth=29;x=13817898618308695410;p=15613990314799267841;
depth=29;x=15589368141424458430;p=15910894009111805953;
depth=30;x=535050801089448450;p=17061046145649737729;
depth=28;x=2163236610507737378;p=18195224058105167873;

There was a print of the solution if depth>=28 (but searched only the p=k*2^29+1 form).

These were known, but with wrong depth:
p=3348317433*2^32+1=14380913871361671169 has depth=30 and not 28
p=3635415415*2^32+1=15613990314799267841 has depth=29 and not 28.

Basically it is a pretty trivial problem, used only one thread, but with multithread and a built in sieve to find the primes in the k*2^29+1 sequence would solve the problem well in under one hour. Ofcourse used c, not the superslow Mathematica.

You can even use only https://en.wikipedia.org/wiki/Tonell...anks_algorithm to find x=x0, where x^2==2 mod p, then for this x0 solve x^2==x0 mod p etc., here it is not interesting if you choose x0 or (-x0), and the same is true for each depth.

A roughly 3-4 times faster algorithm, what used: it is known that x^e==b mod p is solvable iff b^((p-1)/gcd(p-1,e))==1 mod p (or b==0 mod p trivial case). If you want a high depth say x^(2^28)==2 mod p, then use the above condition for e=2^28;b=2; (obviously if the true depth is higher then this condition still holds).

[preda]

Originally Posted by Mysticial

A 31-bit prime works perfectly with Victor Shoup's NTT algorithm. I have an implementation of this (albeit real-numbers instead of Gaussian). And it's memory-bound on my 7900X with AVX512 - though the implementation isn't memory-optimal.

Please, I'd like to know what is "Shoup's NTT algorithm". Specifically, it appears there's a "Shoup's trick" (optimization) -- what is it? Asking to make sure I don't miss it.

[Mysticial]

Originally Posted by preda

Please, I'd like to know what is "Shoup's NTT algorithm". Specifically, it appears there's a "Shoup's trick" (optimization) -- what is it? Asking to make sure I don't miss it.

It's the same thing. I think I'm the only person who uses the term "Shoup's algorithm" when in reality, it's "just" a method (albeit a very power one) to perform a multiply-mod.

To compute: c = a * b mod p

B = the word size (such as 2^32 or 2^64)

a ∈ [0, B)
b ∈ [0, p)
p ∈ [1, B/2)  //  Can be any number. Doesn't need to be prime.

//  Precompute:
b' = floor(b * B / p)

//  Shoup's multiply-mod
q = floor(a * b' / B);   //  B x B integer multiply keeping the top half.
t0 = (a * b) mod B;      //  B x B integer multiply keeping the bottom half.
t1 = (q * p) mod B;      //  B x B integer multiply keeping the bottom half.
c = (t0 - t1) mod B;     //  Integer subtract ignoring overflow.

//  Output
c ≡ a * b mod p
c ∈ [0, 2p)

So the output is either fully reduced, or reduced + p. To get it back into the range: c ∈ [0, p) just do a compare and conditional-move.

Operation Count:

• 1 upper-half integer multply
• 2 lower-half integer multiplies
• 1 integer subtract
• 1 integer compare
• 1 conditional move

The one catch is that you need to be able to precompute b' to be efficient. This is possible for NTTs since b will be a twiddle factor. Otherwise, b' can be computed on the fly with a variant of this trick. But that makes the entire multiply-mod operation about 3x more expensive.

[ATH]

No new depth 29 in p=k*2^28+1 between 2^62 and 2^63, only:
depth 28, p=4998875103153356801
depth 28, p=6249316677578653697

These and most of those R. Gerbicz mentioned has "infinite" depth root-of-two, where you can keep taking the square root until you eventually get back to 2 and start over on the same values. Is that useful at all? or only of the root-of-unity is deep as well?
What is the condition of the prime for this to occur? I thought at first p=7 (mod 8) but that is not correct. Is there anyway to predict the length of the loop until it repeats? It must be at most ord_p(2) but often is it much less.

[R. Gerbicz]

Originally Posted by ATH

No new depth 29 in p=k*2^28+1 between 2^62 and 2^63, only:
depth 28, p=4998875103153356801
depth 28, p=6249316677578653697

These and most of those R. Gerbicz mentioned has "infinite" depth root-of-two, where you can keep taking the square root until you eventually get back to 2 and start over on the same values. Is that useful at all? or only of the root-of-unity is deep as well?
What is the condition of the prime for this to occur? I thought at first p=7 (mod 8) but that is not correct. Is there anyway to predict the length of the loop until it repeats? It must be at most ord_p(2) but often is it much less.

See, there are two cases: the depth<f or depth="infinity", where 2^f||p-1 (this notation means that 2^f is an exact primepower divisor, so 2^f divides, but 2^(f+1) doesn't divide p-1). You can see this quickly, because in the condition (p-1)/gcd(p-1,2^e) is the same if e>=f.

This also means that you have no chance to find depth=29 in the k*2^28+1 sequence, because the k=even case explored by me, and for the odd case depth<28 or infinite.

[ATH]

When I did 2^62-2^63 for k*2^28+1 I also searched even k's, so I verified all the k*2^29+1, k*2^30+1, etc. results.

Now I ran the 2^63-2^64 search as well verifying and finding some new depth 28 primes.

Here is the complete list 2^62-2^64, my search criteria was root-of-two depth >=28. Root-of-two ">32" means the ones where you can keep taking the square root until you eventually reach 2.


EDIT: Ran k*2^28+1 from 2^28 to 2^62 and added 6 new entries to the list.

p			root-of-two	root-of-unity

57bit:
134442072127045633	>32		28

59bit:
569381225639182337	28		30

62bit:
2685240594304860161	>32		28
3207733922356002817	>32		28
4295962111625920513	>32		28
4588742434128658433	29		31

63bit:
4998875103153356801	>32		28
6249316677578653697	>32		28
7049705358222163969	>32		29
7128230911014862849	>32		29
7296706819391488001	>32		28
7301107626761256961	>32		29
7612253723611365377	>32		29
7931690348588302337	>32		29
8262758439982202881	>32		28
8585113616301162497	>32		29
8732723475115933697	>32		28

64bit:
9657578508943097857	>32		29
9857256047176056833	>32		29
9983467653877989377	29		30
10358097392821075969	28		32
10807913186691383297	>32		28
11191256683596742657	>32		28
11975441345180336129	>32		28
12046712698210091009	>32		28
12539167860364148737	>32		28
12548406757803687937	28		30
12552799470770716673	>32		28
13933613281743732737	>32		28
14251526830518435841	28		32
14358224048763174913	28		32
14380913871361671169	30		32
15613990314799267841	29		32
15910894009111805953	29		30
17061046145649737729	30		31
18195224058105167873	28		29