# Second part of the factorial sum divisibility question

Which primes $p$ divide the sum of factorials $1! + 2! + 3! + 4! + 5! + \cdots + (p-1)!$? This is related to my previous question.

-
Have you ever considered beginning one of your contributions "I'm just a simple caveman. Your prime numbers frighten and confuse me. But one thing I DO know is,..." etc... ? –  Stephen Mar 25 '13 at 1:38
Up to 100,000, the only values of $p$ (not necessarily prime) that satisfy this are $3, 9, 11, 33, 99$. –  Yoni Rozenshein Mar 25 '13 at 15:09
What's motivating the question? Fun or work (or both)? Surely this isn't some demented homework problem! –  Douglas B. Staple Mar 26 '13 at 12:12
@Thus, please don't tell that the margin of math.stackexchange.com is too small for your beautiful proof... –  vonbrand Mar 26 '13 at 16:23
I must somehow correct it! –  user67878 Mar 26 '13 at 17:10

Note that heuristically, if these numbers were equidistributed mod $p$ we would expect each one to be zero (i.e., expect $p$ to divide the result) with probability $\approx 1/p$; presuming that all the values are independent (which seems a reasonable assumption), we shoud expect to have $\sum_{p\lt n}\frac{1}{p}\approx M+\ln \ln n$ of them less than $n$, where $M$ is the Meissel-Mertens Constant — this sum is approximately 2.887 for $n=10^6$, so while we would heuristically 'expect' another one in that range it's not surprising that there aren't any, and despite the apparent lack of any more solutions, because of the divergence of the series the 'expected' number of solutions is still infinite!

Added: I put together my own software implementation in C++, confirming the results (that only $p=3$ and $p=11$ work) for $p\lt 5\times 10^5$ but also looking at the distribution of results $\mod p$ to see if any patterns arose. I computed the value of $\dfrac{\sum_{1}^{p-1}i!\bmod p}{p}$ (the division by $p$ serving to 'normalize' values into the range $[0, 1)$ ) for all $p$ in that range and then grouped them in bins; these are the plots with 100 bins and 256 bins. (These bin counts should be small enough relative to the prime values in consideration that there shouldn't be any substantial aliasing effects in binning the data.) The results are a bit scattershot but there doesn't appear to be any skew towards particular values (e.g. $\frac{p-1}2$) which would show up as a spike in the bins. If I have more time I may take a look at the square of the sum and particularly even the inverse of the sum $\pmod p$, but since the latter would actually require me writing a quick GCD algorithm to find the inverse it'll have to wait for later.

(Note: the sharp drop at the last element is an artifact of my adding a 0 value to the end of the bin data so my spreadsheet didn't automatically normalize the range; it's not an actual value.)

-
Interesting. Reminds me of an estimation of the number of Wieferich primes... –  Gottfried Helms Mar 26 '13 at 6:27
+1 for your heuristic. –  Douglas B. Staple Mar 26 '13 at 12:59
+1 But "we would expect each one to be prime with prob..." Shouldn't that be "to be a multiple of $p$"? –  leonbloy Mar 26 '13 at 18:27
@leonbloy Whoops, good catch - fixed now. –  Steven Stadnicki Mar 26 '13 at 18:40
Also, $M + \ln \ln n \ge 4$ for $n$ around $10^{18}$, so perhaps finding the next prime satisfying the property will require us to look at all primes up to that much. –  ShreevatsaR Mar 27 '13 at 10:42

I have a few small analytic comments and a relatively aggressive numerical search.

Firstly, I checked all odd $p$ (not necessarily prime) up to $10^6$ and found only the same five as Yoni above, namely 3, 9, 11, 33, and 99. This took ~1 hour using 1 core of my MacBook Pro, with an early unoptimized version of my code.

Secondly, of course we have: $$p\nmid\sum_{n=1}^{p-1}n!\quad \forall p\in\mathbb{N}, p \textrm{ even}$$ because all of the terms in the sum are even except the first one, so the sum is always odd.

Thirdly, a numerical search is sped-up slightly using the identity: $$\sum_{n=1}^{p-1} n! = 1+2[1+3[1+4[...[1+(p-1)]]]]$$ I've written software to check primes $p$ for $p|\sum_{n=1}^{p-1} n!$ using the above identity, and set it up to use multiple threads in an embarrassingly-parallel fashion. Using this software I've shown: $$\textrm{3 and 11 are the only primes }p<10^7\textrm{ that divide }\sum_{n=1}^{p-1} n!$$ If you want to compile it yourself then you'll have to link with the primesieve library.

#include <stdio.h>
#include <stdlib.h>
#include <primesieve/soe/PrimeSieve.h>

// To compile: g++ -O3 factorialSum.cpp -lprimesieve -o factorialSum
//         or: g++ -O3 -static factorialSum.cpp -lprimesieve -o factorialSum

// We perform arithmetic at u1 precision; p is restricted to u2 precision
typedef uint_fast64_t u1;
typedef uint_fast32_t u2;

// This evaluates $\sum_{n=1}^{p-1} n!$ based on the definition.
// We don't actually use this in our program below.
u2 factorialSum_modp_method1(u2 p){
u1 factorial_modp=1;
u1 factorialSum_modp=0;
for(u1 n=1; n<=p-1; n++){
factorial_modp*=n;
factorial_modp%=p;
factorialSum_modp+=factorial_modp;
factorialSum_modp%=p;
}
return factorialSum_modp;
}

// This evaluates $\sum_{n=1}^{p-1} n! = 1+2[1+3[1+4[...[1+(p-1)]]]]$ via the RHS of the expression.
// This is somewhat faster than 'method1' above.
u2 factorialSum_modp_method2(u2 p){
u1 factorialSum_modp=1;
for(u1 n=p-1; n>1; n--)
factorialSum_modp=(factorialSum_modp*n+1)%p;
return factorialSum_modp;
}

// This prints out p iff p divides $\sum_{n=1}^{p-1} n!$.
// We pass this as an arument to PrimeSieve.generatePrimes().
void check_factorialSum_modp(u2 p){
if(factorialSum_modp_method2(p)==0)
printf("%u\n", (unsigned)p);
}

void error_incorrectUsage(char *args[]){
printf("       numThreads is a natural number, 1, 2, ...\n");
printf("       blockSize is the interleave size between threads\n");
printf("       pMin and pMax are integers, the upper and lower limits on candidate primes p\n\n");
printf("e.g.:  %s 0 1 10000 0 100000\n", args[0]);
printf("       will check [0, 100000] for primes p that divide $\\sum_{n=1}^{p-1} n!$ \n");
printf("       (expression written LaTeX).\n\n");
printf("In order to use multiple threads, run multiple instances with different values of 'thread'.\n");
exit(1);
}

int main(int nargs, char *args[]){

// Parse command line arguments and check for proper usage
if( !(nargs==6) ) error_incorrectUsage(args);
char *error1, *error2, *error3, *error4, *error5;
u2 thread     = strtoul(args[1], &error1, 10);
u2 numThreads = strtoul(args[2], &error2, 10);
u2 blockSize  = strtoul(args[3], &error3, 10);
u2 pMin       = strtoul(args[4], &error4, 10);
u2 pMax       = strtoul(args[5], &error5, 10);
if(*error1 || *error2 || *error3 || *error4 || *error5) error_incorrectUsage(args);

// Sieve intervals [pLower, pUpper] that belong to this thread
PrimeSieve ps;
u2 pUpper=pLower+blockSize;
ps.generatePrimes(pLower, pUpper, check_factorialSum_modp);
}

}

-
Thank you kindly it will be very interesting to see whether there are more or not. –  user58512 Mar 26 '13 at 3:51
You're welcome. I had fun looking at this problem. Computing up to $10^7$ took 45 minutes using 16 threads on an 8-core Opteron. Pushing the limit up to $10^8$ would take about three days on that machine, which we can do if you're using this for research. Otherwise I think we've learned all we're going to learn from computation. –  Douglas B. Staple Mar 26 '13 at 17:44

The following is more a comment than a usable answer, but too long for the comment-box

Another oberservation which has a surprising symmetry, which seems to occur only if the modulus is a prime (and this is thus an observation in its own right).

We define the lower triangular matrix of the Eulerian numbers:

$\qquad \qquad \displaystyle E= \Tiny \left[ \begin{array} {rrrrr} 1 & . & . & . & . & . \\ 1 & 0 & . & . & . & . \\ 1 & 1 & 0 & . & . & . \\ 1 & 4 & 1 & 0 & . & . \\ 1 & 11 & 11 & 1 & 0 & . \\ 1 & 26 & 66 & 26 & 1 & 0 \end{array} \right]$

We observe, that the row sums accumulate to the factorials so $E \cdot V(1) = F$ where $V(1)$ is the column-vector of ones, gives the column-vector of factorials $F$. If we now add the rows up to the index of the prime $p$ in question and take this modulo $p$ then we have another expression of the sum. If we now change the order of computation: we compute first the column-sums (modulo p) then we get a rowvector with symmetric entries, seemingly iff p is prime. For instance, for $p=13$ (and note, that in Pari/GP the index begins at 1 and not at 0)

$\qquad \qquad V(1)^T_{p+1} ~ \cdot E_{p+1,p+1} \pmod p = \\\qquad \qquad \left[ \begin{array} {rrrrrrrrrrrrr} 1 & 2 & 3 & 11 & 6 & 12 & 5 & 12 & 6 & 11 & 3 & 2 & 1 & 0 \end{array} \right]$

and then we sum the resulting vector modulo(p) and subtract 1 because we have also the 0! in the sum.
I don't see actually what this might give us - and I can't go deeper in this - but it may also be significant (and possibly on the way to a -more- closed form) that the columns of the Eulerian triangle are expressible as sequences of leading terms of geometric series and their derivatives and so might have a closed form for their leading segments.
Anyway -perhaps we can make the resulting row and its modular sum simpler by some pattern.

Here are the resulting vectors for the first few primes: $\Tiny \begin{matrix} 2: & 1 & 1 & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ 3: & 1 & 2 & 1 & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ 5: & 1 & 2 & 3 & 2 & 1 & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ 7: & 1 & 2 & 3 & 1 & 3 & 2 & 1 & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ 11: & 1 & 2 & 3 & 10 & 8 & 8 & 8 & 10 & 3 & 2 & 1 & . & . & . & . & . & . & . & . & . & . & . & . & . \\ 13: & 1 & 2 & 3 & 11 & 6 & 12 & 5 & 12 & 6 & 11 & 3 & 2 & 1 & . & . & . & . & . & . & . & . & . & . & . \\ 17: & 1 & 2 & 3 & 13 & . & 13 & 15 & 9 & 3 & 9 & 15 & 13 & . & 13 & 3 & 2 & 1 & . & . & . & . & . & . & . \\ 19: & 1 & 2 & 3 & 14 & 15 & 1 & 7 & 8 & 15 & 10 & 15 & 8 & 7 & 1 & 15 & 14 & 3 & 2 & 1 & . & . & . & . & . \\ 23: & 1 & 2 & 3 & 16 & 9 & . & 15 & 12 & 18 & 12 & 6 & 17 & 6 & 12 & 18 & 12 & 15 & . & 9 & 16 & 3 & 2 & 1 & . \\ 29: & 1 & 2 & 3 & 19 & 24 & 9 & 18 & . & 28 & 23 & 7 & 24 & 22 & 14 & 6 & 14 & 22 & 24 & 7 & 23 & 28 & . & 18 & 9 \end{matrix}$

[update] If we decompose the columns of the Eulerian matrix in their component-series we're getting involved with the leading segments of geometric series and their derivatives - and that all modulo p. Because in the end we need the sum of the column-sums I've expressed the final partial sum in terms of its (truncated) series-components, here for the example prime p=5 . We get as the final result the sum of all elements of the following matrix $$\small \begin{bmatrix} 1\cdot 5^0 & 0\cdot 4^0 & 0\cdot 3^0 & 0\cdot 2^0 & 0\cdot 1^0 \\ 1\cdot 5^1 & -1\cdot 4^1 & 0\cdot 3^1 & 0\cdot 2^1 & 0\cdot 1^1 \\ 1\cdot 5^2 & -2\cdot 4^2 & 1\cdot 3^2 & 0\cdot 2^2 & 0\cdot 1^2 \\ 1\cdot 5^3 & -3\cdot 4^3 & 3\cdot 3^3 & -1\cdot 2^3 & 0\cdot 1^3 \\ 1\cdot 5^4 & -4\cdot 4^4 & 6\cdot 3^4 & -4\cdot 2^4 & 1\cdot 1^4\\ 1\cdot 5^5 & -5\cdot 4^5 & 10\cdot 3^5 & -10\cdot 2^5 & 5\cdot 1^5 \end{bmatrix} \pmod 5$$ where -because p is a prime- we have systematic zeros and congruences expressible in negative numbers $\pmod p$ and get this:

$$\small \left[ \begin{array} {rrrrrrrrrrr} 1 & . & .& . & . \\ . & 1\cdot 1^1 & . & . & . \\ .& -2\cdot 1^2 & 1\cdot 2^2 & . & . \\ . & 3\cdot 1^3 & -3\cdot 2^3 & 1\cdot 3^3 & . \\ . & -4\cdot 1^4 & 6\cdot 2^4 & -4\cdot 3^4 & 1\cdot 4^4\\ . & . & . & . & . \end{array} \right] \pmod 5$$ which can be generalized to other prime moduli p in the obvious way.
The idea with all this is to be of some advantage for a better/more intuitive expression of the sum perhaps by closed forms of the truncated geometric series (modulo p). Still we can recognize the row-sums as factorials (while modulo p), but the change of order of computation and the hopefully advantageous closed forms of the column sums might be helpful.

We can even more make from the modularity. The following represntation can analoguously be made even for higher primes p, but I use just the example before. The above matrix can be rewritten in complete congruence as $$\small \left[ \begin{array} {rrrrrrrrrrr} 1 & . & .& . & . \\ . & 1\cdot 1^1 & . & . & . \\ .& 3\cdot 1^2 & 1\cdot 2^2 & . & . \\ . & 3\cdot 1^3 & 2\cdot 2^3 & 1\cdot 3^3 & . \\ . & 1\cdot 1^4 & 1\cdot 2^4 & 1\cdot 3^4 & 1\cdot 4^4\\ . & . & . & . & . \end{array} \right] \pmod 5$$ and those simple binomials give then as columnsums $$\small \left[ \begin{array} {rrrrrrrrrrr} 1 & 1^1\cdot 2^3 & 2^2\cdot 3^2 & 3^3\cdot 4^1 & 1\\ \end{array} \right] \pmod 5$$ and if I put together each two terms reading from left from right then this is
$$\small \begin{array} {rrrrrrrrrrr} 2 \cdot 1 & + 2 \cdot 2^3 & + 2^2\cdot 2^2 \\ \end{array} \equiv 4\pmod 5$$ where by the definition in the OP we had to decrement by 1 to delete the $0!$ term.

If we write the odd prime $p=2q+1$ then the general formula (including the $0!$) becomes now: $$S_p \equiv \sum_{k=0}^{p-1} k! \equiv 2 \sum _{k=0}^{q-1} k^k (k+1)^{p-1-k} + (q \cdot -q)^q \pmod p$$

which might be finally accessible by a bit better knowledge about primitive roots modulo p.

-

Gottfried's answer and Steven's replies bring up the point that there's a choice to be made in a direct numerical search: either you can compute $(p-1)!$ $(\textrm{mod }p)$ separately for each $p$, or you can store $(p-1)!$ and then compute $p! = p(p-1)!$ to test the next $p$. This same tradeoff comes up in other similar problems, such as the search for Brown Numbers.

In order to check every $p<N$, the first approach requires $O(N^2)$ multiplies of numbers of order $p$; the second approach requires $N$ multiplies, but each multiply is very expensive, being between a number of order $p$ and one of order $p!$. Another consideration is that the second solution, storing $(p-1)!$ un-reduced, requires memory linear in $p$. I don't really want to see what happens when $(p-1)!$ no longer fits in cache. Furthermore, it's not obvious how to parallelize the second algorithm.

My impression is that if you want to check all $p<N$ and not just the primes $p$ it's fastest to store $(p-1)!$ for small $N \lesssim 10^6$, but for larger $N$ it's faster to compute $(p-1)!$ $(\textrm{mod }p)$ separately for each $p$. If you want to only check primes $p$, then the "first approach" gets a $O(\log(n))$ speedup and the "second approach" hardly benefits at all, meaning that the first approach is faster: see the "Update:" section below the code for an explanation.

Here is a C code that uses the GMP to implement this "second approach"; I implemented the "first approach" in my first answer above.

#include <gmp.h>
#include <time.h>

#define NMAX 100000

int main(){

// Print begin time & date
time_t beginTime;
struct tm * timeinfo;
time ( &beginTime );
timeinfo = localtime ( &beginTime );
gmp_printf("# Begin time/date is: %s", asctime (timeinfo) );

// Stores n! and sum(n!)
mpz_t cur_factorial;
mpz_t sum_factorial;
mpz_t sum_factorial_modn;
mpz_init_set_ui(cur_factorial, 1);
mpz_init_set_ui(sum_factorial, 1);
mpz_init(sum_factorial_modn);

for(unsigned long int n=2; n<=NMAX; n++){

// Check to see if n divides sum(m=1..n-1) m!
if( mpz_divisible_ui_p(sum_factorial, n) )
gmp_printf("%u\n", n);

// Update n!
mpz_mul_ui(cur_factorial, cur_factorial, n);

// Update sum(n!)

}

// Print end time & date; total run time
time_t endTime;
time ( &endTime );
timeinfo = localtime ( &endTime );
gmp_printf("# End time/date is:   %s", asctime (timeinfo) );
gmp_printf("# Total runtime of %.0f seconds.\n", difftime(endTime, beginTime) );

// Release and return
mpz_clear(cur_factorial);
mpz_clear(sum_factorial);
}


Update:

This code implementing the "second approach" is strictly slower than the one in my previous answer. The above code takes ~9s to for an upper limit to $10^5$, which only takes ~5s using the code in my previous answer. A search up to $10^6$ takes ~19 min with this code, compared to ~8 min with the other one. Part of the reason for this is that, with the other code, you get a huge speedup if you only test primes $p$. This code is set up to check all $n$, even the pointless even ones. This might not seem like a fair comparison, but modifying the above code to only check odd $n$ doesn't give any speed improvement; the point is that you anyway have to calculate $(n-1)!$ for every $n$ with this method, so you've already paid the cost; the divisibility check doesn't take any time compared with the multiplies.

-

In addition to @Douglas post, here is a very short (and I think: efficient) implementation in Pari/GP - however, I didn't run it to high N myself: [Update] I added an improved computation method at the end[/update]

MaxN=101   \\ set some upper limit for the computation
S=[1,1]    \\ contains current factorial and current sum
for(k=3,MaxN, S *= [k-1,k-1;0,1]; if(isprime(k),print([k,S[1] % k,S[2]%k])))


Result :

p | (p-1)! %p | S %p
--+-----------+-----
[3, 2, 0]
[5, 4, 3]
[7, 6, 5]
[11, 10, 0]
[13, 12, 9]
[17, 16, 12]
[19, 18, 8]
[23, 22, 20]
[29, 28, 16]
[31, 30, 1]
[37, 36, 4]
[41, 40, 3]
...


Up to the 1000'th prime MaxN=7919

MaxN = prime(1000)
...
...
[7919, 7918, 2882]


this needed 281 msec.

A Version using the modulo-computation is the following

gettime()
{for(j=1,1000,p=prime(j);
S=[1,1];
for(k=3,p,
S*=[k-1,k-1;0,1];
S=S % p
);
print([p,S[1] % p,S[2]%p])
)}
gettime()

%1124 = 9219   \\ this is msecs


[update] The following is only relevant, if we want a routine, which produces the list over all primes up to some limit.
It has then some advantage to do a compromise between the first method where we need all factorials and sums only once but which needs many bits for the representation of the numbers and the complete modulo-computation, which operates on small numbers but must do the whole run of computing factorials (residuals) for each prime such that we have a double loop.
If we compute and store the required factorials and partial sums in chunks of, say, $t=10$ terms, such that we have a list of precomputed values $$\small \begin{array} {r|l|l|} & f & S \\ \hline \\ c_1=& 10! & \sum_{k=1}^{10} k! \\ c_2=& {20!\over 10!} & {\sum_{k=11}^{20} k! \over 10!} \\ c_3=& {30!\over 20!} & {\sum_{k=21}^{30} k! \over 20!} \\ \ldots \end{array}$$ then we can compose our results for each prime by modular multiplication of that constants and a short remainder. Finetuning the parameter $t$ for higher allows then to find a compromise between storage/size of the stored numbers and number of steps in the inner loop, which contains the modular computation per prime.
Using some empirical data I have done an estimation with a quadratic curve in Excel, which says, for the list of primes up to 1000000 I need 90 min with blocksize $t=20$ , 30 min with $t=100$ and 15 min with $t=500$ .
(If it is wished I can append the Pari/GP code here)

-
This isn't likely to be very efficient for long because (AFAICT) it computes the exact values of the current factorial and current sum, and these will become very large numbers very quickly. It should be much better to work mod $p$ for each prime $p$ individually. –  Steven Stadnicki Mar 26 '13 at 5:44
@Steven: yes, you may be correct. However, to compute the list up to some MaxPrime this needs then time quadratic in the list length. Up to the 400'th prime (MaxN=2741) that modulo-method needed 3478 msec. I didn't check, whether the ratio of the timings change with higher MaxN/MaxPrime –  Gottfried Helms Mar 26 '13 at 6:12
That's a very good point - there's a lot to be said for the linear time. If I'm figuring right, though, note that your multiplications and additions will start to take time proportional to the length of the values which are roughly proportional to $n\log n$, so your solution is also quadratic in the list length in terms of bit-operations. –  Steven Stadnicki Mar 26 '13 at 6:25
@Steven: agreed! –  Gottfried Helms Mar 26 '13 at 6:29
@Steven You know what... I've encountered exactly this tradeoff between $n^2$ multiplies and reductions mod $n$ versus $n$ multiplies without reduction in the factorial part of the search from Brown Numbers. It came up somewhere on mersenneforum.org. The answer isn't at all obvious; I'm going to post an answer addressing this. –  Douglas B. Staple Mar 26 '13 at 12:08

This is what I thought about, it should be some comment, but this exceeds the character limit for comment.

If you take Wilson Theorem that

$(p-1)!\equiv -1 \mod{p}$, then

$(p-2)!\equiv 1 \mod{p}$.

If we take inverse of $a$ in $F_p$ as $\dfrac{1}{a}$, we know that $\dfrac{1}{a} + \dfrac{1}{b} = \dfrac{a+b}{ab}$ in $F_p$, it can taken as usual addition.

Then

$(p-3)! \equiv -\dfrac{1}{2}\mod{p}$

$(p-4)! \equiv -\dfrac{1}{2}\cdot\left(-\dfrac{1}{3}\right)\mod{p}$

$\cdots$

$1! \equiv (-1)^{p-3}\dfrac{1}{(p-2)!} \mod{p}$

Thus

$1!+2!+\cdots(p-3)! \equiv -\sum_{k=0}^{p-2}(-1)^k\dfrac{1}{k!}\mod{p}$

And we state without proof.

Lemma:

Consider the fraction $\dfrac{p}{q}$ with $q = m!$, then the closest number to $\exp(-1)$ with form $\dfrac{p}{q}$ should satisfy

$$\dfrac{p}{q} = \sum_{j=0}^m (-1)^j\dfrac{1}{j!}$$

Then

$-\sum_{k=0}^{p-2}(-1)^k\dfrac{1}{k!} \equiv - \mathrm{floor}\left[\dfrac{(p-2)!}{\exp(1)}\right]\mod{p}$

However, this cannot be checked for large $p$, since the constant $\exp(1)$ is not accurate enough. I simply checked for smaller primes.

However, I felt that there should be other solutions than 3 and 11, otherwise, this can be a very interesting property for the constant.

-
Unfortunately, this doesn't really provide much help because it's very hard to translate back from the 'analytical' description to anything useful about the residues; even as simple a question as 'is $\frac{1}{n} \bmod p$ less than or greater than $\frac{p}{2}$?' doesn't really (AFAIK) have a clear large-scale 'analytical' structure to it. –  Steven Stadnicki Mar 26 '13 at 21:07
Exactly, thus I said this is rather a comment than an answer. However , whether there are more solutions are still unclear. –  Yimin Mar 26 '13 at 21:20
As far as for me - I like the nice try and would like if we could do more progress with an ansatz like this... (I've come across such an ansatz sometimes elsewhere and got stuck for the same reason of (principally) lack of accuracy. For instance in a treatize on the lucas-lehmer-test: go.helms-net.de/math/expdioph/lucasLehmer.pdf ). The additional argument of @Steven makes it even more complicated... –  Gottfried Helms Mar 28 '13 at 6:46