GWNUM program

[paulunderwood]

This is another call for a simple 50ish line program to be written with GWNUM to implement Tony Reix's Digraph for 2^n-3:

for(n=4, 10000, N=2^n-3; S0=8; x=S0; for(i=1, n-1, x=Mod(x^2-2,N)); if(x==S0, print1(n,", ")))

To reiterate, the code should look something like:

• include GWNUM
• read some parameter(s) from a file or a decimal string
• initialize a few things
• loop over bits including mults, squarings, additions etc. including (special) mod reductions
• test a condition
• report result to screen and file

The reason is to write a small program to illustrate the use of GWNUM. I know we could whittle down the LLR source code but I though some one could just write a program.

[henryzz]

There are two possibilities. A pfgw script could do what you want. Another option is talking to Jean about adding it to llr.
Both would be much easier than a new program.

[danaj]

Originally Posted by henryzz

There are two possibilities. A pfgw script could do what you want. Another option is talking to Jean about adding it to llr.
Both would be much easier than a new program.

From my position, the specific program Paul asked for was an *example*. The point isn't to solve the Reix test, but to illustrate writing a very simple gwnum program. It's "Hello World" so someone familiar with GMP could get a good idea what to do.

I, for instance, am more interested in doing a BPSW test, then adding a nextprime / prevprime framework. I'm not interested in calling the PFGW executable for each value (Niceky's cglp4 does this), and I don't think a PFGW script will do this properly.

I admit there is a "Dear Lazyweb" prefix to this request.

[paulunderwood]

I have hacked it! To get it running I needed to "make" llr64 in order to get gwnum.a in place. The following script is put in the root directory of the LLR source tree and compiled with:

gcc s2.c gwnum/gwnum.a gwnum/gwnum.ld -lm  -lpthread -lstdc++

s2.c:

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

FILE *fd;
char string[1000000];
int i, a = 0; // We want jacobi(a^2-4,n)==-1
unsigned long int LEN;
unsigned long int para_a[1] = {0}, para_b[1] = {2}; // n == 3 (mod 4), b = a+2
unsigned long int init_s[1] = {1}, init_t[1] = {2};
giant  n,  p,  s,  t,      z; // p = n+1
gwnum wn, wp, ws, wt, wx, wz; // wx, wz temps
gwhandle *gwdata;

int main ( int argc, char *argv[] ) {
  fd = fopen ( argv[1], "r" );
  fscanf ( fd, "%s", string );
  fclose ( fd );
  
  LEN = strlen ( string ); LEN = ( LEN >> 2 ) + 8; 
  n = newgiant ( LEN ); p = newgiant ( LEN ); s = newgiant ( LEN );
  t = newgiant ( LEN ); z = newgiant ( LEN );

  gwdata = (gwhandle*) malloc ( sizeof ( gwhandle ) );
  gwinit (gwdata);
  ctog ( string, n );
  gwsetup_general_mod_giant( gwdata, n );
  wn = gwalloc ( gwdata ); wp = gwalloc ( gwdata ); ws = gwalloc ( gwdata );
  wt = gwalloc ( gwdata ); wx = gwalloc ( gwdata ); wz = gwalloc ( gwdata );

  itog ( 1, p ); addg ( n, p );
  LEN = bitlen ( p );
  binary64togw ( gwdata, init_s, 1, ws );
  binary64togw ( gwdata, init_t, 1, wt );
  for ( i = LEN-2; i >= 0; i-- ) {
    binary64togw ( gwdata, para_a, 1, wz );
    gwsafemul ( gwdata, ws, wz );
    gwaddquick ( gwdata, wt, wz );
    gwaddquick ( gwdata, wt, wz );
    gwsafemul ( gwdata, ws, wz );
    gwsub3quick ( gwdata, wt, ws, wx );
    gwaddquick ( gwdata, ws, wt );
    gwsafemul ( gwdata, wx, wt );
    gwcopy ( gwdata, wz, ws );
    if ( bitval ( p, i ) ) {
      binary64togw ( gwdata, para_b, 1, wz );
      gwsafemul ( gwdata, ws, wz );
      gwaddquick ( gwdata, wt, wz );
      gwaddquick ( gwdata, wt, wt );
      gwsubquick ( gwdata, ws, wt );
      gwcopy ( gwdata, wz, ws );
    }       
  }
  itog ( 2*a+5, z );
  gwtogiant ( gwdata, wt, t );
  subg ( z,  t );
  if ( gwiszero ( gwdata, ws ) &&  isZero ( t ) )
    printf ( "Likely prime!\n" );
  else
    printf ( "Not prime :(\n" );
}
//gwtogiant( gwdata, ws, s );
//gtoc( s, string, 100 );
//printf( "%s\n", string );
//gcc s2.c gwnum/gwnum.a gwnum/gwnum.ld -lm  -lpthread -lstdc++
//http://www.mersenneforum.org/showpost.php?p=388342&postcount=1

Any ideas on how to speed up what is already a fast script -- thanks to GW and JP and the late RC -- will be welcome.

At the moment the program requires a number (==3 (mod 4)) to be in a file. Will I need to write my own jacobi() function?

[Prime95]

Originally Posted by paulunderwood

Any ideas on how to speed up what is already a fast script

gwset_general_mod is the slowest way to use gwnum. gwnum does better testing numbers of the form k*b^n+c. My understanding is you are testing numbers where k=1, b=2, c=-3.

The last in a series of gwadds or gwsubs should not be of the "quick" variety.

Pull the binarytogw out of the loop, only do it once. Consider pre-FFTing the value (with gwfft) then use gwfftmul instead of gwsafemul.

The last gwcopy to wz from ws seems pointless.

[paulunderwood]

Originally Posted by Prime95

gwset_general_mod is the slowest way to use gwnum. gwnum does better testing numbers of the form k*b^n+c. My understanding is you are testing numbers where k=1, b=2, c=-3.

Sorry for the confusion; This program implements my 2 selfridge combined Fermat PRP and Lucas PRP test, not Tony's idea. It is meant to be illustrative of what can be done with gwnum.

The last in a series of gwadds or gwsubs should not be of the "quick" variety.

As far as I can see, since the final bit of p=n+1 is zero, the if is not entered into and so the mod is handled by the multiplication functions.

Pull the binarytogw out of the loop, only do it once. Consider pre-FFTing the value (with gwfft) then use gwfftmul instead of gwsafemul.

Sound advice.

The last gwcopy to wz from ws seems pointless.

I think it is needed because I am altering ws and wt in each loop. However, with some clever trickery I might be able to eliminate it.

[Prime95]

Originally Posted by paulunderwood

As far as I can see, since the final bit of p=n+1 is zero, the if is not entered into and so the mod is handled by the multiplication functions.

If you do an addquick or subquick before a gwmul (or gwsafemul, gwsquare, etc.) then you may get a roundoff error in the multiplication.

[paulunderwood]

Originally Posted by Prime95

If you do an addquick or subquick before a gwmul (or gwsafemul, gwsquare, etc.) then you may get a roundoff error in the multiplication.

I now see what you wrote. Thanks for being patient. I have not yet used gwmul() etc.

I found gwiszero() to be buggy.

The new code for s2.c is:

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

FILE *fd;
char string[1000000];
float a = 0.0, b = 2.0; // We want minimum a: jacobi(a^2-4,n)==-1, b = a+2
int i, LEN;
unsigned init_s[1] = {1}, init_t[1] = {2};
giant  n,  p,   s,  t,      z; // p = n+1
gwnum          ws, wt, wx, wz; // wx, wz temps
gwhandle *gwdata;

int main ( int argc, char *argv[] ) {
  fd = fopen ( argv[1], "r" );
  fscanf ( fd, "%s", string );
  fclose ( fd );
  
  LEN = strlen ( string ); LEN = ( LEN >> 2 ) + 8; 
  n = newgiant ( LEN ); p = newgiant ( LEN );
  s = newgiant ( LEN ); t = newgiant ( LEN ); z = newgiant ( LEN );

  gwdata = (gwhandle*) malloc (sizeof ( gwhandle ) );
  gwinit ( gwdata );
  ctog ( string, n );
  gwsetup_general_mod_giant( gwdata, n );
  ws = gwalloc ( gwdata ); wt = gwalloc ( gwdata );
  wx = gwalloc ( gwdata ); wz = gwalloc ( gwdata );

  binarytogw ( gwdata, init_s, 1, ws );
  binarytogw ( gwdata, init_t, 1, wt );
  itog ( 1, p ); addg ( n, p );
  LEN = bitlen ( p );
  for ( i = LEN-2; i >= 0; i-- ) {
    gwcopy ( gwdata, ws, wz );
    gwsmallmul ( gwdata, a, wz );
    gwaddquick ( gwdata, wt, wz );
    gwadd ( gwdata, wt, wz );
    gwsafemul ( gwdata, ws, wz );
    gwsub3 ( gwdata, wt, ws, wx );
    gwadd ( gwdata, ws, wt );
    gwsafemul ( gwdata, wx, wt );
    if ( bitval ( p, i ) ) {
      gwcopy ( gwdata, wz, wx );
      gwsmallmul ( gwdata, b, wx );
      gwadd ( gwdata, wt, wx );
      gwaddquick ( gwdata, wt, wt );
      gwsub ( gwdata, wz, wt );
      gwcopy ( gwdata, wx, ws );
    } else gwcopy ( gwdata, wz, ws );      
  }
  gwtogiant ( gwdata, ws, s );
  gwtogiant ( gwdata, wt, t );
  itog ( 2*(int)a+5, z );
  subg ( t, 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 s2.c gwnum/gwnum.a gwnum/gwnum.ld -lm -lpthread -lstdc++
//http://www.mersenneforum.org/showpost.php?p=388342&postcount=1

[paulunderwood]

I have changed the interface slightly; You now must enter the value of "a" -- minimal a: jacobi(a^2-4,n)==-1 -- on the command line. The interface is not important to this exercise. This is my final version.

I have re-written the main loop, getting rid of some gwcopy().

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

FILE *fd;
char string[1000000];
float a , b;
int i, LEN;
unsigned init_s[1] = {1}, init_t[1] = {2};
giant  n,  p,   s,  t,      z; // p = n+1
gwnum          ws, wt, wx, wz; // wx, wz temps
gwhandle *gwdata;

int main ( int argc, char *argv[] ) {
  if( argc != 3 ) {
    printf( "Usage: s2 file_name a (where a is min: jacobi(a^2-4,n)==-1)\n" );
    return;
  }
  fd = fopen ( argv[1], "r" );
  fscanf ( fd, "%s", string );
  fclose ( fd );
  a = (float)atoi ( argv[2] ); b = a + 2.0;

  LEN = strlen ( string ); LEN = ( LEN >> 2 ) + 8; 
  n = newgiant ( LEN ); p = newgiant ( LEN );
  s = newgiant ( LEN ); t = newgiant ( LEN ); z = newgiant ( LEN );

  gwdata = (gwhandle*) malloc (sizeof ( gwhandle ) );
  gwinit ( gwdata );
  ctog ( string, n );
  gwsetup_general_mod_giant( gwdata, n );
  ws = gwalloc ( gwdata ); wt = gwalloc ( gwdata );
  wx = gwalloc ( gwdata ); wz = gwalloc ( gwdata );

  binarytogw ( gwdata, init_s, 1, ws );
  binarytogw ( gwdata, init_t, 1, wt );
  itog ( 1, p ); addg ( n, p );
  LEN = bitlen ( p );
  gwcopy ( gwdata, ws, wz );
  for ( i = LEN-2; i >= 0; i-- ) {
    gwcopy ( gwdata, wz, ws );
    gwsmallmul ( gwdata, a, wz );
    gwaddquick ( gwdata, wt, wz );
    gwadd ( gwdata, wt, wz );
    gwsafemul ( gwdata, ws, wz );
    gwsub3 ( gwdata, wt, ws, wx );
    gwadd ( gwdata, ws, wt );
    gwsafemul ( gwdata, wx, wt );
    if ( bitval ( p, i ) ) {
      gwcopy ( gwdata, wz, ws );
      gwsmallmul ( gwdata, b, wz );
      gwadd ( gwdata, wt, wz );
      gwaddquick ( gwdata, wt, wt );
      gwsubquick ( gwdata, ws, wt );
    }     
  }
  gwtogiant ( gwdata, wz, s );
  gwtogiant ( gwdata, wt, t );
  itog ( 2*(int)a+5, z );
  subg ( t, 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/showpost.php?p=388342&postcount=1

[paulunderwood]

I added gwset_num_threads ( gwdata, 4 ); to the above program and ran it against this 971486-bit number. The percentage usage was heavily dependent on how many other tasks I was running. After killing a few, like Iceweasel, I noticed that the CPU usage went up to 255% of an idle 4770k (no HT). Can I expect more CPU usage from big mega-primes/PRPs? Here is the result (with some processes being killed throughout the run):

Likely prime!

real	29m42.842s
user	72m11.292s
sys	2m36.888s

[paulunderwood]

Originally Posted by paulunderwood

I added gwset_num_threads ( gwdata, 4 ); to the above program and ran it against this 971486-bit number. The percentage usage was heavily dependent on how many other tasks I was running. After killing a few, like Iceweasel, I noticed that the CPU usage went up to 255% of an idle 4770k (no HT). Can I expect more CPU usage from big mega-primes/PRPs?

To answer my own question: I tested Anand Nair's mega-PRP (1042896 decimal digits) Mersenne co-factor ((2^3464473-1)/604874508299177) using gwsetup() for special mod for which I change the interface and added modg ( n, s ); after the main loop to handle the fraction. The computation used 224% of CPU and resulted in:

Likely prime!

real	123m28.272s
user	269m14.144s
sys	7m14.316s

As this test is doing two multiplications per loop, would it be possible/advisable to thread them? if so, how?