20220716, 15:03  #1 
Jul 2022
49_{10} 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
5^{2} 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^{5}×3×71 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
61_{8} 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
7^{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
19·41 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
31_{16} 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
9936_{10} 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
7^{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
26D0_{16} Posts 

20220717, 10:40  #11 
Jul 2022
7^{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 