mersenneforum.org  

Go Back   mersenneforum.org > Extra Stuff > Miscellaneous Math

Reply
 
Thread Tools
Old 2022-07-16, 15:03   #1
User140242
 
Jul 2022

22·13 Posts
Default 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?
User140242 is offline   Reply With Quote
Old 2022-07-16, 15:27   #2
Denial140
 
Dec 2021

52 Posts
Default

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.
Denial140 is offline   Reply With Quote
Old 2022-07-16, 15:59   #3
kriesel
 
kriesel's Avatar
 
"TF79LL86GIMPS96gpu17"
Mar 2017
US midwest

22×5×73 Posts
Default

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
kriesel is online now   Reply With Quote
Old 2022-07-16, 16:07   #4
User140242
 
Jul 2022

22·13 Posts
Default

The language is python and // is integer division. The algorithm is the same as in the maths page but reversed.
User140242 is offline   Reply With Quote
Old 2022-07-16, 19:28   #5
User140242
 
Jul 2022

1101002 Posts
Default

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
User140242 is offline   Reply With Quote
Old 2022-07-16, 20:50   #6
Viliam Furik
 
Viliam Furik's Avatar
 
"Viliam Furík"
Jul 2018
Martin, Slovakia

19×41 Posts
Default

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.
Viliam Furik is offline   Reply With Quote
Old 2022-07-16, 22:15   #7
User140242
 
Jul 2022

22·13 Posts
Default

The proof is simple in general if

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

then if 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<D then Ri=2*Ri+1 or also Ri+1=Ri/2 then Ri is even

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


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

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

Last fiddled with by User140242 on 2022-07-16 at 22:25
User140242 is offline   Reply With Quote
Old 2022-07-16, 22:42   #8
Batalov
 
Batalov's Avatar
 
"Serge"
Mar 2008
Phi(4,2^7658614+1)/2

2×3×1,657 Posts
Default

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.
Batalov is offline   Reply With Quote
Old 2022-07-16, 23:01   #9
User140242
 
Jul 2022

22×13 Posts
Default

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
User140242 is offline   Reply With Quote
Old 2022-07-16, 23:06   #10
Batalov
 
Batalov's Avatar
 
"Serge"
Mar 2008
Phi(4,2^7658614+1)/2

233268 Posts
Default

Quote:
Originally Posted by User140242 View Post
there is no need to take 23 steps
Exactly! 5 steps is enough, if done properly.
Batalov is offline   Reply With Quote
Old 2022-07-17, 10:40   #11
User140242
 
Jul 2022

3416 Posts
Default

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 <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;
}

Last fiddled with by User140242 on 2022-07-17 at 10:46
User140242 is offline   Reply With Quote
Reply

Thread Tools


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


Mon Oct 3 15:16:29 UTC 2022 up 46 days, 12:45, 1 user, load averages: 1.48, 1.47, 1.44

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2022, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.

≠ ± ∓ ÷ × · − √ ‰ ⊗ ⊕ ⊖ ⊘ ⊙ ≤ ≥ ≦ ≧ ≨ ≩ ≺ ≻ ≼ ≽ ⊏ ⊐ ⊑ ⊒ ² ³ °
∠ ∟ ° ≅ ~ ‖ ⟂ ⫛
≡ ≜ ≈ ∝ ∞ ≪ ≫ ⌊⌋ ⌈⌉ ∘ ∏ ∐ ∑ ∧ ∨ ∩ ∪ ⨀ ⊕ ⊗ 𝖕 𝖖 𝖗 ⊲ ⊳
∅ ∖ ∁ ↦ ↣ ∩ ∪ ⊆ ⊂ ⊄ ⊊ ⊇ ⊃ ⊅ ⊋ ⊖ ∈ ∉ ∋ ∌ ℕ ℤ ℚ ℝ ℂ ℵ ℶ ℷ ℸ 𝓟
¬ ∨ ∧ ⊕ → ← ⇒ ⇐ ⇔ ∀ ∃ ∄ ∴ ∵ ⊤ ⊥ ⊢ ⊨ ⫤ ⊣ … ⋯ ⋮ ⋰ ⋱
∫ ∬ ∭ ∮ ∯ ∰ ∇ ∆ δ ∂ ℱ ℒ ℓ
𝛢𝛼 𝛣𝛽 𝛤𝛾 𝛥𝛿 𝛦𝜀𝜖 𝛧𝜁 𝛨𝜂 𝛩𝜃𝜗 𝛪𝜄 𝛫𝜅 𝛬𝜆 𝛭𝜇 𝛮𝜈 𝛯𝜉 𝛰𝜊 𝛱𝜋 𝛲𝜌 𝛴𝜎𝜍 𝛵𝜏 𝛶𝜐 𝛷𝜙𝜑 𝛸𝜒 𝛹𝜓 𝛺𝜔