20220716, 15:03  #1 
Jul 2022
3×5^{2} Posts 
Mp mod D
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).
Code:
def is_pow_onemod(Mp,D): # return true if 2^Mp = 1 (mod D) R=1 for i in range(0,Mplen(bin(D)[2:])+1): if R%2==1: R=1+D//2+R//2 else: R/=2 if R==1<<(Mpi1): return True else: return False 
20220716, 15:27  #2 
Dec 2021
43 Posts 
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.

20220716, 15:59  #3 
"TF79LL86GIMPS96gpu17"
Mar 2017
US midwest
2·29·127 Posts 
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^{32} 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. Last fiddled with by kriesel on 20220716 at 16:03 
20220716, 16:07  #4 
Jul 2022
1001011_{2} Posts 
The language is python and // is integer division. The algorithm is the same as in the maths page but reversed.

20220716, 19:28  #5 
Jul 2022
3×5^{2} Posts 
example
is_pow_onemod(11,23) return True in python bin(23)[2:]=10111 5 bit number step 115+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=117 then return true Last fiddled with by User140242 on 20220716 at 19:31 
20220716, 20:50  #6 
"Viliam Furík"
Jul 2018
Martin, Slovakia
13·61 Posts 
Well, it only uses division by 2, which can be done as a rightshift, 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 leftshift.

20220716, 22:15  #7 
Jul 2022
3·5^{2} Posts 
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} = 1 + D // 2 + R // 2 then we find R_{2} and so on in the end since (2^{b1})%D = 2^{b1} with b = number of bits of D proceeding with the algorithm if R_{Pb+1} = 2^{b1} then 2^{P}=1 (mod D) Last fiddled with by User140242 on 20220716 at 22:25 
20220716, 22:42  #8 
"Serge"
Mar 2008
Phi(4,2^7658614+1)/2
2^{2}·5·503 Posts 
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_{2}(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. 
20220716, 23:01  #9 
Jul 2022
3·5^{2} Posts 
The algorithm is correct
D=47 b=6 bit step 236+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^(b1) there is no need to take 23 steps 
20220716, 23:06  #10 
"Serge"
Mar 2008
Phi(4,2^7658614+1)/2
2^{2}·5·503 Posts 

20220717, 10:40  #11 
Jul 2022
3·5^{2} Posts 
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 rightshift etc. Code:
#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)Pb; 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 <<(b1)) return true; return false; } Last fiddled with by User140242 on 20220717 at 10:46 