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

2022-07-20, 23:20   #24
chalsall
If I May

"Chris Halsall"
Sep 2002

24×5×7×19 Posts

Quote:
 Originally Posted by Batalov Just 2 cents for the road. These hints might be useful for the general process of prototyping anything. Anything that you might think of tomorrow.
Sound counsel.

2022-07-21, 18:30   #25
User140242

Jul 2022

72 Posts

Quote:
 Originally Posted by Batalov Just to give you, perhaps, a sandwich for the road (because this road is long), there are certain elements that can be used for 1) prototyping, 2) refactoring into something intermediate, 3) deep refactoring, with custom building blocks.

In the meantime I created this function that allows you to find the factor (~ 42 bits) of 13379827 or the factor (~ 46 bits) of 668267879 in less than 0.7 s.

Code:
uint64_t XY_modulo_D_LN(uint64_t X,uint64_t Y,uint64_t D,uint64_t D_max)
{
// D<2^54 D>D_max
uint64_t mod_D;
if (X<=D_max && Y<=D_max )
{
mod_D=(X*Y)%D;
}
else if (X<=D_max || Y<=D_max )
{
uint64_t X1,N_t;
uint64_t RD=(uint64_t)D%D_max;
uint64_t KD=(uint64_t)D/D_max;
if (Y<=D_max)
{
X1=Y;
Y=X;
X=X1;
}
mod_D=(X*(Y%D_max))%D_max;
N_t=X*(Y/D_max)+(X*(Y%D_max))/D_max;
mod_D=(mod_D+(D_max*(N_t%KD))%D)%D;
mod_D=(mod_D+D-((N_t/KD)*RD)%D)%D;

}
else
{
uint64_t RD=(uint64_t)D%D_max;
uint64_t KD=(uint64_t)D/D_max;
uint64_t Rx=(uint64_t)X%D_max;
uint64_t Ry=(uint64_t)Y%D_max;
uint64_t Kx=(uint64_t)X/D_max;
uint64_t Ky=(uint64_t)Y/D_max;
uint64_t mod_D_t,N_t,R_t;

mod_D=(Rx*Ry)%D;

mod_D_t=((D_max*(D_max%KD))%D)%D;
mod_D_t=(mod_D_t+D-((D_max/KD)*RD)%D)%D;

N_t=mod_D_t;
R_t=Kx;
if (N_t>=(D_max/R_t))
{
mod_D_t=(R_t*(N_t%D_max))%D_max;
N_t=R_t*(N_t/D_max)+(R_t*(N_t%D_max))/D_max;
mod_D_t=(mod_D_t+(D_max*(N_t%KD))%D)%D;
mod_D_t=(mod_D_t+D-((N_t/KD)*RD)%D)%D;
}
else
mod_D_t=mod_D_t*R_t;
N_t=mod_D_t;
R_t=Ky;
if (N_t>=(D_max/R_t))
{
mod_D_t=(R_t*(N_t%D_max))%D_max;
N_t=R_t*(N_t/D_max)+(R_t*(N_t%D_max))/D_max;
mod_D_t=(mod_D_t+(D_max*(N_t%KD))%D)%D;
mod_D_t=(mod_D_t+D-((N_t/KD)*RD)%D)%D;
}
else
mod_D_t=mod_D_t*R_t;
N_t=mod_D+mod_D_t;
mod_D=N_t;
if (N_t>=D_max)
{
mod_D=(N_t)%D_max;
mod_D=(mod_D+(D_max*((N_t/D_max)%KD))%D)%D;
mod_D=(mod_D+D-(((N_t/D_max)/KD)*RD)%D)%D;
}

if (Rx>=KD)
{
mod_D_t=((D_max*(Rx%KD))%D)%D;
mod_D_t=(mod_D_t+D-((Rx/KD)*RD)%D)%D;
}
else
mod_D_t=(Rx*D_max)%D;
N_t=mod_D_t;
R_t=Ky;
if (N_t>=(D_max/R_t))
{
mod_D_t=(R_t*(N_t%D_max))%D_max;
N_t=R_t*(N_t/D_max)+(R_t*(N_t%D_max))/D_max;
mod_D_t=(mod_D_t+(D_max*(N_t%KD))%D)%D;
mod_D_t=(mod_D_t+D-((N_t/KD)*RD)%D)%D;
}
else
mod_D_t=mod_D_t*R_t;
N_t=mod_D+mod_D_t;
mod_D=N_t;
if (N_t>=D_max)
{
mod_D=(N_t)%D_max;
mod_D=(mod_D+(D_max*((N_t/D_max)%KD))%D)%D;
mod_D=(mod_D+D-(((N_t/D_max)/KD)*RD)%D)%D;
}

if (Ry>=KD)
{
mod_D_t=((D_max*(Ry%KD))%D)%D;
mod_D_t=(mod_D_t+D-((Ry/KD)*RD)%D)%D;
}
else
mod_D_t=(Ry*D_max)%D;
N_t=mod_D_t;
R_t=Kx;
if (N_t>=(D_max/R_t))
{
mod_D_t=(R_t*(N_t%D_max))%D_max;
N_t=R_t*(N_t/D_max)+(R_t*(N_t%D_max))/D_max;
mod_D_t=(mod_D_t+(D_max*(N_t%KD))%D)%D;
mod_D_t=(mod_D_t+D-((N_t/KD)*RD)%D)%D;
}
else
mod_D_t=mod_D_t*R_t;
N_t=mod_D+mod_D_t;
mod_D=N_t;
if (N_t>=D_max)
{
mod_D=(N_t)%D_max;
mod_D=(mod_D+(D_max*((N_t/D_max)%KD))%D)%D;
mod_D=(mod_D+D-(((N_t/D_max)/KD)*RD)%D)%D;
}
}

return mod_D;
}

bool is_2pow_onemod(uint64_t P,uint64_t D)
{
// return true if 2^P = 1 (mod D)

uint64_t x=2;
uint64_t y=1;
uint64_t D_max=(uint64_t)1<<30;
if (D<=D_max)
{
while (P > (uint64_t)1)
{
if (P%2 == 1)
y = (x * y) %D;
x = (x * x) %D;
P >>= 1;
}
return (((x*y)%D) == 1);
}
else
{
while (P > (uint64_t)1)
{
if (P%2 == 1)
y=XY_modulo_D_LN(x,y,D,D_max);
x=XY_modulo_D_LN(x,x,D,D_max);
P >>= 1;
}
return (XY_modulo_D_LN(x,y,D,D_max) == (uint64_t)1);
}

}

2022-07-21, 19:23   #26
science_man_88

"Forget I exist"
Jul 2009
Dumbassville

24·3·52·7 Posts

Quote:
 Originally Posted by User140242 I will try to follow your advice. In the meantime I created this function that allows you to find the factor (~ 42 bits) of 13379827 or the factor (~ 46 bits) of 668267879 in less than 0.7
I'm not fluent in C but did your code use Mp prime ? If so D=2k(Mp)+1 for some k. You only have to test D with certain k values. Nevermind I see the 1 mod 8 versus 7 mod 8 talk earlier in the thread didn't read it fully.

Last fiddled with by science_man_88 on 2022-07-21 at 20:00

2022-07-21, 21:17   #27
science_man_88

"Forget I exist"
Jul 2009
Dumbassville

840010 Posts

Quote:
 Originally Posted by science_man_88 I'm not fluent in C but did your code use Mp prime ? If so D=2k(Mp)+1 for some k. You only have to test D with certain k values. Nevermind I see the 1 mod 8 versus 7 mod 8 talk earlier in the thread didn't read it fully.
you can check $k\not\equiv p\pmod 3$ because that's when 2kp+1 doesn't divide by 3.

2022-07-22, 08:21   #28
User140242

Jul 2022

72 Posts

Quote:
 Originally Posted by science_man_88 you can check $k\not\equiv p\pmod 3$ because that's when 2kp+1 doesn't divide by 3.

Here you find the final version https://gist.github.com/user140242/0...d2738788bb1cae.

A wheel sieve is used which has been set to eliminate all multiples of small primes less than min(bW, 524288) including multiples of 3.

2022-09-04, 15:44   #29
User140242

Jul 2022

1100012 Posts

To get the mind busy.

Several times I have been told in this thread that the best solution to calculate (2 ^ p-1) module D is the one described in the math page.

As described, given a certain number P, a number of steps equal to the number of bits of P converted into binary are required.

Example for p=86000063 you need 27 steps.

In reality, different solutions can be adopted to achieve the same result with the same complexity, for example:

fix b = 6 then 2^b=64

P=P%64 +P/64 *64 = P%64 +b1 *64

b1=b1%64 +b1/64 *64 = P%64 +b2 *64 then b2=b1/64

.....

bm+1 =bm%64 +bm/64 *64

continue until bm/64 <64

the remainder must be stored in a vector [P%64, b1%64, ... , bm%64]

then we start from

mod=2^(bm/64)

the product is repeated for b=6 times

mod=(mod*mod)%D
mod=(mod*mod)%D
mod=(mod*mod)%D
mod=(mod*mod)%D
mod=(mod*mod)%D
mod=(mod*mod)%D

then it is calculated

mod=(mod* 2^(bm%64) )%D

the whole is repeated for all the remainders in the vector.

Practical example

P=86000063 and D=61920045361

Quote:
 P = 86000063 = 63 + 1343750 *64 b1 = 1343750 = 6 + 20996 *64 b2 = 20996 = 4 + 328 *64 b3 = 328 = 8 + 5 *64 remainders vector [63, 6, 4, 8] b3 / 64 = 5
then

Quote:
 mod=2^5=32 --- mod=(mod*mod)%D = 1024 mod=(mod*mod)%D = 1048576 mod=(mod*mod)%D = 46870856639 mod=(mod*mod)%D = 2154143653 mod=(mod*mod)%D = 52612687576 mod=(mod*mod)%D = 22360353543 mod=(mod* 2^8 )%D = 27606333796 --- mod=(mod*mod)%D = 4234002796 mod=(mod*mod)%D = 34887580973 mod=(mod*mod)%D = 39592247589 mod=(mod*mod)%D = 51169879776 mod=(mod*mod)%D = 26016070733 mod=(mod*mod)%D = 29657271891 mod=(mod* 2^4 )%D = 41076032729 --- mod=(mod*mod)%D = 22120162447 mod=(mod*mod)%D = 45609483711 mod=(mod*mod)%D = 49416360515 mod=(mod*mod)%D = 58025736930 mod=(mod*mod)%D = 25169243724 mod=(mod*mod)%D = 29919015099 mod=(mod* 2^6 )%D = 57215605506 --- mod=(mod*mod)%D = 31050105354 mod=(mod*mod)%D = 2960155768 mod=(mod*mod)%D = 17389932407 mod=(mod*mod)%D = 28504056210 mod=(mod*mod)%D = 40422916975 mod=(mod*mod)%D = 9761495609 mod=(mod* 2^63 )%D = 1
in this way 28 steps are carried out but it is not necessary to distinguish the two cases on the basis of the value of the bits of P in each single step.

 2022-09-15, 13:21 #30 User140242   Jul 2022 72 Posts In reference to post #29 maybe this: Find factors of a Mersenne number might be of interest to someone.

All times are UTC. The time now is 23:08.

Sun Sep 25 23:08:55 UTC 2022 up 38 days, 20:37, 0 users, load averages: 1.49, 1.22, 1.16