20170714, 09:25  #1 
Banned
"Luigi"
Aug 2002
Team Italia
3×5×317 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^%llu1)+1 is a factor of MM%d !!! \n", k, exponent, exponent); } else { gmp_printf("\n2*%llu*(2^%llu1)+1 is not a factor of MM%d ... \n", k, exponent, exponent); } I have been told that "For the DMers divisibility it is easier to use exponentiation of Mod(2, k*2p+1(2*k1)) to 2^p (no need to multby2 it each iteration) and compare the result to 2." Now, I read about sqrmodbased 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: 
20170714, 12:47  #2 
Sep 2002
Database er0rr
CBD_{16} Posts 
Adding one to MM will make it 2^p  you calculate it here:
Code:
mpz_pow_ui(MM, base, exponent); 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 20170714 at 13:02 
20170714, 13:37  #3  
Sep 2002
Database er0rr
3·1,087 Posts 
You might also look at:
Quote:


20170714, 14:14  #4 
Banned
"Luigi"
Aug 2002
Team Italia
3×5×317 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 Code:
"For the DMers divisibility it is easier to use exponentiation of Mod(2, k*2p+1(2*k1)) to 2^p (no need to multby2 it each iteration) and compare the result to 2." Now, I can check primality of such huge factors using pfgw's P1 test on PRP factors 2kp+1 where p is the size of a Mersenne prime > M(34), in other words p has more than 110 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 Top5000 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^expo1))+1 DIMS rstr OPENFILEAPP r_file,results.txt POWMOD result,base,2^expo1,mers WRITE r_file,result So, the final question is: Is it possible to write a GMP  gwlibbased  handtuned 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 
20170714, 15:20  #5 
"Robert Gerbicz"
Oct 2005
Hungary
2×3×227 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 <stdio.h> #include <time.h> #include "gmp.h" int dmm(mpz_ptr result,mpz_t k,unsigned int p){ // set result=2^(2^p) mod (k*(2^p1)+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==k1 modulo mod // result==a1*(k1)+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; } 
20170714, 16:15  #6 
Banned
"Luigi"
Aug 2002
Team Italia
3×5×317 Posts 

20170714, 20:30  #7 
Sep 2002
Database er0rr
CBD_{16} Posts 
This can be written in GWNUM using a special modulo...
We are calculating Mod(2^(2^p1)1, k*(2^p1)+1)==0 <=> Mod(2^(2^p1), k*(2^p1)+1)==1 <=> Mod(2, k*(2^p1)+1)^(2^p1)==1 <=> Mod(2, k*(2^p1)+1)^(2^p)==2 <=> Mod(2, k*2^p  (k1))^(2^p)==2 Thus we can use GWNUM's special modulo over k*2^p  (k1) 
20170714, 20:52  #8  
Banned
"Luigi"
Aug 2002
Team Italia
1293_{16} Posts 
Quote:
Unfortunately, I (still) didn't get acquainted with GWNUM. Maybe it's time to start learning something new... 

20170714, 20:58  #9 
Sep 2002
Database er0rr
3×1,087 Posts 
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.

20170714, 21:00  #10 
Banned
"Luigi"
Aug 2002
Team Italia
3·5·317 Posts 

20170714, 22:42  #11 
Sep 2002
Database er0rr
3·1,087 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, (k1)); 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 gwset_num_threads ( gwdata, 4 ); //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++ Last fiddled with by paulunderwood on 20170714 at 22:49 
Thread Tools  
Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Prime95 Source Code required  Qazi  PrimeNet  1  20170203 15:40 
Triple check required  houding  PrimeNet  14  20151221 09:34 
New hardware  help required  PageFault  Software  2  20140414 15:23 
New Algorithm  Tester(s) required  bearnol  Factoring  29  20060816 04:57 
Which version of SHLWAPI.dll is required?  gowen72  Software  1  20021101 00:36 