x^2+x+1

[jnml]

From the department of useless things:

Let \(p\) be an odd prime and \(M_p\) a Mersenne number \(2^p-1\).
Let's have a generating function \(x^2+x+1\) that for some \(x\) is equal to \(M_p\).
Example \(\{x, p\}\) are \(\{2, 3\}\), \(\{5, 5\}\) and \(\{90, 13\}\).
I don't know if the list is complete and haven't tried to search above \(p=23\).
Let's say we want to perform the Euler probable prime test, \(a^{n-1 \over 2} \equiv \pm 1 \pmod n\). Multiplying by \(a\) we get \(a^{{n-1 \over 2}+1} \equiv \pm a \pmod n\).
Plug in \(M_p\) for \(n\), we have \(a^{2^{p-1}} \equiv \pm a \pmod {M_p}\).
The exponent of \(a\) is a power of two, so the modular exponentiation reduces to repeated squaring of \(a \pmod {M_p}\).
Lets chose \(x\) for \(a\), provided it's coprime to \(M_p\). For the pairs listed above the modulus \(M_p\) is \(x^2+x+1\). So \(x\) must be coprime to \(x^2+x+1\), which it is, IINM. (\(x \equiv 1 \pmod {x^2+x+1}\))
First step: Squaring of \(x\) is \(x^2\) and that is \(-x-1 \pmod {x^2+x+1}\).
Second step: Squariing of \(-x-1\) is \(x^2+2x+1\) and that is \(x \pmod {x^2+x+1}\).
So we have a cycle of length two, an even number. \(p-1\) is also even so we conclude the final result of the repeated squaring is \(x\) without the need to actually compute it.
In other words we know analytically that any \(M_p\) generated by this particular function will pass the Euler probable prime test.
My intuition is that the pair list above is probably complete, sadly.
The reason for posting here is that maybe someone can come with a different generating function/functions that demonstrate similar cycles that can be verified analytically.
Maybe this generating function is the only one existing, I don't know and I don't have time today to think about it more.
But in the case there exist more/many/infinitely many generators with [reasonably short] cycles, the above described approach might become useful.

If I made a terrible mistake in this post, I'm very sorry. I didn't use paper and pencil so it would be 100% my fault.

----

The background of this is partially rooted in this previous thread.

[jnml]

I should probably add that any cycle of length \(n\), such that \(n|p-1\) is what will work.

[jnml]

A little exhaustive search for cycles of maximum length 2, limited to \(p \le 31\) (line 11), outputs:

$ time cycles 
M3: x=2, cycle=[4 2], len=2, x0=[{2 1}]
M3: x=3, cycle=[2 4], len=2, x0=[{3 1}]
M5: x=5, cycle=[25 5], len=2, x0=[{5 1}]
M5: x=6, cycle=[5 25], len=2, x0=[{2 1} {3 1}]
M7: x=19, cycle=[107 19], len=2, x0=[{19 1}]
M7: x=20, cycle=[19 107], len=2, x0=[{2 2} {5 1}]
M13: x=90, cycle=[8100 90], len=2, x0=[{2 1} {3 2} {5 1}]
M13: x=91, cycle=[90 8100], len=2, x0=[{7 1} {13 1}]
M17: x=21905, cycle=[109165 21905], len=2, x0=[{5 1} {13 1} {337 1}]
M17: x=21906, cycle=[21905 109165], len=2, x0=[{2 1} {3 2} {1217 1}]
M19: x=69492, cycle=[454794 69492], len=2, x0=[{2 2} {3 1} {5791 1}]
M19: x=69493, cycle=[69492 454794], len=2, x0=[{69493 1}]
M29: x=23884363, cycle=[408885410 489326097], len=2, x0=[{2803 1} {8521 1}]
M29: x=47544814, cycle=[408885410 489326097], len=2, x0=[{2 1} {23772407 1}]
M29: x=56556324, cycle=[489326097 408885410], len=2, x0=[{2 2} {3 2} {11 1} {251 1} {569 1}]
M29: x=109481575, cycle=[408885410 489326097], len=2, x0=[{5 2} {7 1} {625609 1}]
M29: x=127985501, cycle=[489326097 408885410], len=2, x0=[{7 2} {19 1} {23 1} {43 1} {139 1}]
M29: x=180910752, cycle=[408885410 489326097], len=2, x0=[{2 5} {3 1} {11 1} {171317 1}]
M29: x=189922262, cycle=[489326097 408885410], len=2, x0=[{2 1} {89 1} {1066979 1}]
M29: x=261351439, cycle=[489326097 408885410], len=2, x0=[{277 1} {433 1} {2179 1}]
M31: x=634005911, cycle=[1513477735 634005911], len=2, x0=[{7 1} {11 1} {8233843 1}]
M31: x=634005912, cycle=[634005911 1513477735], len=2, x0=[{2 3} {3 1} {109 1} {242357 1}]

real	8m2,825s
user	8m14,783s
sys	0m7,165s
$

In this limited set of data I'm observing, with no intent to generalize:

• There's at most one cycle of length 2 for any \(M_p\) with odd prime \(p\). Note that, for example, cycles [25 5] and [5 25] are the same cycle, just "entered" at different indices.
• If \(M_p\) is prime then

  - such 2-cycle exists
  - there are exactly two entry points
  - the entry points are next to each other
  

in #1 it was mistakenly stated that the very short list of \(\{x, p\}\) pairs generated by \(x^2+x+1\) is probably complete.
The mistake is in the fact, that the value of the generating function has to be considered \(\pmod {M_p}\).

Let's try another exhaustive (I hope) search:

? forprime(p=3,31,m=2^p-1;for(x=sqrtint(m)-1,m,y=Mod(x,m);if(y^2+y+1==m,print("{",x,",",p,"}");break)))
{2,3}
{5,5}
{19,7}
{90,13}
{21905,17}
{69492,19}
{634005911,31}
? ##
  ***   last result: cpu time 6min, 41,454 ms, real time 6min, 41,457 ms.
?

(The inner loop upper bound, \(m\), is probably unnecessarily high, but I haven't yet tried to really derive it, so this is a guessed "safe" value that should not miss any \(x\), I hope)

As you can see, the produced x values are those that can be shown to pass the Euler PRP test in the arrangement as in #1 and the same values that enter the 2-cycle above.

And as we're here in misc. math, let me ask a misc. math question:

Does \(x^2+x+1 \pmod {M_p}\) generate more (all?) Mersenne primes and possibly only Mersenne primes?

[LaurV]

First, inside of that "if", y is a Mod data type, while m is an int, so the y^2+y+1 is never m, but it is zero. You can compare with 0, to save a Mod() operation (the m on the right side of the equal sign is implicitly Mod()-ed to get zero). It is not significant, but let me be syntax nazi for a moment (well, your code will become few nanoseconds faster, hihi).

Then, your "conjecture" can be true even if you don't require p to be prime. It has to do with triangular numbers and multiples of mersenne numbers (i.e. every perfect number is also triangular). When you move from some x to its successor, you compute (x+1)^2+(x+1)+1, so in fact, you compute x^2+2x+1+x+1+1, i.e. (x^2+x+1)+(2x+2), i.e. in fact, you just add 2(x+1) to the former residue. As x grows, you just build a very HUGE triangular number (well, its double plus one). Your m will divide this, if it is prime. Try re-writing your one-liner like:

forprime(p=3,31,m=1<<p-1;z=m<<1;t=Mod(3,m);n=1;s=2;while(s<z,if(t==0,print("{"n","p"}");break);n++;t+=s+=2))
{2,3}
{5,5}
{19,7}
{90,13}
{21905,17}
{69492,19}
{634005911,31}
time = 12min, 59,642 ms.

(this has no multiplication, and no square root, and starting from 1, not from the square root, n is only to show, not used in calculation; no optimization was intended, as it is, it is even slower than yours, but it is functional equivalent - note that time comparison means nothing, my computer may have different speed, but this version is indeed slower than yours - it can however be optimized to be much faster, but that would be futile anyhow - of course, this "result" is not of any practical importance, we can't find primes like this - it would be much MUCH too slow).

[jnml]

Thanks for all the insights, LaurV!

I'll possibly come back on this later, if I make some progress. At the moment I'll only comment on the "unusuable for tests".
I agree as this runs out of steam quickly for \(p > 31\). But that code was an exhaustive search, not an attempt of a test algorithm.
OTOH, I would not rule out future usability of any piece of knowledge.

I only find it "interesting" if we can possibly claim, provided it's not an already known fact in disguise, that \(M_p\) is prime, or can be prime,
if/iff there exists \(x\) such that \(x^2+x+1 \equiv 0 \pmod{M_p}\). There's one more root to the quadratic but I messed up the formula.

It would be nice if it can be proven.

[Dr Sardonicus]

Originally Posted by jnml

<snip>
I only find it "interesting" if we can possibly claim, provided it's not an already known fact in disguise, that \(M_p\) is prime, or can be prime,
if/iff there exists \(x\) such that \(x^2+x+1 \equiv 0 \pmod{M_p}\).
<snip>

All that is required for the quadratic to be solvable (mod N) is that all prime factors of N be congruent to 1 (mod 3).

The smallest composite N = 2^p - 1 with this property is N = 2^37 - 1 = 137438953471. The solutions are:

Mod(7348910325, 137438953471)
Mod(74005089038, 137438953471)
Mod(63433864432, 137438953471)
Mod(130090043145, 137438953471)

[jnml]

Originally Posted by Dr Sardonicus

All that is required for the quadratic to be solvable (mod N) is that all prime factors of N be congruent to 1 (mod 3).

The smallest composite N = 2^p - 1 with this property is N = 2^37 - 1 = 137438953471. The solutions are:

Mod(7348910325, 137438953471)
Mod(74005089038, 137438953471)
Mod(63433864432, 137438953471)
Mod(130090043145, 137438953471)

Thanks for the reply! Can you please let me know how did you get the solutions?

[Dr Sardonicus]

Originally Posted by jnml

Thanks for the reply! Can you please let me know how did you get the solutions?

I used Pari-GP's factormod() to solve x^2 + x + 1 = 0 modulo each of the prime factors, then used chinese() to solve the congruence modulo their product 2^37 - 1.

There have been similar threads to this part of the Forum, though they were years ago.

I used 2^37 - 1 similarly as an example with x^2 + 6 if memory serves.

The fact that -3 is a quadratic residue mod M_p if M_p is prime is actually of secondary importance. More significant is that there's actually a simple formula for a square root of a quadratic residue r modulo a prime P congruent to 3 (mod 4), r^((P+1)/4). This is especially nice if P is one less than a power of 2.

Raising both sides of that formula to the 4th power gets rid of any pesky minus signs. Applying it to the present case, we see than if M_p is prime, we have

(-3)^(M_p + 1) == (-3)^2 (mod M_p), which may be rewritten

3^(M_p + 1) == 9 (mod M_p).

This is, if my understanding is correct, the PRP test currently in use, and usually applied after "enough" trial factorization and other standard methods have failed to find any factors.

AFAIK there are no known examples of composite M_p which have "passed" the PRP test.

[jnml]

Thanks again for your insights!

Originally Posted by Dr Sardonicus

I used 2^37 - 1 similarly as an example with x^2 + 6 if memory serves.

FTR: I think your memory serves you well ;-)