Mersenne Semiprimes

[ET_]

Originally Posted by Batalov

With a small code patch (like the one described earlier), we can use AVX or FMA3 FFT length 288K PRP test. Much faster. It only takes a couple hours on 2-4 threads.

This is because the tested value is a Phi(p*q, 2) where 2^p-1 is prime; earlier I used this code for Phi(3*p, -b), and now we can use it for any Phi(p*q, 2) including Phi(p^2, 2).

Is a similar patch also appliable to double Mersenne factors?

[Batalov]

Originally Posted by ET_

Is a similar patch also applicable to double Mersenne factors?

Double Mersenne factor checks already use the optimal (not general) FFT. It is exactly as if it were a 2-PRP test which is cut short (less iterations: 2^p-1 instead of 2*k*(2^p-1)+1).

When we use gwnum-based tools to run ECM, P-1 or PRP on a (k,b,n,c) number divided by some reasonably small value e (which is provided in "..." field), then all calculations can be done with special mod on (k,b,n,c) basis. Only the last step requires general mod N and then we have the final residue.

In contrast, when we are considering a (k,b,n,c) number divided by something large, then all gwnum based tools by default prepare the N value and then treat it as a general number - which leads to ~2 times larger zero-padded FFT size and using general mod step. So, that's ~2x slower due to FFT size and additional ~2x slower due to general mod. The patch ad hoc reduces this situation to the one above.

[Batalov]

Originally Posted by Batalov

(less iterations: 2^p-1 instead of 2*k*(2^p-1)+1).

Sloppy phrasing there. What I meant was of course log₂(2^p-1), or simply p steps, instead of log₂(2*k*(2^p-1)+1) steps with an appropriate bit pattern from N.
For the D-Mers divisibility it is easier to use exponentiation of Mod(2, k*2^p+1-(2*k-1)) to 2^p (no need to mult-by-2 it each iteration) and compare the result to 2.

[JeppeSN]

Originally Posted by Batalov

[Wed Feb  3 01:19:54 2016 UTC]
M5202961/M2281 is not prime.  RES64: 2D2764F7AB292478. We8: 00000000,00000000
[Wed Feb  3 04:36:31 2016 UTC]
M19562929/M4423 is not prime.  RES64: 83DA690553D9EE3C. We8: 00000000,00000000
...
[Work thread Feb 3 05:02] Starting PRP test of M98823481/M9941 using FMA3 FFT length 5376K, Pass1=896, Pass2=6K, 18 threads
...
[Work thread Feb 3 05:47] Iteration: 200000 / 98813540 [0.20%], ms/iter:  2.215, ETA: 60:41:15
[Work thread Feb 3 05:48] Iteration: 210000 / 98813540 [0.21%], ms/iter:  2.216, ETA: 60:41:25

That's not bad for a "M98823481" candidate size! I am going to run these last two, too.

Serge, that is great. You are going to gain instant fame if $\Phi_{p^2}(2) = M_{p^2}/M_p$ is prime for either $p=9941$ or $p=11213$ /JeppeSN

[LaurV]

He can at most find a very large PRP. It doesn't seem to me like a real primality test. And he is already famous! (Respect! )

[Batalov]

Yeah, a large PRP is nice (with almost negligible probability though) ...but a large prime that I'd really like to find would be a b^2[SUP]n+1[/SUP] - b^2[SUP]n[/SUP] +1 which is provable. Or any other Phi(3ⁿ*2^m, b) (a bundle of 3's still leads to a B^2 - B + 1 form -- which is {evidently from its appearance} provable by N-1 method with 150% proof), known as P.I.E.S. primes.

Phi(Kⁿ*2^m, b), even if found, would also remain a PRP (with K>=5). Same with any other higher cyclotomic.

[axn]

Originally Posted by Batalov

Phi(Kⁿ*2^m, b), even if found, would also remain a PRP (with K>=5). Same with any other higher cyclotomic.

Isn't it possible to do a CHG proof for K=5, since you get 25% factorization of N-1 for free (plus algebraic factors of the remainder from which some ECM factorizations might also be obtainable)?

[Batalov]

Small ones could be proven but not the 388340+ digit numbers (the UTM cutoff).
Even to get to 26% we need another 3883 digits of factors, but a 26% proof would be nearly impossible at that size. (I did do a 29.08% proof once ...at that size. I couldn't believe my eyes when the first step of the proof came through overnight. I think a high 28-ish% proof could work -- but that's more than 10000 digits of extra factors to collect, almost impossible.)

[Prime95]

Originally Posted by Batalov

Here is the modified patch:

Patch added to prime95 source code. You'll need to put "PhiExtensions=1" in prime.txt.

[Batalov]

I'll clean it up some more.
Now that I've used it in action for a while, I noticed that divg() might be not the best way to do it (even though simple to write); it uses quite a bit of RAM which many users might not have. For M98823481/M9941, >55G is used at peak. Then the memory use wanes down for the rest of the run, so I gather that this is the initial divg(). Instead I can code the summation of (-b)^kq, k=0..p-1.

Also, of course, it is only for b=2, for two reasons: a) the array of known M-primes should not be used, and b) the denominator must be itself divided by (b-1) for the second division.

[Batalov]

[Work thread Feb 6 00:44] 
M98823481/"1" is not prime.  RES64: 333069AAE77863BA. We8: 00000000,00000000
...
[Work thread Feb 6 00:44] 
Starting PRP test of M125731369/"1" using FMA3 FFT length 6720K, Pass1=896, Pass2=7680, 18 threads
Killed

So, there, I precog'd the future out of memory death (the node has 64Gb); I did get lucky with divg() still fitting in memory for the previous candidate, but this one didn't fit. So, just in time to rewrite the code with summation; here, b=2, all terms with the plus sign, so I can code it easily from the small terms up. Ok, then...
EDIT: this patch-2.0 works much better:

        *N = allocgiant ((bits >> 5) + 5);
        if (*N == NULL) return (OutOfMemory (thread_num));

        if(w->k == 1.0 && w->b == 2 && w->c == -1) { /*=== this input means Phi(n,2) with n semiprime ===*/
            int i,k,q, knownSmallMers[] = {3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217,
               4253, 4423, 9689, 9941, 11213, 19937, 999999999}; /* for now, just the cases where w->n = p * q, and 2^q-1 is prime */
            for(i=0; (q=knownSmallMers[i]) < w->n || q*q <= w->n;i++) if((w->n%q) == 0) {
            giant tmp = allocgiant ((bits >> 5) + 5);
            if (!tmp) return (OutOfMemory (thread_num));
            ultog (1, tmp);
            ultog (1, *N);
            gshiftleft (w->n-w->n/q, *N);
            for(k=2; k < q; k++) {
              gshiftleft (w->n/q, tmp);
              addg (tmp, *N);
            }
            iaddg (1, *N);
            if(q != w->n/q) {
              ultog (w->b, tmp);
              power (tmp, w->n/q);
              iaddg (w->c, tmp);
              divg (tmp, *N);
            }
            free (tmp);
            /* w->forced_fftlen = w->n; */ /*=== too late to do this here. Moved before gwsetup() --SB. */
            if(!w->known_factors || !strcmp(w->known_factors,"1")) {
                  p = sprintf(buf, "M%d", w->n/q); if(q != w->n/q) p += sprintf(buf+p, "/M%d", q);
                  w->known_factors = malloc(p+1);
                  memcpy(w->known_factors, buf, p+1);
            }
            return (0);
          }
        }

        ultog (w->b, *N);