mersenneforum.org Help required
 Register FAQ Search Today's Posts Mark Forums Read

 2017-07-14, 09:25 #1 ET_ Banned     "Luigi" Aug 2002 Team Italia 476310 Posts Help required I wrote the following piece of GMP code to calculate the result of MMp mod kMp+1. All GMP values are of type mpz_t. All integer values are of type unsigned long int. The code works as expected, I left variables' definition and initializations outside the code for clarity. Code:  k2 = 2*k; // Mp definition mpz_pow_ui(MM, base, exponent); mpz_sub_ui(MM, MM, 1); // 2kMp+1 definition mpz_init_set(mod, MM); mpz_mul_ui(mod, mod, k2); mpz_add_ui(mod, mod, 1); // Calculate MMp mod kMp+1 mpz_powm(result, base, MM, mod); if (mpz_cmp_ui(result,1) == 0) { gmp_printf("2*%llu*(2^%llu-1)+1 is a factor of MM%d !!! \n", k, exponent, exponent); } else { gmp_printf("\n2*%llu*(2^%llu-1)+1 is not a factor of MM%d ... \n", k, exponent, exponent); } I need a hint on how to proceed to get a quicker, more optimized code. I have been told that "For the D-Mers divisibility it is easier to use exponentiation of Mod(2, k*2p+1-(2*k-1)) to 2^p (no need to mult-by-2 it each iteration) and compare the result to 2." Now, I read about sqrmod-based powmod, invmod, left/right side exponentiation, Montgomery multiplication/exponentiation, and am a bit worried and confused about the right method to use, and why. I am asking for some basic/technical help on how to proceed :smille:
 2017-07-14, 12:47 #2 paulunderwood     Sep 2002 Database er0rr 1101000010012 Posts Adding one to MM will make it 2^p -- you calculate it here: Code: mpz_pow_ui(MM, base, exponent); Create another variable MMp1. So you can replace the above with: Code: mpz_pow_ui(MMp1, base, exponent); mpz_sub_ui(MM, MMp1, 1); .... Code: mpz_powm(result, base, MMp1, mod); .... Code: if (mpz_cmp_ui(result,2) == 0) { Last fiddled with by paulunderwood on 2017-07-14 at 13:02
2017-07-14, 13:37   #3
paulunderwood

Sep 2002
Database er0rr

47×71 Posts

You might also look at:

Quote:
 Function: void mpz_mul_2exp (mpz_t rop, const mpz_t op1, mp_bitcnt_t op2) Set rop to op1 times 2 raised to op2. This operation can also be defined as a left shift by op2 bits.
Where op1 would be 1 in your case

 2017-07-14, 14:14 #4 ET_ Banned     "Luigi" Aug 2002 Team Italia 11×433 Posts I will try to be more specific. I know that Code: PARI/GP formula to check divisibility of a 2kp+1 factor for a factor Q. 2 * 584988 * (Mod(2, 100021968071941 )^74207281 - 1) + 1 == 0 2 * k * (Mod(2, Q)^p - 1) + 1 == 0 It should translate in Code: "For the D-Mers divisibility it is easier to use exponentiation of Mod(2, k*2p+1-(2*k-1)) to 2^p (no need to mult-by-2 it each iteration) and compare the result to 2." Am I right up to here? Now, I can check primality of such huge factors using pfgw's P-1 test on PRP factors 2kp+1 where p is the size of a Mersenne prime > M(34), in other words p has more than 1-10 million digits. Once the factor is a proven prime, I use the code on post #1 to test if the prime factor divides the double Mersenne number MMp. If it does, we have a new composite double Merssenne number, and a possible prime factor that will get into the Top-5000 prime database on the first 10 positions. The GMP code on post #1 works fine, at least with the known factors of double Mersennes. But it is also terribly slow. I used the following pfgw script to test MMp dividibility: Code:  SCRIPT DIM expo, 110503 DIM base, 2 DIM k, 1567437 DIM result, 0 DIM mers, (2*k*(2^expo-1))+1 DIMS rstr OPENFILEAPP r_file,results.txt POWMOD result,base,2^expo-1,mers WRITE r_file,result Sadly, such script stops working when mers becomes "too big", and writes a false positive inside the result file. My GMP code works fine, but it is quite slower. So, the final question is: Is it possible to write a GMP || gwlib-based || hand-tuned code that is both reliable and faster than my GMP code? And if it is, how? I am sorry to say that I didn't understand how the hints showed by Paul could possibly speed up the code
 2017-07-14, 15:20 #5 R. Gerbicz     "Robert Gerbicz" Oct 2005 Hungary 7·197 Posts Roughly 3 times faster than your routine, here you can even use k value that is larger than 2^64. Try it out for example your 3134874 110503 case (here 2*1567437=3134874). Code: #include #include #include "gmp.h" int dmm(mpz_ptr result,mpz_t k,unsigned int p){ // set result=2^(2^p) mod (k*(2^p-1)+1) // return 1 for a divisor unsigned int i,j; mpz_t lo,hi,a1,a0,k1,mod; mpz_init(lo); mpz_init(hi); mpz_init(a1); mpz_init(a0); mpz_init(k1); mpz_init(mod); mpz_sub_ui(k1,k,1); mpz_mul_2exp(mod,k,p); mpz_sub(mod,mod,k1); mpz_set_ui(result,2); for(i=1;i<=p;i++){ mpz_pow_ui(result,result,2); for(;;){ mpz_fdiv_q_2exp(hi,result,p); mpz_fdiv_r_2exp(lo,result,p); // result=hi*2^p+lo=(a1*k+a0)*2^p+lo mpz_fdiv_qr(a1,a0,hi,k); if(mpz_cmp_ui(a1,1)<=0)break; // k*2^p==k-1 modulo mod // result==a1*(k-1)+a0*2^p+lo mpz_mul_2exp(result,a0,p); mpz_addmul(result,a1,k1); mpz_add(result,result,lo); if(mpz_cmp(result,mod)<0)break; } } mpz_mod(result,result,mod); mpz_clear(lo); mpz_clear(hi); mpz_clear(a1); mpz_clear(a0); mpz_clear(k1); mpz_clear(mod); return (mpz_cmp_ui(result,2)==0); } int main(void){ int ans; unsigned int p; mpz_t r,k; mpz_init(k); mpz_init(r); while(gmp_scanf("%Zd%u",k,&p)!=EOF){ time_t sec=time(NULL); ans=dmm(r,k,p); gmp_printf("result=%Zd;ans=%d;time=%ld sec.\n",r,ans,time(NULL)-sec); } mpz_clear(k); mpz_clear(r); return 0; }
2017-07-14, 16:15   #6
ET_
Banned

"Luigi"
Aug 2002
Team Italia

112338 Posts

Quote:
 Originally Posted by R. Gerbicz Roughly 3 times faster than your routine, here you can even use k value that is larger than 2^64. Try it out for example your 3134874 110503 case (here 2*1567437=3134874).
Thank you Robert

Now I have something to study for the next 3 weeks.

 2017-07-14, 20:30 #7 paulunderwood     Sep 2002 Database er0rr 47·71 Posts This can be written in GWNUM using a special modulo... We are calculating Mod(2^(2^p-1)-1, k*(2^p-1)+1)==0 <=> Mod(2^(2^p-1), k*(2^p-1)+1)==1 <=> Mod(2, k*(2^p-1)+1)^(2^p-1)==1 <=> Mod(2, k*(2^p-1)+1)^(2^p)==2 <=> Mod(2, k*2^p - (k-1))^(2^p)==2 Thus we can use GWNUM's special modulo over k*2^p - (k-1)
2017-07-14, 20:52   #8
ET_
Banned

"Luigi"
Aug 2002
Team Italia

11·433 Posts

Quote:
 Originally Posted by paulunderwood This can be written in GWNUM using a special modulo... We are calculating Mod(2^(2^p-1)-1, k*(2^p-1)+1)==0 <=> Mod(2^(2^p-1), k*(2^p-1)+1)==1 <=> Mod(2, k*(2^p-1)+1)^(2^p-1)==1 <=> Mod(2, k*(2^p-1)+1)^(2^p)==2 <=> Mod(2, k*2^p - (k-1))^(2^p)==2 Thus we can use GWNUM's special modulo over k*2^p - (k-1)
Now, this is what I needed to fully understand the theory behind the formula, thank you Paul!

Unfortunately, I (still) didn't get acquainted with GWNUM. Maybe it's time to start learning something new...

2017-07-14, 20:58   #9
paulunderwood

Sep 2002
Database er0rr

47·71 Posts

Quote:
 Originally Posted by ET_ Now, this is what I needed to fully understand the theory behind the formula, thank you Paul! Unfortunately, I (still) didn't get acquainted with GWNUM. Maybe it's time to start learning something new...
You are welcome. I can knock something together over the weekend -- I will then leave it to you to do the I/O and interface.

2017-07-14, 21:00   #10
ET_
Banned

"Luigi"
Aug 2002
Team Italia

11×433 Posts

Quote:
 Originally Posted by paulunderwood You are welcome. I can knock something together over the weekend -- I will then leave it to you to do the I/O and interface.
You are too kind
Meanwhile I will try to download and compile the gwnum library as a diligent student.

2017-07-14, 22:42   #11
paulunderwood

Sep 2002
Database er0rr

333710 Posts

Here you go:

Code:
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "gwnum/giants.h"
#include "gwnum/gwnum.h"
#include "gwnum/gwcommon.h"

unsigned init_wa[1] = {2};
gwnum  wa;
gwhandle *gwdata;

int is_a_divisor ( unsigned long p, unsigned long k ) {
int i;
//setup mod
gwsetup ( gwdata, k, 2, p, -(k-1));
wa = gwalloc ( gwdata );
// initiaize wa to 2
binarytogw ( gwdata, init_wa, 1, wa );
//do some slow rounds
gwset_square_carefully_count ( gwdata, 30 );
// main loop
gwstartnextfft ( gwdata, 1 );
for( i = 1; i < p; i++ ) gwsquare ( gwdata, wa );
gwstartnextfft ( gwdata, 0 );
gwsquare ( gwdata, wa );
// return
gwsmalladd ( gwdata, -2.0, wa );
i = gwiszero ( gwdata, wa );
gwdone( gwdata );
return ( i );
}
// MAIN
int main ( int argc, char *argv[] ) {
if( argc != 3 ) {
printf( "Usage: dMM p k\n" );
exit ( -1 );
}
// set up gwdata
gwdata = (gwhandle*) malloc (sizeof ( gwhandle ) );
gwinit ( gwdata );
// 4 threads - could be an ARGV
//assume k is even!!
if ( is_a_divisor ( atoll ( argv[1] ), atoll ( argv[2] ) ) )
printf( "a divisor!\n");
else
printf( "does not divide!\n");
}
//gcc -o dMM dMM.c gwnum/gwnum.a gwnum/gwnum.ld -lm -lpthread -lstdc++
Attached Files
 dMM.c.zip (783 Bytes, 36 views)

Last fiddled with by paulunderwood on 2017-07-14 at 22:49

 Similar Threads Thread Thread Starter Forum Replies Last Post Qazi PrimeNet 1 2017-02-03 15:40 houding PrimeNet 14 2015-12-21 09:34 PageFault Software 2 2014-04-14 15:23 bearnol Factoring 29 2006-08-16 04:57 gowen72 Software 1 2002-11-01 00:36

All times are UTC. The time now is 00:46.

Thu Aug 13 00:46:26 UTC 2020 up 26 days, 20:33, 0 users, load averages: 1.80, 1.64, 1.53