mersenneforum.org Mp mod D
 Register FAQ Search Today's Posts Mark Forums Read

 2022-07-16, 15:03 #1 User140242   Jul 2022 22·13 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,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?
 2022-07-16, 15:27 #2 Denial140   Dec 2021 29 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.
 2022-07-16, 15:59 #3 kriesel     "TF79LL86GIMPS96gpu17" Mar 2017 US midwest 22×52×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 232 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 2022-07-16 at 16:03
 2022-07-16, 16:07 #4 User140242   Jul 2022 22·13 Posts The language is python and // is integer division. The algorithm is the same as in the maths page but reversed.
 2022-07-16, 19:28 #5 User140242   Jul 2022 22·13 Posts 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 Last fiddled with by User140242 on 2022-07-16 at 19:31
 2022-07-16, 20:50 #6 Viliam Furik     "Viliam Furík" Jul 2018 Martin, Slovakia 24×72 Posts 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.
 2022-07-16, 22:15 #7 User140242   Jul 2022 22·13 Posts The proof is simple in general if Ri+1=(2k)%D Ri=(2k+1)%D=(2*Ri+1)%D = (2*Ri+1=D Ri=2*Ri+1-D or also Ri+1=(D+Ri)/2 but D is odd then Ri is odd then Ri+1=1+D//2+Ri//2 otherwise if 2*Ri+1
 2022-07-16, 22:42 #8 Batalov     "Serge" Mar 2008 Phi(4,2^7658614+1)/2 234018 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 log2(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.
 2022-07-16, 23:01 #9 User140242   Jul 2022 22·13 Posts 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
2022-07-16, 23:06   #10
Batalov

"Serge"
Mar 2008
Phi(4,2^7658614+1)/2

5×1,997 Posts

Quote:
 Originally Posted by User140242 there is no need to take 23 steps
Exactly! 5 steps is enough, if done properly.

 2022-07-17, 10:40 #11 User140242   Jul 2022 22×13 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 right-shift etc. Code: #include #include #include 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; } Last fiddled with by User140242 on 2022-07-17 at 10:46

All times are UTC. The time now is 16:02.

Sat Nov 26 16:02:13 UTC 2022 up 100 days, 13:30, 1 user, load averages: 1.49, 1.14, 1.03