mersenneforum.org  

Go Back   mersenneforum.org > Extra Stuff > Programming

Reply
 
Thread Tools
Old 2017-07-14, 09:25   #1
ET_
Banned
 
ET_'s Avatar
 
"Luigi"
Aug 2002
Team Italia

3×5×317 Posts
Exclamation 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:
ET_ is offline   Reply With Quote
Old 2017-07-14, 12:47   #2
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

CBD16 Posts
Default

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
paulunderwood is offline   Reply With Quote
Old 2017-07-14, 13:37   #3
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

3·1,087 Posts
Default

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
paulunderwood is offline   Reply With Quote
Old 2017-07-14, 14:14   #4
ET_
Banned
 
ET_'s Avatar
 
"Luigi"
Aug 2002
Team Italia

3×5×317 Posts
Default

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
ET_ is offline   Reply With Quote
Old 2017-07-14, 15:20   #5
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

2×3×227 Posts
Default

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^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;
}
R. Gerbicz is offline   Reply With Quote
Old 2017-07-14, 16:15   #6
ET_
Banned
 
ET_'s Avatar
 
"Luigi"
Aug 2002
Team Italia

3×5×317 Posts
Default

Quote:
Originally Posted by R. Gerbicz View Post
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.
ET_ is offline   Reply With Quote
Old 2017-07-14, 20:30   #7
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

CBD16 Posts
Default

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)
paulunderwood is offline   Reply With Quote
Old 2017-07-14, 20:52   #8
ET_
Banned
 
ET_'s Avatar
 
"Luigi"
Aug 2002
Team Italia

129316 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
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...
ET_ is offline   Reply With Quote
Old 2017-07-14, 20:58   #9
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

3×1,087 Posts
Default

Quote:
Originally Posted by ET_ View Post
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.
paulunderwood is offline   Reply With Quote
Old 2017-07-14, 21:00   #10
ET_
Banned
 
ET_'s Avatar
 
"Luigi"
Aug 2002
Team Italia

3·5·317 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
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.
ET_ is offline   Reply With Quote
Old 2017-07-14, 22:42   #11
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

3·1,087 Posts
Default

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
  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++
Attached Files
File Type: zip dMM.c.zip (783 Bytes, 34 views)

Last fiddled with by paulunderwood on 2017-07-14 at 22:49
paulunderwood is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
Prime95 Source Code required Qazi PrimeNet 1 2017-02-03 15:40
Triple check required houding PrimeNet 14 2015-12-21 09:34
New hardware - help required PageFault Software 2 2014-04-14 15:23
New Algorithm - Tester(s) required bearnol Factoring 29 2006-08-16 04:57
Which version of SHLWAPI.dll is required? gowen72 Software 1 2002-11-01 00:36

All times are UTC. The time now is 21:40.

Thu Jul 2 21:40:03 UTC 2020 up 99 days, 19:13, 0 users, load averages: 1.63, 2.30, 2.42

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2020, 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.