mersenneforum.org  

Go Back   mersenneforum.org > Extra Stuff > Programming

Reply
 
Thread Tools
Old 2015-08-29, 22:25   #1
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

65638 Posts
Default GWNUM program

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:

Code:
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.

Last fiddled with by paulunderwood on 2015-08-29 at 22:26
paulunderwood is offline   Reply With Quote
Old 2015-08-31, 08:39   #2
henryzz
Just call me Henry
 
henryzz's Avatar
 
"David"
Sep 2007
Cambridge (GMT/BST)

22×1,433 Posts
Default

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.
henryzz is online now   Reply With Quote
Old 2015-08-31, 12:44   #3
danaj
 
"Dana Jacobsen"
Feb 2011
Bangkok, TH

16128 Posts
Default

Quote:
Originally Posted by henryzz View Post
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.
danaj is offline   Reply With Quote
Old 2015-09-07, 05:29   #4
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

11×313 Posts
Cool

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:

Code:
gcc s2.c gwnum/gwnum.a gwnum/gwnum.ld -lm  -lpthread -lstdc++
s2.c:
Code:
#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?

Last fiddled with by paulunderwood on 2015-09-07 at 05:51
paulunderwood is offline   Reply With Quote
Old 2015-09-07, 13:51   #5
Prime95
P90 years forever!
 
Prime95's Avatar
 
Aug 2002
Yeehaw, FL

2·3·29·41 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
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.
Prime95 is online now   Reply With Quote
Old 2015-09-07, 14:04   #6
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

1101011100112 Posts
Default

Quote:
Originally Posted by Prime95 View Post
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.

Quote:
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.

Quote:
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.

Quote:
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.

Last fiddled with by paulunderwood on 2015-09-07 at 14:10
paulunderwood is offline   Reply With Quote
Old 2015-09-07, 19:28   #7
Prime95
P90 years forever!
 
Prime95's Avatar
 
Aug 2002
Yeehaw, FL

2×3×29×41 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
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.
Prime95 is online now   Reply With Quote
Old 2015-09-08, 06:05   #8
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

11×313 Posts
Default

Quote:
Originally Posted by Prime95 View Post
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:
Code:
#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

Last fiddled with by paulunderwood on 2015-09-08 at 06:25
paulunderwood is offline   Reply With Quote
Old 2015-09-08, 18:45   #9
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

11×313 Posts
Default

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().

Code:
#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 is offline   Reply With Quote
Old 2016-04-13, 18:53   #10
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

D7316 Posts
Default

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):

Code:
Likely prime!

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

Last fiddled with by paulunderwood on 2016-04-13 at 19:13
paulunderwood is offline   Reply With Quote
Old 2016-04-14, 16:33   #11
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

344310 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
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:

Code:
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?
paulunderwood 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 20:31.

Thu Oct 22 20:31:16 UTC 2020 up 42 days, 17:42, 0 users, load averages: 1.47, 1.93, 1.98

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.