A new (but not really powerful) Mersenne theorem

Tony Reix · 2002-10-04, 07:20

smile
Hello, some years ago I've found and proven the following theorem about Mersenne numbers

For each Mersenne number Mq = 2^q - 1 (q prime > 4),
and for each pair (a,b) that is a solution of Mq = ab,
there does exist a unic pair (x,y) such that Mq = (8x)^2 - (3qy)^2 (I).

It doesn't seem it may help factorizing Mersenne numbers or proving their non-primality. But who knows ?

My original proof was quite long and inelegant.
I've recently received a shorter and more elegant proof.
It is hereafter. (As a game, try to prove the theorem before looking at the proof !)

-

Code:

-----------------------------------------------------------------------
Theorem
--------

For each Mersenne number Mq = 2^q - 1 (q prime > 4),
and for each pair (a,b) that is a solution of Mq = ab,
there does exist a unic pair (x,y) such that

                        Mq = (8x)^2 - (3qy)^2  (I).

-----------------------------------------------------------------------

Proof
------

Lets have      Mq = 2^q - 1 = ab   (with q odd prime)
and            A  = (a+b)/2
                B  = (a-b)/2
So, for each pair (a,b) is associated a unic pair (A,B).

Moreover       Mq = A^2 - B^2    [ = ((a+b)^2 - (a-b)^2)/4 = 4ab/4 = ab . ]

So we must demonstrate 
        8 | A
        3 | B
        q | B

Since 2 = -1 (mod 3), then 2^(2p+1) = (-1)^(2p+1) = -1 (mod 3),
and thus (with q = 2p+1)  ab = Mq = 2^(2p+1) - 1 = -2 = 1 (mod 3).
Since 1*1 = 2*2 = 1 (mod 3), then we have a = b (mod 3),
and thus 3 | a-b and 3 | B.

Since every prime divisor of Mq is equal to 1 (mod q),
thus we have a = b = 1 (mod q) , and q | a-b and then q | B.

Since every prime divisor of Mq is equal to +/- 1 (mod 8),
thus we have b = +/- 1 (mod 8) , and b^2 = 1 (mod 16).

Since (with q > 4) Mq = -1 (mod 16), then  ab = -1 (mod 16),
and thus 2bA = ab + b^2 = -1 + 1 = b(a+b) = 0 (mod 16).

Finally, since b is odd, that entails a+b = 0 (mod 16),
and 16 | a+b  , and thus 8 | A .

CQFD

-----------------------------------------------------------------------

Examples
---------

2^11 - 1  = 2047 =  1 * 2047 = (8*128)^2 - (3*11*31)^2 = (8X)^2 - (3qY)^2
                 = 23 *   89 = (8*  7)^2 - (3*11* 1)^2 = (8x)^2 - (3qy)^2

(X,Y)=(128,31) is the trivial solution of (I) for q=11 .
(x,y)=(  7, 1) is the non-trivial solution of (I) for q=11 .

2^23 - 1  = 8388607 = 47 * 178481 = (8*11158)^2 - (3*23*1293)^2
(x,y)=(11158,1293)  is the non-trivial solution of (I) for q=23 .

(x,y)=(47626997813,1894638183) is the non-trivial solution of (I) for q=67 .

-----------------------------------------------------------------------

Properties
-----------

Let say (X,Y) = (2^(q-4),(2^(q-1)-1)/(3q)) , the trivial pair.
Hereafter, the pair (x,y) represents any solution of (I), except (X,Y).

x and y have special properties

                        x =     X =   1/8   (mod q)
                        x =     X           (mod q^2)
                        x =     X = +/- 1   (mod 3)
                        x = +/- X           (mod 9)
if q = -1 (mod 6) then x =         +/- 2   (mod 9)
if q = +1 (mod 6) then x =         +/- 1   (mod 9)

                        y =     Y =   1     (mod 2)
                                Y = -1/(3q) (mod 2^5)
                        y = +/- Y           (mod 2^5)

---------------------------------------------------------------------

Maybeso · 2002-11-06, 23:34

Hello Tony,

After reading your elegent theorem 8) , I wondered exactly how powerful (or not) it was. The phrase that caught my attention was "... there does exist a unique pair ...". Since the trivial pair (X,Y) are maximum values for both x and y, Mq could be proven prime by failing to find a valid pair (x<X, y<Y).

I started with a program that took x == +/- X (mod 9) from 1 to X and solved for y. It was immediately obvious that completely solving for y for 2X/9 = 2^(q-5)/9 values of x would take too long.
I switched to finding x for y == +/-Y (mod 64) from 1 to M(q-1)/3q, M(q-1)/96q values. But again, the program was too slow.

Then a breakthrough :idea: I wasn't looking for y an integer, but y^2 a square

After rethinking my approach, here is the strategy of the current program:

From the Reix Theorem ;) : Mq = (8x)^2 – (3qy)^2 (1).

64x^2 = 2^q – 1 + (3qy)^2
x^2 = (2^q + (3qy)^2 – 1)/64
x^2 = 2^(q–6) + (9q^2y^2 – 1)/64 (2).

[code:1]Let [P] equal the first 40 odd primes.
For each p in [P], precalculate the quadratic residues.
Then for each q,
For each p in [P],
precalculate 9q^2 (mod 64p),
precalculate 2^(q-6) (mod p) using power mod.
Find smallest m such that: 9q^2m^2 == 1 (mod 64). then m == Y (mod 64).
For each y = 32n +/- m, up to Y,
For each p in [P],
z == 2^(q-6) (mod p) + [9q^2 (mod 64p) * y^2 (mod 64p) – 1]/64 (mod p).
If z isn't a quadratic residue of p, skip to next y since z isn't a square.
If z is a square, then Mq is composite.
If no squares are found, then Mq is prime. [/edit][/code:1]

This technique has several advantages.[list]The inner loop exits quickly for most non-squares.
The inner loop consists of 2 multiplies, a divide, 3 mods, and 2 adds.The other operations are already done.
All of the values remain small except y, and y can be handled as n (mod 64p)*32 + m.
We don't mess with the two largest/hardest calculations – x^2, and sqrt(x^2).[/list:u]However, it still has problems.[list]The value of Y gets very big very quickly; this means that y, n do too.
It still has to check M(q-1)/96q values.
Composite Mq's are not actually factored. That requires finding x^2, x, and then (a, b).[/list:u]In order for this program to be viable, the search space of y must be reduced, perhaps by further manipuling the properties of (x, y). Have you proven any corollaries or lemmas that might help? ;) Is there a quicker way to test whether an integer is a square? Most that I can think of require the entire integer, rather than a modulo.

This will give you an idea of how fast y grows:[code:1]
M( 7) is prime. Y = 30
M( 13) is prime. Y = 105
M( 17) is prime. Y = 1285
M( 19) is prime. Y = 4599
M( 31) is prime. Y = 11545611
M( 61) is prime. Y = 6300117511512825 [/code:1]

Tony Reix · 2002-11-07, 10:54

Hello Maybeso,
I'm very happy someone is using my theorem for trying to find some new way of proving primality of Mq !
I think I've already tried to do that, but I stopped quickly. I think I didn't try your way, but I haven't still understood all your work (I'm not a mathematician, only an engineer in computer science).

I have some comments
- I didn't realize that (X,Y) were maximum values for (x,y). Maybe you should prove it.
- As you say, your test is useful only if the number of y,n to be tested don't increase too fast. I'm afraid it will, with big values for q. It may take years or more than the life of the universe. Nevertheless, you can try using properties of x,y
At the end of my note, I gave some properties I've found between (X,Y) and (x,y), but without any proof. I think your algorithm already uses some of these properties. Maybe not all of them. What about
x = 1/8 (mod q)
x = X (mod q^2)
x = +/- 1 (mod 3)
?
If you provide me with some email address, I will send you a .dvi file describing my (inelegent) proof plus the (maybe also inelegent, but correct) proof for the (x,y) properties.

About proving that z is a square, I can only think about testing the last (1 or more) digits (mod 10...)
1 0 1 4 5 6 9 (6/10)
2 00 01 04 09 16 21 24 25 29 36 41 44 49 56 61 64 69 76 81 84 89 96 (22/100)
3 000 ... ( 159/1000)
4 0000 ... (1044/10000)
eliminating nearly 9 numbers out of 10 with 4 digits.

Please let me know if you find faster algorithms. (and provide me your email. Mine is Tony.Reix@bull.net).

Regards,

Tony

Maybeso · 2002-11-08, 19:50

Hello Tony,

My email address is also in my profile. mailto:Bruce_Leenstra@adp.com

I haven't proven that (X, Y) are maximum values, but here is a first try.
8x = A = (a + b)/2, and a = Mq/b, so x = (Mq/b + b)/16

The maximum integer value for the right hand side is when b = 1 or b = Mq, so
x = (Mq + 1)/16 = 2^q / 2^4 = 2^(q-4) = X

Similarly,
3qy = B = (a - b)/2, y = (Mq/b - b)/6q

again, the maximum occurs when b = 1, so
y = (Mq - 1)/6q = 2*M(q - 1)/6q = (2^(q-1) - 1)/3q = Y

therefore the trivial pair (X, Y) is the maximum integer solution to
Mq = (8x)^2 - (3qy)^2

(X, Y) is trivial because it gives the solution (a, b) = (Mq, 1).

A curious observation: Normally I have my program stop after the first pair is found, but I noticed that on Mq with multiple factors, it was finding the largest one. So I had it continue to Y, and it listed the factors in reverse order. Then I realised that since y is proportional to the difference between the factors, it would find the smallest *difference* first. Which woud be p1p2p3 - p4, and then p1p2p4 - p3, and so on.

ie. for q = 29, I get (1103*233, 2089), (2089*233, 1103), (2089*1103, 233); in that order.

Thanks for the additional properties you listed, they should help reducing the number of y I need to test. Although they relate to x, I can use substitution and exclude a few y values. I'm trying to derive rules for y itself; the best would be of the form y == +/- c (mod q^d). That would filter more y as q got larger.

And uh, as for not understanding the math, um, what are the rules for x = 1/8 (mod q)? or more generally x = a^(-n) (mod q) ??
I know that x = a (mod q) ==> x^n = a^n (mod q), but does that extend to negative n?

My first method to prove z a square was the digit checking that you suggested, (which is the quadratic residues of 10^n). But it didn't screen out all non-squares, so I switched to a sequential list of primes. I thought that primes would be more efficient. Now I'll write a program to generate the residue set for every number up to 10000, then select those with the smallest sets. (edit) turns out 300 is overkill (/edit)

The good thing is, I can make my sequence as long as I want. The program only goes thru all of them when z is a square -- and that means Mq is composite, so its done. It doesn't matter if a valid value of y takes a while, as long as invalid ones exit quickly.
Using the prime list, I kept track of how far in z made it before failing:
While testing M(59), one value of z was a square modulo the first 29 primes but not the 30th! The average depth is very small. The deepest for a given q is around 15 - 22.

My programs are in pascal, so they're longer than they would be in c, but I find it easier to debug. If you want to see snippets of the interesting parts, let me know. I'm moving this conversation off of the forum, or do you think we should continue posting it there?

Regards,

Bruce

Maybeso · 2002-11-12, 20:50

Hi Tony,

I posted my previous email to the forum, (I confess, I made one or two changes).
From now on I will post them all.

Excellent paper!
I'm going to take time to study it - the formal terminology is a little confusing, but after a while it is ok. The effort is worth it, I find that the more rigorous logic helps me close the gaps in my understanding.

Quote:

Originally Posted by Tony

Bruce wrote:
{ But it didn't screen out all non-squares, so I switched to a sequential
{ list of primes.

This is the part I don't understand.
What is the aim of these primes ?

My apologies, I should have said a sequential list of primes and their quadratic residues. On the first pass the inner loop calculates the right-hand-side of my equation in small pieces (mod 3). If it is not a residue mod 3, it can't be a square and the program drops out and tries with another value of y. If it is a residue, the program moves on to 5, then 7, ... etc.

I had thought that primes were better, but it turns out that #residues(p1p2p3p4...) <= sum(#residues(pi)). So composites filter out non-squares more quickly, (fewer false hits).

Quote:

Originally Posted by Tony

Do you have figures about the time needed for proving that some prime q doesn't build a prime Mq ?
So you can build some picture of the time needed for bigger q. Did you try with big q ?

No, I don't have any actual numbers. Right now it takes > 1 hour to prove M(59) composite and overnight to prove M(61) prime. :( It would take just as long for a composite Mq = a*b, q big, if (a - b) > 2^61.

That is why I am returning to y^2 = f(x^2): I have to check 2 in 32 y values, but now using x = X (mod q^2) I think I can drop it down to only 2 in 9q^2 x values. 2X/9q^2 = 2^(q-3)/9q^2 < 2Y/32 ~= 2^(q-5)/3q
Obviously it would be great to drop below sqrt(Mq) ~ 2^(q/2).

Quote:

Originally Posted by Tony

Bruce wrote:
{ ... what are the rules for x = 1/8 (mod q)? ...
In my paper, I give a way for computing it quickly :
1/8 (mod q) = [(8 - (q mod 8 ))*q + 1] / 8

Ah, this is great! Now I can precalculate 1/9q^2 (mod r), r in [N & residues]
--> No more mixing (mod 9q^2*r) with (mod r) to handle divisions in the inner loop.

Shrink that inner loop!

Regards,

Bruce

Tony Reix · 2002-11-13, 17:19

Hi Bruce,

About my formula for 1/8
1/8 (mod q) = [(8 - (q mod 8 ))*q + 1] / 8
I'm afraid I can't remember how I found it and if I proved it ...
Just now, I've experimented with all 9999 first odd primes and that works well. But experimenting is not a proof ....

The formula
1/n (mod q) = [(n - (q mod n))*q +1] / n
works well also for n 2 3 4 6 8 12 24 , with q prime > 3 .
But it fails with all other n (5 7 9 10 11 13 ...) .
I don't understand why ...

Tony

ewmayer · 2002-11-13, 20:07

The short answer:

We desire 2^(-1) = x, defined by 2*x == 1 mod q.
For q odd, (q+1) is even, hence x = (q+1)/2.
Simply cube x (modulo q) to get 2^(-3) mod q.

The long answer:

Remember, a multiplicative inverse for all integers in [1,q-1]
exists only if arithmetic modulo q defines a field, i.e. if q is prime.
However, for composite q we still will have mult. inverses of some subset
of the integers < q, as above, with 2 - the formula still works if, e.g.
q = 9. If an inverse of some general number is needed (and no closed-form
expression for it can be found), we can use the extended GCD algorithm to
find it: x = n^(-1) mod q satisfies n*x == 1 mod q. Taking a GCD(n,q)
amounts to finding a pair of integers a and b such that a*n + b*q = 1.
Modulo q, the b*q term drops out, and thus a*n == 1 mod q, i.e. the
a-coefficient is our desired modular inverse x.
If all we are interested in is GCD(n,q) we don't bother to accumulate
the coefficients a and b, since these aren't needed at the end - that's
the regular GCD algorithm. If we want a mod inverse, we need to do some
extra work to get a and b at the end - that's the extended GCD algorithm.
Here's a small example: to find 2^(-1) modulo 7, we calculate as follows:

gcd(2,7):
7/2 = 3, 7-3*2 = 1, so now have gcd(2,1).
In matrix form this is
[code:1]
|2| = |+1 0|*|2|
|1| |-3 +1| |7|
[/code:1]
Next, we need gcd(2,1):
2/1 = 2, 2-2*1 = 0, and we're done, at least as far as the value of the
gcd is concerned - gcd(7,2) is the smaller of the 2 inputs, hence = 1.
In matrix form this is
[code:1]
|0| = |+1 -2|*|2|
|1| | 0 +1| |1|
[/code:1]
To get the final coeeficients matrix, simply multiply the 2 above matrices
together:
[code:1]
|0| = |+1 -2|*|2| = |+1 -2|*|+1 0|*|2| = |+7 -2|*|2|
|1| | 0 +1| |1| | 0 +1| |-3 +1| |7| |-3 +1| |7|
[/code:1]
It's the row of coefficients (a, b) = (-3, +1) that give a linear combination
summing to one that interests us here - if (a*2 + b*7) = 1, writing this mod 7
eliminates the b-term, yielding a*2 == 1 mod 7, hence a = 2^(-1) == -3 == +4 mod 7.
(And for q = 7, we verify that (q+1)/2 == 4.)

p.s: If q has some special form (besides merely being odd), there may
be further simplifications we can use. For example, if q is a Mersenne
number, q = 2^p - 1, there are simple closed-form expressions for the
inverse powers of 2 mod q, based on the fact that 2^p == 1 mod q.
For example, our formula for 2^(-1) mod q simply becomes

2^(-1) == (q+1)/2 = 2^(p-1). Any desired inverse power of 2 can be gotten
by powering this result, i.e. for any integer k,
2^(-k) == [2^(p-1)]^k = 2^[k*(p-1)] == 2^(-k mod p),
where the result of the mod-p is understood to be nonnegative.

Tony Reix · 2002-11-15, 14:52

Hi ewmayer,
Thanks for explaining how to compute n^(-p) (mod q).
I think I read that long time ago, but I forgot ...
Your solution for 2(^-3) (mod q) is simpler than my formula.
But, for n prime, one has to use the extended GCD algorithm,
as you said.
I've played a little more with 1/n and found the following:

q is prime.
1 < n < q
Let's call: n^(-1) = 1/n .
Let say :
n*1/n = 1 (mod q) = 1 + k*q .
Thus k < q .
n*1/n (mod n) = 0
(1 + k*q) (mod n) = 1 + k*(q mod n)
Thus: k*(q mod n) = -1 (mod n)
If n not prime, then some (q mod n) have an inverse (mod n).
If n prime, then all (q mod n) have an inverse (mod n).
And thus:
k = -1/(q mod n) (mod n)

Thus, finding k such that:
n*1/n = 1 + k*q
doesn't depend on q, but on: (q mod n).
That means that k can be computed for all b ( b = (q mod n) with 1<b<n )
independently to any q, put in some table, and then it can be used for any q.

As an example, for n = 29, with c = (b*k+1)/29, we have:
[code:1]
b: 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
k: 28 14 19 07 23 24 04 18 16 26 21 12 20 02 27 09 17 08 03 13 11 25 05 06 22 10 15 01
c: 01 01 02 01 04 05 01 05 05 09 08 05 09 01 14 05 10 05 02 09 08 19 04 05 19 09 14 01
[/code:1]

We can check that:
(q mod 29)*k = 29*K-1
Also, we have (1<b<n):
k[b] + k[29-b] = 29

Thus computing 1/29 (mod 877) is:
877 mod 29 = 7
By the table : k = 4
1/29 = (4*877 + 1)/29 = 121

For someone studying 1/n (mod q) with n known and for many q, that seems faster.
Tony

Tony Reix · 2002-11-15, 15:01

I forgot to say that, for n = 8, the table is:
[code:1]
1 2 3 4 5 6 7
7 0 5 0 3 0 1
[/code:1]

Since, with q prime, (q mod 8) is odd, then it's easy to say that the formula I gave is proven:
[code:1]
8-1 = 7
8-3 = 5
8-5 = 3
8-7 = 1

1/n = (q*(8 - (q mod8)) + 1 ) / 8
[/code:1]
Tony

Maybeso · 2002-11-15, 18:57

Hi Tony,

Yes, your formula: 1/8 (mod q) = [(8 - (q mod 8 ))*q + 1] / 8
works perfectly for q prime. However, when I tried to use it to find
1/9q^2 (mod n), q prime,
it failed on many values of q & n, as both of you would expect.

Fortunately, I only need to find 200+ workable values for n, so I made a table and picked n with inverses for all values of q.

Bruce

ps. I haven't brought the timings for my program down by much. :(

Tony Reix · 2002-11-19, 08:15

Hi,
After some more search about computing inverses, I've found:
[code:1]
We want to compute 1/p (mod q).
0 < p < q
If: p * 1/p ~= 1 (mod q) (where 1/p is the inverse of p)
with: 0 < 1/p < q (1),
then there is some k such that: p * 1/p = 1 + k*q .
We have: (1) ==> 0 < p * 1/p < pq ==> -1 < p * 1/p - 1 < pq - 1 < pq
and thus: 0 =< k=(p * 1/p - 1)/q < p .
Since: p * 1/p (mod p) ~= 0 (mod p)
and since: 1 + k*q (mod p) ~= 1 + k*(q mod p) (mod p)
then : k ~= -1/(q mod p) (mod p)
and, since: 0= < k < p, finaly: k = p - 1/(q mod p)
and 1/p = (1 + q*(p - 1/(q mod p)))/p .

If we repeat the same operation with:
(q',p') = (p,(q mod p))
till p' == 0 (or 1),
the values of q' and p' are the rest of the GCD algorithm:
a = q
b = p
c = a mod b
d = b mod c
e = c mod d
etc.

Then:
1/b mod a = (1 + a*(b - (1 + b*(c - (1 + c*(d - (1 + d*(e - 1))/e))/d))/c))/b
After simplifying, we get:
1/b mod a = a + (1 - a*(1 - b*(1 - c*( 1 - d)/e )/d )/c )/b (I)
or
1/b mod a = a + (cde + a(-de + be -bc +bcd))/(bcde)
= a + (cde +a(b(e - c(1-d)) - de)/(bcde)

Example:
a = q = 61
b = p = 13
c = 61 mod 13 = 9
d = 13 mod 9 = 4
e = 9 mod 4 = 1

1/13 = 61 + (1 - 61*(1 - 13*(1 - 9)/4)/9)/13 = 61 - 14 = 47 = 1 + 10*61

With the optimized GCD (|r| =< (b-1)/2 , where a = q*b + r) :
a = q = 61
b = p = 13
c = 61 mod 13 = -4
d = 13 mod 4 = 1

1/13 = 61 + (1 - 61*(1 - 13/(-4))/13 = 61 - 14 = 47
[/code:1]

I'm not sure that the formula (I) can be derived from the usual
computation of the (a,b) coefficients of the ax + by = 1 Bezoux
formula. And I'm not sure it helps. Who knows ?

[code:1]
1/b mod a = a + (1 - a*(1 - b*(1 - c*( 1 - d*(.....))/e )/d )/c )/b (I)
[/code:1]

2002-11-07, 10:54	#3
Tony Reix Oct 2002 33₁₀ Posts	Using my theorem as a proof ? why not ! Hello Maybeso, I'm very happy someone is using my theorem for trying to find some new way of proving primality of Mq ! I think I've already tried to do that, but I stopped quickly. I think I didn't try your way, but I haven't still understood all your work (I'm not a mathematician, only an engineer in computer science). I have some comments - I didn't realize that (X,Y) were maximum values for (x,y). Maybe you should prove it. - As you say, your test is useful only if the number of y,n to be tested don't increase too fast. I'm afraid it will, with big values for q. It may take years or more than the life of the universe. Nevertheless, you can try using properties of x,y At the end of my note, I gave some properties I've found between (X,Y) and (x,y), but without any proof. I think your algorithm already uses some of these properties. Maybe not all of them. What about x = 1/8 (mod q) x = X (mod q^2) x = +/- 1 (mod 3) ? If you provide me with some email address, I will send you a .dvi file describing my (inelegent) proof plus the (maybe also inelegent, but correct) proof for the (x,y) properties. About proving that z is a square, I can only think about testing the last (1 or more) digits (mod 10...) 1 0 1 4 5 6 9 (6/10) 2 00 01 04 09 16 21 24 25 29 36 41 44 49 56 61 64 69 76 81 84 89 96 (22/100) 3 000 ... ( 159/1000) 4 0000 ... (1044/10000) eliminating nearly 9 numbers out of 10 with 4 digits. Please let me know if you find faster algorithms. (and provide me your email. Mine is Tony.Reix@bull.net). Regards, Tony

2002-11-08, 19:50	#4
Maybeso Aug 2002 Portland, OR USA 274₁₀ Posts	The Reix Theory: Improving the algorithm. Hello Tony, My email address is also in my profile. mailto:Bruce_Leenstra@adp.com I haven't proven that (X, Y) are maximum values, but here is a first try. 8x = A = (a + b)/2, and a = Mq/b, so x = (Mq/b + b)/16 The maximum integer value for the right hand side is when b = 1 or b = Mq, so x = (Mq + 1)/16 = 2^q / 2^4 = 2^(q-4) = X Similarly, 3qy = B = (a - b)/2, y = (Mq/b - b)/6q again, the maximum occurs when b = 1, so y = (Mq - 1)/6q = 2M(q - 1)/6q = (2^(q-1) - 1)/3q = Y therefore the trivial pair (X, Y) is the maximum integer solution to Mq = (8x)^2 - (3qy)^2 (X, Y) is trivial because it gives the solution (a, b) = (Mq, 1). A curious observation: Normally I have my program stop after the first pair is found, but I noticed that on Mq with multiple factors, it was finding the largest one. So I had it continue to Y, and it listed the factors in reverse order. Then I realised that since y is proportional to the difference between the factors, it would find the smallest difference* first. Which woud be p1p2p3 - p4, and then p1p2p4 - p3, and so on. ie. for q = 29, I get (1103233, 2089), (2089233, 1103), (2089*1103, 233); in that order. Thanks for the additional properties you listed, they should help reducing the number of y I need to test. Although they relate to x, I can use substitution and exclude a few y values. I'm trying to derive rules for y itself; the best would be of the form y == +/- c (mod q^d). That would filter more y as q got larger. And uh, as for not understanding the math, um, what are the rules for x = 1/8 (mod q)? or more generally x = a^(-n) (mod q) ?? I know that x = a (mod q) ==> x^n = a^n (mod q), but does that extend to negative n? My first method to prove z a square was the digit checking that you suggested, (which is the quadratic residues of 10^n). But it didn't screen out all non-squares, so I switched to a sequential list of primes. I thought that primes would be more efficient. Now I'll write a program to generate the residue set for every number up to 10000, then select those with the smallest sets. (edit) turns out 300 is overkill (/edit) The good thing is, I can make my sequence as long as I want. The program only goes thru all of them when z is a square -- and that means Mq is composite, so its done. It doesn't matter if a valid value of y takes a while, as long as invalid ones exit quickly. Using the prime list, I kept track of how far in z made it before failing: While testing M(59), one value of z was a square modulo the first 29 primes but not the 30th! The average depth is very small. The deepest for a given q is around 15 - 22. My programs are in pascal, so they're longer than they would be in c, but I find it easier to debug. If you want to see snippets of the interesting parts, let me know. I'm moving this conversation off of the forum, or do you think we should continue posting it there? Regards, Bruce

2002-11-15, 14:52	#8
Tony Reix Oct 2002 3×11 Posts	Pre-compute k for finding 1/n (mod q) Hi ewmayer, Thanks for explaining how to compute n^(-p) (mod q). I think I read that long time ago, but I forgot ... Your solution for 2(^-3) (mod q) is simpler than my formula. But, for n prime, one has to use the extended GCD algorithm, as you said. I've played a little more with 1/n and found the following: q is prime. 1 < n < q Let's call: n^(-1) = 1/n . Let say : n1/n = 1 (mod q) = 1 + kq . Thus k < q . n1/n (mod n) = 0 (1 + kq) (mod n) = 1 + k(q mod n) Thus: k(q mod n) = -1 (mod n) If n not prime, then some (q mod n) have an inverse (mod n). If n prime, then all (q mod n) have an inverse (mod n). And thus: k = -1/(q mod n) (mod n) Thus, finding k such that: n1/n = 1 + kq doesn't depend on q, but on: (q mod n). That means that k can be computed for all b ( b = (q mod n) with 1<b<n ) independently to any q, put in some table, and then it can be used for any q. As an example, for n = 29, with c = (bk+1)/29, we have: [code:1] b: 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 k: 28 14 19 07 23 24 04 18 16 26 21 12 20 02 27 09 17 08 03 13 11 25 05 06 22 10 15 01 c: 01 01 02 01 04 05 01 05 05 09 08 05 09 01 14 05 10 05 02 09 08 19 04 05 19 09 14 01 [/code:1] We can check that: (q mod 29)k = 29K-1 Also, we have (1<b<n): k[b] + k[29-b] = 29 Thus computing 1/29 (mod 877) is: 877 mod 29 = 7 By the table : k = 4 1/29 = (4877 + 1)/29 = 121 For someone studying 1/n (mod q) with n known and for many q, that seems faster. Tony

2002-11-15, 15:01	#9
Tony Reix Oct 2002 3·11 Posts	Proof for 1/8 formula. I forgot to say that, for n = 8, the table is: [code:1] 1 2 3 4 5 6 7 7 0 5 0 3 0 1 [/code:1] Since, with q prime, (q mod 8) is odd, then it's easy to say that the formula I gave is proven: [code:1] 8-1 = 7 8-3 = 5 8-5 = 3 8-7 = 1 1/n = (q*(8 - (q mod8)) + 1 ) / 8 [/code:1] Tony

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
2 holes in bizarre theorem about composite Mersenne Numbers	wildrabbitt	Math	120	2016-09-29 21:52
Størmer's theorem	grandpascorpion	Factoring	0	2006-09-10 04:59
Proof of LL theorem	Kevin	Math	53	2006-02-21 14:54
Powerful Numbers programming challenge	geoff	Programming	31	2005-02-20 18:25
Fermat,s Theorem	devarajkandadai	Miscellaneous Math	3	2004-06-05 10:15

2002-11-06, 23:34	#2
Maybeso Aug 2002 Portland, OR USA 2×137 Posts	Hello Tony, After reading your elegent theorem 8) , I wondered exactly how powerful (or not) it was. The phrase that caught my attention was "... there does exist a unique pair ...". Since the trivial pair (X,Y) are maximum values for both x and y, Mq could be proven prime by failing to find a valid pair (x<X, y<Y). I started with a program that took x == +/- X (mod 9) from 1 to X and solved for y. It was immediately obvious that completely solving for y for 2X/9 = 2^(q-5)/9 values of x would take too long. I switched to finding x for y == +/-Y (mod 64) from 1 to M(q-1)/3q, M(q-1)/96q values. But again, the program was too slow. Then a breakthrough :idea: I wasn't looking for y an integer, but y^2 a square After rethinking my approach, here is the strategy of the current program: From the Reix Theorem ;) : Mq = (8x)^2 – (3qy)^2 (1). 64x^2 = 2^q – 1 + (3qy)^2 x^2 = (2^q + (3qy)^2 – 1)/64 x^2 = 2^(q–6) + (9q^2y^2 – 1)/64 (2). [code:1]Let [P] equal the first 40 odd primes. For each p in [P], precalculate the quadratic residues. Then for each q, For each p in [P], precalculate 9q^2 (mod 64p), precalculate 2^(q-6) (mod p) using power mod. Find smallest m such that: 9q^2m^2 == 1 (mod 64). then m == Y (mod 64). For each y = 32n +/- m, up to Y, For each p in [P], z == 2^(q-6) (mod p) + [9q^2 (mod 64p) * y^2 (mod 64p) – 1]/64 (mod p). If z isn't a quadratic residue of p, skip to next y since z isn't a square. If z is a square, then Mq is composite. If no squares are found, then Mq is prime. [/edit][/code:1] This technique has several advantages.[list]The inner loop exits quickly for most non-squares. The inner loop consists of 2 multiplies, a divide, 3 mods, and 2 adds.The other operations are already done. All of the values remain small except y, and y can be handled as n (mod 64p)*32 + m. We don't mess with the two largest/hardest calculations – x^2, and sqrt(x^2).[/list:u]However, it still has problems.[list]The value of Y gets very big very quickly; this means that y, n do too. It still has to check M(q-1)/96q values. Composite Mq's are not actually factored. That requires finding x^2, x, and then (a, b).[/list:u]In order for this program to be viable, the search space of y must be reduced, perhaps by further manipuling the properties of (x, y). Have you proven any corollaries or lemmas that might help? ;) Is there a quicker way to test whether an integer is a square? Most that I can think of require the entire integer, rather than a modulo. This will give you an idea of how fast y grows:[code:1] M( 7) is prime. Y = 30 M( 13) is prime. Y = 105 M( 17) is prime. Y = 1285 M( 19) is prime. Y = 4599 M( 31) is prime. Y = 11545611 M( 61) is prime. Y = 6300117511512825 [/code:1]

2002-11-13, 17:19	#6
Tony Reix Oct 2002 21₁₆ Posts	Hi Bruce, About my formula for 1/8 1/8 (mod q) = [(8 - (q mod 8 ))q + 1] / 8 I'm afraid I can't remember how I found it and if I proved it ... Just now, I've experimented with all 9999 first odd primes and that works well. But experimenting is not a proof .... The formula 1/n (mod q) = [(n - (q mod n))q +1] / n works well also for n 2 3 4 6 8 12 24 , with q prime > 3 . But it fails with all other n (5 7 9 10 11 13 ...) . I don't understand why ... Tony

2002-11-13, 20:07	#7
ewmayer ∂²ω=0 Sep 2002 República de California 103×113 Posts	The short answer: We desire 2^(-1) = x, defined by 2x == 1 mod q. For q odd, (q+1) is even, hence x = (q+1)/2. Simply cube x (modulo q) to get 2^(-3) mod q. The long answer: Remember, a multiplicative inverse for all integers in [1,q-1] exists only if arithmetic modulo q defines a field, i.e. if q is prime. However, for composite q we still will have mult. inverses of some subset of the integers < q, as above, with 2 - the formula still works if, e.g. q = 9. If an inverse of some general number is needed (and no closed-form expression for it can be found), we can use the extended GCD algorithm to find it: x = n^(-1) mod q satisfies nx == 1 mod q. Taking a GCD(n,q) amounts to finding a pair of integers a and b such that an + bq = 1. Modulo q, the bq term drops out, and thus an == 1 mod q, i.e. the a-coefficient is our desired modular inverse x. If all we are interested in is GCD(n,q) we don't bother to accumulate the coefficients a and b, since these aren't needed at the end - that's the regular GCD algorithm. If we want a mod inverse, we need to do some extra work to get a and b at the end - that's the extended GCD algorithm. Here's a small example: to find 2^(-1) modulo 7, we calculate as follows: gcd(2,7): 7/2 = 3, 7-32 = 1, so now have gcd(2,1). In matrix form this is [code:1] \|2\| = \|+1 0\|\|2\| \|1\| \|-3 +1\| \|7\| [/code:1] Next, we need gcd(2,1): 2/1 = 2, 2-21 = 0, and we're done, at least as far as the value of the gcd is concerned - gcd(7,2) is the smaller of the 2 inputs, hence = 1. In matrix form this is [code:1] \|0\| = \|+1 -2\|\|2\| \|1\| \| 0 +1\| \|1\| [/code:1] To get the final coeeficients matrix, simply multiply the 2 above matrices together: [code:1] \|0\| = \|+1 -2\|\|2\| = \|+1 -2\|\|+1 0\|\|2\| = \|+7 -2\|\|2\| \|1\| \| 0 +1\| \|1\| \| 0 +1\| \|-3 +1\| \|7\| \|-3 +1\| \|7\| [/code:1] It's the row of coefficients (a, b) = (-3, +1) that give a linear combination summing to one that interests us here - if (a2 + b7) = 1, writing this mod 7 eliminates the b-term, yielding a2 == 1 mod 7, hence a = 2^(-1) == -3 == +4 mod 7. (And for q = 7, we verify that (q+1)/2 == 4.) p.s: If q has some special form (besides merely being odd), there may be further simplifications we can use. For example, if q is a Mersenne number, q = 2^p - 1, there are simple closed-form expressions for the inverse powers of 2 mod q, based on the fact that 2^p == 1 mod q. For example, our formula for 2^(-1) mod q simply becomes 2^(-1) == (q+1)/2 = 2^(p-1). Any desired inverse power of 2 can be gotten by powering this result, i.e. for any integer k, 2^(-k) == [2^(p-1)]^k = 2^[k(p-1)] == 2^(-k mod p), where the result of the mod-p is understood to be nonnegative.

2002-11-15, 18:57	#10
Maybeso Aug 2002 Portland, OR USA 2×137 Posts	Hi Tony, Yes, your formula: 1/8 (mod q) = [(8 - (q mod 8 ))*q + 1] / 8 works perfectly for q prime. However, when I tried to use it to find 1/9q^2 (mod n), q prime, it failed on many values of q & n, as both of you would expect. Fortunately, I only need to find 200+ workable values for n, so I made a table and picked n with inverses for all values of q. Bruce ps. I haven't brought the timings for my program down by much. :(

2002-11-19, 08:15	#11
Tony Reix Oct 2002 3×11 Posts	Hi, After some more search about computing inverses, I've found: [code:1] We want to compute 1/p (mod q). 0 < p < q If: p * 1/p ~= 1 (mod q) (where 1/p is the inverse of p) with: 0 < 1/p < q (1), then there is some k such that: p * 1/p = 1 + kq . We have: (1) ==> 0 < p 1/p < pq ==> -1 < p * 1/p - 1 < pq - 1 < pq and thus: 0 =< k=(p * 1/p - 1)/q < p . Since: p * 1/p (mod p) ~= 0 (mod p) and since: 1 + kq (mod p) ~= 1 + k(q mod p) (mod p) then : k ~= -1/(q mod p) (mod p) and, since: 0= < k < p, finaly: k = p - 1/(q mod p) and 1/p = (1 + q(p - 1/(q mod p)))/p . If we repeat the same operation with: (q',p') = (p,(q mod p)) till p' == 0 (or 1), the values of q' and p' are the rest of the GCD algorithm: a = q b = p c = a mod b d = b mod c e = c mod d etc. Then: 1/b mod a = (1 + a(b - (1 + b(c - (1 + c(d - (1 + d(e - 1))/e))/d))/c))/b After simplifying, we get: 1/b mod a = a + (1 - a(1 - b(1 - c( 1 - d)/e )/d )/c )/b (I) or 1/b mod a = a + (cde + a(-de + be -bc +bcd))/(bcde) = a + (cde +a(b(e - c(1-d)) - de)/(bcde) Example: a = q = 61 b = p = 13 c = 61 mod 13 = 9 d = 13 mod 9 = 4 e = 9 mod 4 = 1 1/13 = 61 + (1 - 61(1 - 13(1 - 9)/4)/9)/13 = 61 - 14 = 47 = 1 + 1061 With the optimized GCD (\|r\| =< (b-1)/2 , where a = qb + r) : a = q = 61 b = p = 13 c = 61 mod 13 = -4 d = 13 mod 4 = 1 1/13 = 61 + (1 - 61(1 - 13/(-4))/13 = 61 - 14 = 47 [/code:1] I'm not sure that the formula (I) can be derived from the usual computation of the (a,b) coefficients of the ax + by = 1 Bezoux formula. And I'm not sure it helps. Who knows ? [code:1] 1/b mod a = a + (1 - a(1 - b(1 - c( 1 - d*(.....))/e )/d )/c )/b (I) [/code:1]