Mp mod D

[User140242]

I don't know if it fits into the topic but some time ago I developed this algorithm to test if 2^Mp = 1 (mod D).

def is_pow_onemod(Mp,D):
    # return true if 2^Mp = 1 (mod D)
    R=1
    for i in range(0,Mp-len(bin(D)[2:])+1):
        if R%2==1:
           R=1+D//2+R//2
        else:
           R/=2
    if R==1<<(Mp-i-1):
        return True
    else:
        return False

Adequately adapted can't make it faster than mpz_pow?

[Denial140]

This algorithm is detailed in the "Trial Factoring" section of the maths page on mersenne.org. So there is no doubt this is already being used.

[kriesel]

Welcome to the forum. It's not clear to me what // as an operator means or what language that code is. (Please clarify.)

You may find it useful to consult the available reference info before posting. See for example https://www.mersenneforum.org/showpo...23&postcount=6 regarding GPU based trial factoring. GPUs are so much more efficient at TF, up to exponent 2³² and TF candidates up to 92 bits (mfaktc) or 95 (mfakto), it's a waste of CPU resources do do such TF on CPUs generally.

[User140242]

The language is python and // is integer division. The algorithm is the same as in the maths page but reversed.

[User140242]

example

is_pow_onemod(11,23) return True

in python bin(23)[2:]=10111 5 bit

number step 11-5+1=7

D//2=11

R=1

1 R=1+11=12
2 R=6
3 R=3
4 R=1+11+1=13
5 R=1+11+6=18
6 R=9
7 R=1+11+4=16


R==16=2^4 and 4=11-7 then return true

[Viliam Furik]

Well, it only uses division by 2, which can be done as a right-shift, and addition and value comparing, so if it works as intended, then I believe it is much faster - at least algorithmically, as it doesn't use multiplication, and barely uses a left-shift.

[User140242]

The proof is simple in general if

R_i+1=(2^k)%D
R_i=(2^k+1)%D=(2*R_i+1)%D = (2*R_i+1<D) ? 2*R_i+1 : 2*R_i+1-D

then if 2*R_i+1>=D R_i=2*R_i+1-D or also R_i+1=(D+R_i)/2 but D is odd then R_i is odd then R_i+1=1+D//2+R_i//2

otherwise if 2*R_i+1<D then R_i=2*R_i+1 or also R_i+1=R_i/2 then R_i is even

To prove that 2^P = 1 (mod D) we proceed in the opposite direction placing R = (2^P)% D = 1
then R is odd and R₁ = 1 + D // 2 + R // 2 then we find R₂ and so on


in the end since (2^b-1)%D = 2^b-1 with b = number of bits of D proceeding with the algorithm

if R_P-b+1 = 2^b-1 then 2^P=1 (mod D)

[Batalov]

You are computing Mod(1/2,D)^7 and compare to 2^4. (where 4 is simply a red herring value. You can change it to 3, to 7, to anything.)

That's the same as as computing Mod(1/2,D)^P by dividing by 2, P times.
This will be slower than doing squarings x log₂(P) times. Even for P=11.

Try your algorithm for p=86000063 and D=61920045361.
Compare to how appropriate algorithm works. How many steps? You will use > 86000000 steps, but only 27 steps are needed.

[User140242]

The algorithm is correct

D=47 b=6 bit

step 23-6+1=18

R=1

R1 =1+23=24
R2 =12
R3 =6
R4 =3
R5 =1+23+1=25
R6 =1+23+12=36
R7 =18
R8 =9
R9 =1+23+4=28
R10 =14
R11 =7
R12 =1+23+3=27
R13 =1+23+13=37
R14 =1+23+18=42
R15 =21
R16 =1+23+10=34
R17 =17
R18 =1+23+8=32=2^5=2^(b-1)


there is no need to take 23 steps

[Batalov]

Originally Posted by User140242

there is no need to take 23 steps

Exactly! 5 steps is enough, if done properly.

[User140242]

I quickly implemented this version in C++ and for p = 86000063 and D = 61920045361 even though I don't have a comparison with the other method it's fast.

Of course you have to optimize for example using right-shift etc.

#include <iostream>
#include <cmath>
#include <stdint.h>

bool is_pow_onemod(uint64_t P , uint64_t D)
{
    
    uint64_t  D_2=D/(uint64_t)2;
    int b = 0;
    uint64_t  N=D;
    while (N)
    {
        b++;
        N >>= 1;
    }
    uint64_t  imax=(uint64_t)P-b;
    uint64_t  R=(uint64_t)1;
    for (uint64_t i=0;i<=imax;i++)
        R=(R%2==1)?(uint64_t)1+R/(uint64_t)2+D_2:R/(uint64_t)2;
    if (R==(uint64_t)1 <<(b-1))
        return true;
    return false;
}