mersenneforum.org  

Go Back   mersenneforum.org > Extra Stuff > Programming

Reply
 
Thread Tools
Old 2016-06-06, 20:50   #12
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

74158 Posts
Default

This is a straight Lucas test (with a line for threads) and it runs in excess of 375% on a 4 core system:

LucasPRP.c:

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

FILE *fd;
char string[5000000];
float fA;
int i, LEN, K, B, E, C;
unsigned init_va[1] = {2}, init_vb[1], A;
giant  n, gva, gvb, z, gD; 
gwnum  wva, wvb;
gwhandle *gwdata;


int Jacobi ( giant D, giant N, int LEN ){
  int j = 1;
  giant temp = newgiant ( LEN );
  giant n = newgiant ( LEN );
  gtog( N , n );
  while ( !isZero ( D ) ) {
    while ( !bitval ( D, 0 ) ) {
      gtogshiftright ( 1, D, D );      
      if(bitval(n,0)&&(( bitlen(n)==2&&bitval(n,1) ) || ( bitlen(n)>=3&&!bitval(n,2)&&bitval(n,1) ) || ( bitlen(n)>=3&&bitval(n,2)&&!bitval(n,1) )) )
        j = -j;
    }
    gtog ( D, temp );
    gtog ( n , D );
    gtog ( temp, n );
    if ( bitlen(D)>=2 && bitval(D,1) && bitval(D,0) && bitlen(n)>=2 && bitval(n,1) && bitval(n,0) )
      j = -j;
      modg ( n, D );
  }
  if ( isone ( n ) )
   return (j);
  else
    return(0);
}


int main ( int argc, char *argv[] ) {
  if( argc != 2 && argc != 6 ) {
    printf( "Usage: lucasPRP file_name (and optionally K B N C)\n" );
    return;
  }
  fd = fopen ( argv[1], "r" );
  fscanf ( fd, "%s", string );
  fclose ( fd );
  if ( argc == 6 ) {
    K = (double)atoi ( argv[2] );
    B = (unsigned long)atoi ( argv[3] );
    E = (unsigned long)atoi ( argv[4] );
    C = (signed long)atoi ( argv[5] );
  }
  LEN = strlen ( string ); LEN = ( LEN >> 2 ) + 8; 
  n = newgiant ( LEN );
  gva = newgiant ( LEN ); gvb = newgiant ( LEN ); z = newgiant ( LEN );
  gD = newgiant ( LEN );
  ctog ( string, n );
  A = 3;
  itog ( A*A-4, gD );
  while( Jacobi( gD, n , LEN ) != -1 ) {
    A++;
    if ( A > 1000000 ) {
      printf("Is this a square number?!\n");
      exit(0);
    }
    itog ( A*A-4, gD );
  }
  printf ( "Lucas testing on x^2 - %d*x + 1 ...\n", A );
  gwdata = (gwhandle*) malloc (sizeof ( gwhandle ) );
  gwinit ( gwdata );
  gwset_num_threads ( gwdata, 4 );
  if( argc == 2 )
    gwsetup_general_mod_giant( gwdata, n );
  else
    gwsetup ( gwdata, K, B, E, C);
  wva = gwalloc ( gwdata ); wvb = gwalloc ( gwdata );

  init_vb[0] = A;
  fA = (float) A;  

  LEN = bitlen ( n );
  binarytogw ( gwdata, init_va, 1, wva );
  binarytogw ( gwdata, init_vb, 1, wvb );
  for ( i = LEN; i >= 0; i-- ) {
    if ( bitval ( n, i ) ) {
      gwsafemul ( gwdata, wvb , wva );
      gwsmalladd ( gwdata, -fA, wva );
      gwsquare ( gwdata, wvb );
      gwsmalladd ( gwdata, -2.0, wvb );
    } else {
      gwsafemul ( gwdata, wva , wvb );
      gwsmalladd ( gwdata, -fA, wvb );
      gwsquare ( gwdata, wva );
      gwsmalladd ( gwdata, -2.0, wva );
    }
  }
  gwtogiant ( gwdata, wva, gva );
  gwtogiant ( gwdata, wvb, gvb );
  modg ( n, gva );
  itog ( A, z );
  subg ( z, gva );
  modg ( n, gvb );
  itog ( 2, z );
  subg ( z, gvb );
  if ( isZero ( gva ) && isZero ( gvb ) )
    printf ( "Is Lucas PRP!\n" );
  else
    printf ( "Not prime :(\n" );
}
//gwtogiant( gwdata, ws, z );
//gtoc( z, string, 100 );
//printf( "s%s\n", string );
//gcc -o lucasPRP lucasPRP.c gwnum/gwnum.a gwnum/gwnum.ld -lm -lpthread -lstdc++
paulunderwood is offline   Reply With Quote
Old 2016-06-06, 21:20   #13
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

74158 Posts
Default

Quote:
/* Prior to calling one of the gwsetup routines, you can tell the library */
/* how many threads it can use to perform a multiply. */

#define gwset_num_threads(h,n) ((h)->num_threads = n)
George, are you planning to thread other GWNUM functions?
paulunderwood is offline   Reply With Quote
Old 2016-06-07, 07:39   #14
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

F0D16 Posts
Default

This version of lucasPRP.c forward-FFTs two variables before multiplying them and squaring one. It is significantly quicker.

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

FILE *fd;
char string[5000000];
float fA;
int i, LEN, K, B, E, C;
unsigned init_va[1] = {2}, init_vb[1], A;
giant  n, gva, gvb, z, gD; 
gwnum  wva, wvb;
gwhandle *gwdata;


int Jacobi ( giant D, giant N, int LEN ){
  int j = 1;
  giant temp = newgiant ( LEN );
  giant n = newgiant ( LEN );
  gtog( N , n );
  while ( !isZero ( D ) ) {
    while ( !bitval ( D, 0 ) ) {
      gtogshiftright ( 1, D, D );      
      if(bitval(n,0)&&(( bitlen(n)==2&&bitval(n,1) ) || ( bitlen(n)>=3&&!bitval(n,2)&&bitval(n,1) ) || ( bitlen(n)>=3&&bitval(n,2)&&!bitval(n,1) )) )
        j = -j;
    }
    gtog ( D, temp );
    gtog ( n , D );
    gtog ( temp, n );
    if ( bitlen(D)>=2 && bitval(D,1) && bitval(D,0) && bitlen(n)>=2 && bitval(n,1) && bitval(n,0) )
      j = -j;
      modg ( n, D );
  }
  if ( isone ( n ) )
   return (j);
  else
    return(0);
}


int main ( int argc, char *argv[] ) {
  if( argc != 2 && argc != 6 ) {
    printf( "Usage: lucasPRP file_name (and optionally K B N C)\n" );
    return;
  }
  fd = fopen ( argv[1], "r" );
  fscanf ( fd, "%s", string );
  fclose ( fd );
  if ( argc == 6 ) {
    K = (double)atoi ( argv[2] );
    B = (unsigned long)atoi ( argv[3] );
    E = (unsigned long)atoi ( argv[4] );
    C = (signed long)atoi ( argv[5] );
  }
  LEN = strlen ( string ); LEN = ( LEN >> 2 ) + 8; 
  n = newgiant ( LEN );
  gva = newgiant ( LEN ); gvb = newgiant ( LEN ); z = newgiant ( LEN );
  gD = newgiant ( LEN );
  ctog ( string, n );
  A = 3;
  itog ( A*A-4, gD );
  while( Jacobi( gD, n , LEN ) != -1 ) {
    A++;
    if ( A > 1000000 ) {
      printf("Is this a square number?!\n");
      exit(0);
    }
    itog ( A*A-4, gD );
  }
  printf ( "Lucas testing on x^2 - %d*x + 1 ...\n", A );
  gwdata = (gwhandle*) malloc (sizeof ( gwhandle ) );
  gwinit ( gwdata );
  gwset_num_threads ( gwdata, 4 );
  if( argc == 2 )
    gwsetup_general_mod_giant( gwdata, n );
  else
    gwsetup ( gwdata, K, B, E, C);
  wva = gwalloc ( gwdata ); wvb = gwalloc ( gwdata );

  init_vb[0] = A;
  fA = (float) A;  

  LEN = bitlen ( n );
  binarytogw ( gwdata, init_va, 1, wva );
  binarytogw ( gwdata, init_vb, 1, wvb );
  for ( i = LEN; i >= 0; i-- ) {
    gwfft ( gwdata, wva, wva );
    gwfft ( gwdata, wvb, wvb );
    if ( bitval ( n, i ) ) {
      gwfftfftmul ( gwdata, wvb , wva, wva );
      gwsmalladd ( gwdata, -fA, wva );
      gwfftfftmul ( gwdata, wvb , wvb, wvb );
      gwsmalladd ( gwdata, -2.0, wvb );
    } else {
      gwfftfftmul ( gwdata, wvb , wva, wvb );
      gwsmalladd ( gwdata, -fA, wvb );
      gwfftfftmul ( gwdata, wva , wva, wva );
      gwsmalladd ( gwdata, -2.0, wva );
    }
  }
  gwtogiant ( gwdata, wva, gva );
  gwtogiant ( gwdata, wvb, gvb );
  modg ( n, gva );
  itog ( A, z );
  subg ( z, gva );
  modg ( n, gvb );
  itog ( 2, z );
  subg ( z, gvb );
  if ( isZero ( gva ) && isZero ( gvb ) )
    printf ( "Is Lucas PRP!\n" );
  else
    printf ( "Not prime :(\n" );
}
//gwtogiant( gwdata, ws, z );
//gtoc( z, string, 100 );
//printf( "s%s\n", string );
//gcc -o lucasPRP lucasPRP.c gwnum/gwnum.a gwnum/gwnum.ld -lm -lpthread -lstdc++
paulunderwood is offline   Reply With Quote
Old 2017-06-11, 19:01   #15
Batalov
 
Batalov's Avatar
 
"Serge"
Mar 2008
Phi(4,2^7658614+1)/2

2×7×683 Posts
Default

Added a few lines for you, so that you could read input from stdin (use filename = "-")
Code:
diff lucasPRP.c.orig lucasPRP.c
9c9
< char string[5000000];
---
> char s[5000000], *string=s;
48,50c48,55
<   fd = fopen ( argv[1], "r" );
<   fscanf ( fd, "%s", string );
<   fclose ( fd );
---
>   if(argv[1][0] == '-') {
>     gets(s);
>   } else {
>     fd = fopen ( argv[1], "r" );
>     fscanf ( fd, "%s", string );
>     fclose ( fd );
>   }
>   if((string=strchr(s,':')) && string[1]==' ') string+=2; else string=s;
Then you can run it like this, for example:
pfgw64 -k -od -q"L(81671)" | ./lucasPRP -
Batalov is offline   Reply With Quote
Old 2017-06-11, 19:41   #16
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

3,853 Posts
Default

Thanks for that Serge. I have made the change to my local copy.

gcc 6.3.0 complains with "lucasPRP.c:(.text+0x261): warning: the `gets' function is dangerous and should not be used." But what the hell!
paulunderwood is offline   Reply With Quote
Old 2017-06-11, 19:48   #17
Batalov
 
Batalov's Avatar
 
"Serge"
Mar 2008
Phi(4,2^7658614+1)/2

2·7·683 Posts
Default

Yes but it's true for fscanf just as well.
Here's how you can do with a bit of cleaning up:
Code:
# diff lucasPRP.c lucasPRP.c.orig
7d6
< #define MAXLEN 5000000
10c9
< char *st, *string;
---
> char string[5000000];
49,57c48,50
<   st = malloc(MAXLEN); if(!st) { printf( "out of memory\n"); exit(-1); }
<   if(!strcmp(argv[1],"-")) {
<     if(!fgets ( st, MAXLEN, stdin )) { printf( "didn't get any input\n"); exit(-1); }
<   } else {
<     fd = fopen ( argv[1], "r" );
<     if(!fgets ( st, MAXLEN, fd )) { printf( "didn't get any input\n"); exit(-1); }
<     fclose ( fd );
<   }
<   if((string=strchr(st,':')) && string[1]==' ') string+=2; else string=st;
---
>   fd = fopen ( argv[1], "r" );
>   fscanf ( fd, "%s", string );
>   fclose ( fd );
Batalov is offline   Reply With Quote
Old 2017-06-11, 20:07   #18
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

385310 Posts
Default

The code does not like square numbers -- needs some more work there.

And is says 2 and 3 are not prime, but the code is for really big numbers

Last fiddled with by paulunderwood on 2017-06-11 at 20:13
paulunderwood is offline   Reply With Quote
Old 2017-06-12, 06:51   #19
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

3,853 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
The code does not like square numbers -- needs some more work there.
Changing to the following line for checking squareness stops an int overflow:

Code:
    if ( A > 32768 ) {
paulunderwood is offline   Reply With Quote
Old 2017-06-18, 18:55   #20
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

74158 Posts
Default s2.c update

This version of s2.c has auto jacobi symbol selection, applies Serge's patch for piping input, and does the crude squareness test. I got rid of some gwcopy and used gwfft and gwswap. I am dismayed by its performance, being almost twice as slow as the above lucasPRP.c program

Code:
time ../../coding/gwnum/lucasPRP M611999-cofactor 1 2 611999 -1
Lucas testing on x^2 - 9*x + 1 ...
Is Lucas PRP!

real	2m21.692s
user	5m2.932s
sys	2m41.120s
Code:
time ../../coding/gwnum/s2 M611999-cofactor 1 2 611999 -1
Testing (x + 2)^n == 2*9 + 5 (mod n, x^2 - 9*x + 1)...
Likely prime!

real	4m14.164s
user	7m18.964s
sys	2m35.140s
s2.c:
Code:
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "gwnum/giants.h"
#include "gwnum/gwnum.h"
#include "gwnum/gwcommon.h"

#define MAXLEN 5000000

FILE *fd;
char *st, *string;
float a , b;
int i, LEN, K, B, E, C, A;
unsigned init_s[1] = {1}, init_t[1] = {2};
giant  n,  p,   s,  t,      z,     gD ; // p = n+1
gwnum          ws, wt, wv, wx, wy, wz; // wx, wz temps
gwhandle *gwdata;

int Jacobi ( giant D, giant N, int LEN ){
  int j = 1;
  giant temp = newgiant ( LEN );
  giant n = newgiant ( LEN );
  gtog( N , n );
  while ( !isZero ( D ) ) {
    while ( !bitval ( D, 0 ) ) {
      gtogshiftright ( 1, D, D );
      if(bitval(n,0)&&(( bitlen(n)==2&&bitval(n,1) ) || ( bitlen(n)>=3&&!bitval(n,2)&&bitval(n,1) ) || ( bitlen(n)>=3&&bitval(n,2)&&!bitval(n,1) )) )
        j = -j;
    }
    gtog ( D, temp );
    gtog ( n , D );
    gtog ( temp, n );
    if ( bitlen(D)>=2 && bitval(D,1) && bitval(D,0) && bitlen(n)>=2 && bitval(n,1) && bitval(n,0) )
      j = -j;
      modg ( n, D );
  }
  if ( isone ( n ) )
   return (j);
  else
    return(0);
}

int main ( int argc, char *argv[] ) {
  if( argc != 2 && argc != 6 ) {
    printf( "Usage: s2 file_name (and optionally K B N C)\n" );
    printf( "file_name can be replace with - if the number is piped in...\n" );
    printf( "E.g:  pfgw64 -k -f0 -od -q\"L(81671)\" | ./s2 -\n" );
    exit ( -1 );
  }

  st = malloc(MAXLEN); if(!st) { printf( "out of memory\n"); exit(-1); }
  if(!strcmp(argv[1],"-")) {
    if(!fgets ( st, MAXLEN, stdin )) { printf( "didn't get any input\n"); exit(-1); }
  } else {
    fd = fopen ( argv[1], "r" );
    if(!fgets ( st, MAXLEN, fd )) { printf( "didn't get any input\n"); exit(-1); }
     fclose ( fd );
  }
  if((string=strchr(st,':')) && string[1]==' ') string+=2; else string=st;

  if( argc == 6 ) {
    K = (double)atoi ( argv[2] );
    B = (unsigned long)atoi ( argv[3] );
    E = (unsigned long)atoi ( argv[4] );
    C = (signed long)atoi ( argv[5] );
  }

  LEN = strlen ( string ); LEN = ( LEN >> 2 ) + 8;
  n = newgiant ( LEN ); p = newgiant ( LEN ); gD = newgiant ( LEN );
  s = newgiant ( LEN ); t = newgiant ( LEN ); z = newgiant ( LEN );
  gwdata = (gwhandle*) malloc (sizeof ( gwhandle ) );
  gwinit ( gwdata );
  gwset_num_threads ( gwdata, 4 );
  ctog ( string, n );
  // Calculate the right Jacobi Sysbol
  if ( bitlen ( n ) >= 2 && bitval( n, 1 ) && bitval( n, 0 ) ) {
    A = 0;
  } else {
    itog ( 6, z );
    ctog ( string, p );
    modg ( z, p );
    ulsubg ( 5, p );
    if isZero ( p ) {
      A = 1;
    } else {
      A = 3;
      itog ( A*A-4, gD );
      while( Jacobi( gD, n , LEN ) != -1 ) {
        A++;
        if ( A > 32768 ) {
          printf("Is this a square number?!\n");
          exit(0);
        }
        itog ( A*A-4, gD );
      }
    }
  }
  printf ( "Testing (x + 2)^n == 2*%d + 5 (mod n, x^2 - %d*x + 1)...\n", A, A );
  a = (float) A;
  b = a + 2.0;
  //  gwsetup_general_mod_giant( gwdata, n );
  if( argc == 2 )
    gwsetup_general_mod_giant( gwdata, n );
  else
    gwsetup ( gwdata, K, B, E, C);
  ws = gwalloc ( gwdata ); wt = gwalloc ( gwdata ); wv = gwalloc ( gwdata );
  wx = gwalloc ( gwdata ); wy = gwalloc ( gwdata ); wz = gwalloc ( gwdata );

  itog ( 1, p ); addg ( n, p );
  LEN = bitlen ( p );
  binarytogw ( gwdata, init_s, 1, ws );
  binarytogw ( gwdata, init_t, 1, wt );
  for ( i = LEN-2; i >= 0; i-- ) {
    gwfft ( gwdata, ws, wx );
    gwfft ( gwdata, wt, wy );
    gwfftaddsub4 ( gwdata, wy, wx, wv, wz );
    gwfftfftmul ( gwdata, wv, wz, wz );
    gwsmallmul ( gwdata, a, ws );
    gwsmallmul ( gwdata, 2, wt );
    gwaddquick ( gwdata, wt, ws );
    gwfftmul ( gwdata, wx, ws );
    gwswap ( wt, wz );
    if ( bitval ( p, i ) ) {
      gwadd3 ( gwdata, wt, wt, wx );
      gwsub3 ( gwdata, wx, ws, wx);
      gwsmallmul ( gwdata, b, ws );
      gwadd ( gwdata, wt, ws );
      gwswap ( wx, wt );
    }
  }
  gwtogiant ( gwdata, ws, s );
  gwtogiant ( gwdata, wt, t );
  modg ( n, s );
  itog ( 2*(int)a+5, z );
  subg ( t, z );
  modg ( n, z );
  if ( isZero ( s ) && isZero ( z ) )
    printf ( "Likely prime!\n" );
  else
    printf ( "Not prime :(\n" );
}
//gwtogiant( gwdata, ws, z );
//gtoc( z, string, 100 );
//printf( "s%s\n", string );
//gcc -o s2 s2.c gwnum/gwnum.a gwnum/gwnum.ld -lm -lpthread -lstdc++
//http://www.mersenneforum.org/showthread.php?p=388342#post388342

Last fiddled with by paulunderwood on 2017-06-18 at 19:46
paulunderwood is offline   Reply With Quote
Old 2017-06-19, 13:53   #21
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

3,853 Posts
Default

(Note: the previous s2.c should say "(x + 2)^(n + 1)")

I ran some timings on L(2316773)...

Code:
time ./pfgw64-threaded -b2 -q"L(2316773)"
PFGW Version 3.8.3.64BIT.20161203.x86_Dev [GWNUM 28.6]

L(2316773) is 2-PRP! (1452.5474s+0.0192s)                                    

real	24m12.643s
user	71m28.216s
sys	1m50.824s
Code:
time ./pfgw64-threaded -tp -q"L(2316773)"
PFGW Version 3.8.3.64BIT.20161203.x86_Dev [GWNUM 28.6]

Primality testing L(2316773) [N+1, Brillhart-Lehmer-Selfridge]                                    
Running N+1 test using discriminant 2, base 1+sqrt(2)                                    
Calling Brillhart-Lehmer-Selfridge with factored part 0.00%
L(2316773) is Lucas PRP! (6905.8126s+0.0106s)                                    

real	115m5.852s
user	296m35.392s
sys	7m42.216s
Code:
time ../../coding/gwnum/lucasPRP Lucas_2316773
Lucas testing on x^2 - 3*x + 1 ...
Is Lucas PRP!

real	51m50.458s
user	148m41.324s
sys	5m29.136s
Code:
time ../../coding/gwnum/s2 Lucas_2316773 
Testing (x + 2)^(n + 1) == 2*1 + 5 (mod n, x^2 - 1*x + 1)...
Not prime :(

real	124m2.242s
user	319m45.840s
sys	8m17.864s
So either my s2.c code is in error or the heat of a very hot summer's day got to my heavily overclocked CPU. I suspect it is because I am not using gwsafemul as George suggested.

I am now trying gwadd instead of gwaddquick...

Last fiddled with by paulunderwood on 2017-06-19 at 14:27
paulunderwood is offline   Reply With Quote
Old 2017-06-19, 14:42   #22
Batalov
 
Batalov's Avatar
 
"Serge"
Mar 2008
Phi(4,2^7658614+1)/2

255A16 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
... I suspect it is because I am not using gwsafemul as George suggested.

I am now trying gwadd instead of gwaddquick...
While reviewing your code, I also had the same suspicion. It is "short and sweet" (without gwsafemul for the first and last 50 iterations) - but here's the price.
Batalov is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
LLR V3.8.2 using gwnum 26.2 is available! Jean Penné Software 25 2010-11-01 15:18
GWNUM? Unregistered Information & Answers 3 2010-09-12 19:52
GWNUM library and llr leizhoucn Programming 2 2007-11-05 09:34
Compiling GMP-ECM With GWNUM tmorrow GMP-ECM 5 2007-04-04 00:39
GWNUM as DLL? Cyclamen Persicum Software 1 2007-01-02 20:53

All times are UTC. The time now is 07:13.


Sun Oct 17 07:13:26 UTC 2021 up 86 days, 1:42, 0 users, load averages: 1.12, 1.47, 1.39

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