Replace isaac with xorshift

No difference in performance. rdrand instruction alone is too slow.

14 files changed, 79 insertions, 271 deletions

diff --git a/Makefile.am b/Makefile.am
index 42056e9..362424f 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -8,9 +8,9 @@ ACLOCAL_AMFLAGS = -I m4
 man1_MANS = flam3-animate.man flam3-genome.man flam3-render.man flam3-convert.man
 bin_PROGRAMS = flam3-animate flam3-genome flam3-render flam3-convert
 lib_LTLIBRARIES = libflam3.la
-include_HEADERS = flam3.h isaac.h isaacs.h rect.c
+include_HEADERS = flam3.h xorshift.h rect.c
 
-libflam3_la_SOURCES = flam3.c filters.c parser.c variations.c interpolation.c palettes.c jpeg.c png.c isaac.c
+libflam3_la_SOURCES = flam3.c filters.c parser.c variations.c interpolation.c palettes.c jpeg.c png.c xorshift.c
 libflam3_la_LDFLAGS = -no-undefined -ljpeg -lpng -lz -lpthread
 
 flam3_genome_SOURCES = flam3-genome.c docstring.c
@@ -28,7 +28,7 @@ flam3_convert_LDADD = libflam3.la -lm
 pkgdata_DATA = flam3-palettes.xml
 
 EXTRA_DIST = rect.c flam3.h palettes.h variations.h interpolation.h parser.h \
-   private.h isaac.h isaacs.h img.h test.flam3 vidres.flam3 \
+   private.h img.h test.flam3 vidres.flam3 \
    flam3.dsw render.dsp animate.dsp genome.dsp convert.dsp \
    mkinstalldirs README.txt COPYING.txt flam3-palettes.xml \
    flam3-animate.man flam3-render.man flam3-genome.man flam3-convert.man 
diff --git a/README.txt b/README.txt
index aaf0f63..e7a30cc 100644
--- a/README.txt
+++ b/README.txt
@@ -70,7 +70,6 @@ ss              1               size scale, multiply size (in pixels) of all fra
 jpeg            NA              jpeg quality for compression, default is native jpeg default
 format          png             "jpg" or "ppm" or "png"
 pixel_aspect    1.0             aspect ratio of pixels (width over height), eg 0.90909 for NTSC
-isaac_seed      random          string to be used in generating random seed.  defaults to time(0)
 seed            random          integer seed for random numbers, defaults to time+pid.  deprecated.
 nthreads        auto            number of threads to use (render and animate)
 verbose         0               if non-zero then print progress meter on stderr
diff --git a/docstring.c b/docstring.c
index 7e6357d..4e72ff9 100644
--- a/docstring.c
+++ b/docstring.c
@@ -68,7 +68,6 @@ static char *the_docstring1 =
 "format          png         jpg or ppm or png\n"
 "pixel_aspect    1.0         aspect ratio of pixels (width over height), eg 0.90909 for NTSC\n"
 "seed            random      integer seed for random numbers, defaults to time+pid\n"
-"isaac_seed      random      character-based seed for iteration loop randomness, defaults to time\n"
 "nthreads        0           number of threads to use for render.  default auto-detects.\n"
 "verbose         0           if non-zero then print progress meter on stderr\n"
 "bits            33          also 32 or 64: sets bit-width of internal buffers (33 means 32-bit floating-point)\n"
diff --git a/flam3-animate.c b/flam3-animate.c
index 3c32628..ccf0b5a 100644
--- a/flam3-animate.c
+++ b/flam3-animate.c
@@ -19,7 +19,6 @@
 
 #include "private.h"
 #include "img.h"
-#include "isaacs.h"
 
 int main(int argc, char **argv) {
   char *ai, *fname;
diff --git a/flam3-genome.c b/flam3-genome.c
index 4551274..9ab08ad 100644
--- a/flam3-genome.c
+++ b/flam3-genome.c
@@ -872,7 +872,7 @@ main(argc, argv)
                int mutmeth;
 
                parent0 = string_to_cp(mutate, &parent0_n);
-               flam3_copy(&selp0, &(parent0[((unsigned)irand(&f.rc))%parent0_n]));
+               flam3_copy(&selp0, &(parent0[(xorshift_step(&f.rc))%parent0_n]));
                flam3_copy(&cp_orig, &selp0);
                aselp0 = &selp0;
                aselp1 = NULL;
@@ -917,8 +917,8 @@ main(argc, argv)
                parent0 = string_to_cp(cross0, &parent0_n);
                parent1 = string_to_cp(cross1, &parent1_n);
 
-               i0 = ((unsigned)irand(&f.rc))%parent0_n;
-               i1 = ((unsigned)irand(&f.rc))%parent1_n;
+               i0 = (xorshift_step(&f.rc))%parent0_n;
+               i1 = (xorshift_step(&f.rc))%parent1_n;
 
                flam3_copy(&selp0, &(parent0[i0]));
                flam3_copy(&selp1, &(parent1[i1]));
diff --git a/flam3-render.c b/flam3-render.c
index ddab377..e3acd64 100644
--- a/flam3-render.c
+++ b/flam3-render.c
@@ -30,7 +30,6 @@
 
 #include "private.h"
 #include "img.h"
-#include "isaacs.h"
 
 
 int calc_nstrips(flam3_frame *spec) {
diff --git a/flam3.c b/flam3.c
index 6a02335..ebf6aa5 100644
--- a/flam3.c
+++ b/flam3.c
@@ -28,6 +28,7 @@
 #include "parser.h"
 #include "filters.h"
 #include "palettes.h"
+#include "xorshift.h"
 #include <limits.h>
 #include <locale.h>
 #include <math.h>
@@ -250,9 +251,9 @@ int flam3_iterate(flam3_genome *cp, int n, int fuse,  double *samples, unsigned
    
 //         fn = xform_distrib[ lastxf*CHOOSE_XFORM_GRAIN + (((unsigned)irand(rc)) % CHOOSE_XFORM_GRAIN)];
       if (cp->chaos_enable)
-         fn = xform_distrib[ lastxf*CHOOSE_XFORM_GRAIN + (((unsigned)irand(rc)) & CHOOSE_XFORM_GRAIN_M1)];
+         fn = xform_distrib[ lastxf*CHOOSE_XFORM_GRAIN + (xorshift_step(rc) & CHOOSE_XFORM_GRAIN_M1)];
       else
-         fn = xform_distrib[ ((unsigned)irand(rc)) & CHOOSE_XFORM_GRAIN_M1 ];
+         fn = xform_distrib[ xorshift_step(rc) & CHOOSE_XFORM_GRAIN_M1 ];
       
       if (apply_xform(cp, fn, p, q, rc)>0) {
          consec ++;
@@ -2474,36 +2475,16 @@ double flam3_random11() {
 
 /* This function must be called prior to rendering a frame */
 void flam3_init_frame(flam3_frame *f) {
-
-   char *ai;
-   char *isaac_seed = args("isaac_seed",NULL);
-   long int default_isaac_seed = (long int)time(0);  
-
-   /* Clear out the isaac state */
-   memset(f->rc.randrsl, 0, RANDSIZ*sizeof(ub4));
-
-   /* Set the isaac seed */
-   if (NULL == isaac_seed) {
-      int lp;
-      /* No isaac seed specified.  Use the system time to initialize. */
-      for (lp = 0; lp < RANDSIZ; lp++)
-         f->rc.randrsl[lp] = default_isaac_seed;
-   } else {
-      /* Use the specified string */
-      strncpy((char *)&f->rc.randrsl,(const char *)isaac_seed, RANDSIZ*sizeof(ub4));
-   }
-
-   /* Initialize the random number generator */
-   irandinit(&f->rc,1);
+	xorshift_seed (&f->rc);
 }
 
 /* returns uniform variable from ISAAC rng */
 double flam3_random_isaac_01(randctx *ct) {
-   return ((int)irand(ct) & 0xfffffff) / (double) 0xfffffff;
+   return ((uint32_t) xorshift_step (ct) & 0xfffffff) / (double) 0xfffffff;
 }
 
 double flam3_random_isaac_11(randctx *ct) {
-   return (((int)irand(ct) & 0xfffffff) - 0x7ffffff) / (double) 0x7ffffff;
+   return (((uint32_t) xorshift_step(ct) & 0xfffffff) - 0x7ffffff) / (double) 0x7ffffff;
 }
 
 int flam3_random_bit() {
@@ -2521,8 +2502,7 @@ int flam3_random_bit() {
 }
 
 int flam3_random_isaac_bit(randctx *ct) {
-   int tmp = irand(ct);
-   return tmp & 1;
+   return xorshift_step(ct) & 1;
 }
 
 static double round6(double x) {
@@ -2882,7 +2862,7 @@ void flam3_mutate(flam3_genome *cp, int mutate_mode, int *ivars, int ivars_n, in
       flam3_random(&mutation, ivars, ivars_n, sym, 2);
       
       /* Which xform do we mutate? */
-      modxf = ((unsigned)irand(rc)) % cp->num_xforms;
+      modxf = (xorshift_step(rc)) % cp->num_xforms;
       
       add_to_action(action,"mutate xform ");
       sprintf(ministr,"%d coefs",modxf);
@@ -2905,8 +2885,8 @@ void flam3_mutate(flam3_genome *cp, int mutate_mode, int *ivars, int ivars_n, in
       
    } else if (mutate_mode == MUTATE_POST_XFORMS) {
    
-      int b = 1 + ((unsigned)irand(rc))%6;
-      int same = ((unsigned)irand(rc))&3; /* 25% chance of using the same post for all of them */
+      int b = 1 + (xorshift_step(rc))%6;
+      int same = (xorshift_step(rc))&3; /* 25% chance of using the same post for all of them */
       
       sprintf(ministr,"(%d%s)",b,(same>0) ? " same" : "");
       add_to_action(action,"mutate post xforms ");
@@ -3014,7 +2994,7 @@ void flam3_mutate(flam3_genome *cp, int mutate_mode, int *ivars, int ivars_n, in
       }
    } else if (mutate_mode == MUTATE_DELETE_XFORM) {
    
-      int nx = ((unsigned)irand(rc))%cp->num_xforms;
+      int nx = (xorshift_step(rc))%cp->num_xforms;
       sprintf(ministr,"%d",nx);
       add_to_action(action,"mutate delete xform ");
       add_to_action(action,ministr);
@@ -3936,16 +3916,11 @@ double flam3_dimension(flam3_genome *cp, int ntries, int clip_to_camera) {
   double bmax[2];
   double d2max;
   int lp;
-  long int default_isaac_seed = (long int)time(0);
   randctx rc;
   int SBS = 10000;
   int i, n1=0, n2=0, got, nclipped;
 
-  /* Set up the isaac rng */
-  for (lp = 0; lp < RANDSIZ; lp++)
-     rc.randrsl[lp] = default_isaac_seed;
-
-  irandinit(&rc,1);
+  xorshift_seed (&rc);
 
   if (ntries < 2) ntries = 3000*1000;
 
@@ -4057,15 +4032,9 @@ double flam3_lyapunov(flam3_genome *cp, int ntries) {
   unsigned short *xform_distrib;
 
   int lp;
-  long int default_isaac_seed = (long int)time(0);
   randctx rc;
 
-  /* Set up the isaac rng */
-  for (lp = 0; lp < RANDSIZ; lp++)
-     rc.randrsl[lp] = default_isaac_seed;
-
-  irandinit(&rc,1);
-
+  xorshift_seed (&rc);
 
   if (ntries < 1) ntries = 10000;
 
diff --git a/flam3.h b/flam3.h
index f46e09e..a66e40e 100644
--- a/flam3.h
+++ b/flam3.h
@@ -22,7 +22,7 @@
 
 #include <stdio.h>
 #include <libxml/parser.h>
-#include "isaac.h"
+#include "xorshift.h"
 
 #if defined(_MSC_VER) /* VC++ */
 #include <windows.h>
diff --git a/isaac.c b/isaac.c
deleted file mode 100644
index 014c4d9..0000000
--- a/isaac.c
+++ /dev/null
@@ -1,110 +0,0 @@
-/*
-------------------------------------------------------------------------------
-isaac.c: By Bob Jenkins.  My random number generator, ISAAC.
-MODIFIED:
-  960327: Creation (addition of randinit, really)
-  970719: use context, not global variables, for internal state
-  980324: make a portable version
-  991209: modified for inclusion with GNU Backgammon by Gary Wong
-  070121: modified for inclusion with flam3 by Erik Reckase
-------------------------------------------------------------------------------
-*/
-
-#include "isaacs.h"
-#include "isaac.h"
-
-#define ind(mm,x)  ((mm)[(x>>2)&(RANDSIZ-1)])
-#define rngstep(mix,a,b,mm,m,m2,r,x) \
-{ \
-  x = *m;  \
-  a = ((a^(mix)) + *(m2++)) & 0xffffffff; \
-  *(m++) = y = (ind(mm,x) + a + b) & 0xffffffff; \
-  *(r++) = b = (ind(mm,y>>RANDSIZL) + x) & 0xffffffff; \
-}
-
-void     isaac(randctx *ctx)
-{
-   register ub4 a,b,x,y,*m,*mm,*m2,*r,*mend;
-   mm=ctx->randmem; r=ctx->randrsl;
-   a = ctx->randa; b = (ctx->randb + (++ctx->randc)) & 0xffffffff;
-   for (m = mm, mend = m2 = m+(RANDSIZ/2); m<mend; )
-   {
-      rngstep( a<<13, a, b, mm, m, m2, r, x);
-      rngstep( a>>6 , a, b, mm, m, m2, r, x);
-      rngstep( a<<2 , a, b, mm, m, m2, r, x);
-      rngstep( a>>16, a, b, mm, m, m2, r, x);
-   }
-   for (m2 = mm; m2<mend; )
-   {
-      rngstep( a<<13, a, b, mm, m, m2, r, x);
-      rngstep( a>>6 , a, b, mm, m, m2, r, x);
-      rngstep( a<<2 , a, b, mm, m, m2, r, x);
-      rngstep( a>>16, a, b, mm, m, m2, r, x);
-   }
-   ctx->randb = b; ctx->randa = a;
-}
-
-
-#define mix(a,b,c,d,e,f,g,h) \
-{ \
-   a^=b<<11; d+=a; b+=c; \
-   b^=c>>2;  e+=b; c+=d; \
-   c^=d<<8;  f+=c; d+=e; \
-   d^=e>>16; g+=d; e+=f; \
-   e^=f<<10; h+=e; f+=g; \
-   f^=g>>4;  a+=f; g+=h; \
-   g^=h<<8;  b+=g; h+=a; \
-   h^=a>>9;  c+=h; a+=b; \
-}
-
-/* if (flag==TRUE), then use the contents of randrsl[] to initialize mm[]. */
-void irandinit(randctx *ctx, word     flag)
-{
-   word i;
-   ub4 a,b,c,d,e,f,g,h;
-   ub4 *m,*r;
-   ctx->randa = ctx->randb = ctx->randc = 0;
-   m=ctx->randmem;
-   r=ctx->randrsl;
-   a=b=c=d=e=f=g=h=0x9e3779b9;  /* the golden ratio */
-
-   for (i=0; i<4; ++i)          /* scramble it */
-   {
-     mix(a,b,c,d,e,f,g,h);
-   }
-
-   if (flag)
-   {
-     /* initialize using the contents of r[] as the seed */
-     for (i=0; i<RANDSIZ; i+=8)
-     {
-       a+=r[i  ]; b+=r[i+1]; c+=r[i+2]; d+=r[i+3];
-       e+=r[i+4]; f+=r[i+5]; g+=r[i+6]; h+=r[i+7];
-       mix(a,b,c,d,e,f,g,h);
-       m[i  ]=a; m[i+1]=b; m[i+2]=c; m[i+3]=d;
-       m[i+4]=e; m[i+5]=f; m[i+6]=g; m[i+7]=h;
-     }
-     /* do a second pass to make all of the seed affect all of m */
-     for (i=0; i<RANDSIZ; i+=8)
-     {
-       a+=m[i  ]; b+=m[i+1]; c+=m[i+2]; d+=m[i+3];
-       e+=m[i+4]; f+=m[i+5]; g+=m[i+6]; h+=m[i+7];
-       mix(a,b,c,d,e,f,g,h);
-       m[i  ]=a; m[i+1]=b; m[i+2]=c; m[i+3]=d;
-       m[i+4]=e; m[i+5]=f; m[i+6]=g; m[i+7]=h;
-     }
-   }
-   else
-   {
-     for (i=0; i<RANDSIZ; i+=8)
-     {
-       /* fill in mm[] with messy stuff */
-       mix(a,b,c,d,e,f,g,h);
-       m[i  ]=a; m[i+1]=b; m[i+2]=c; m[i+3]=d;
-       m[i+4]=e; m[i+5]=f; m[i+6]=g; m[i+7]=h;
-     }
-   }
-
-   isaac(ctx);            /* fill in the first set of results */
-   ctx->randcnt=RANDSIZ;  /* prepare to use the first set of results */
-}
diff --git a/isaac.h b/isaac.h
deleted file mode 100644
index 6554569..0000000
--- a/isaac.h
+++ /dev/null
@@ -1,56 +0,0 @@
-/*
-------------------------------------------------------------------------------
-isaac.h: definitions for a random number generator
-MODIFIED:
-  960327: Creation (addition of randinit, really)
-  970719: use context, not global variables, for internal state
-  980324: renamed seed to flag
-  980605: recommend RANDSIZL=4 for noncryptography.
-  991209: modified for inclusion with GNU Backgammon by Gary Wong
-  070121: modified for inclusion with flam3 by Erik Reckase
-------------------------------------------------------------------------------
-*/
-
-#include "isaacs.h"
-
-#ifndef _ISAAC_H_
-#define _ISAAC_H_
-
-#define RANDSIZL   (4)  /* I recommend 8 for crypto, 4 for simulations */
-#define RANDSIZ    (1<<RANDSIZL)
-
-/* context of random number generator */
-struct randctx
-{
-  ub4 randcnt;
-  ub4 randrsl[RANDSIZ];
-  ub4 randmem[RANDSIZ];
-  ub4 randa;
-  ub4 randb;
-  ub4 randc;
-};
-typedef  struct randctx  randctx;
-
-/*
-------------------------------------------------------------------------------
- If (flag==TRUE), then use the contents of randrsl[0..RANDSIZ-1] as the seed.
-------------------------------------------------------------------------------
-*/
-void irandinit( randctx *r, word flag );
-
-void isaac( randctx *r );
-
-
-/*
-------------------------------------------------------------------------------
- Call irand(/o_ randctx *r _o/) to retrieve a single 32-bit random value
-------------------------------------------------------------------------------
-*/
-#define irand(r) \
-   (!(r)->randcnt-- ? \
-     (isaac(r), (r)->randcnt=RANDSIZ-1, (r)->randrsl[(r)->randcnt]) : \
-     (r)->randrsl[(r)->randcnt])
-
-#endif
-
-
diff --git a/isaacs.h b/isaacs.h
deleted file mode 100644
index 687d8ae..0000000
--- a/isaacs.h
+++ /dev/null
@@ -1,46 +0,0 @@
-/*
-
-------------------------------------------------------------------------------
-
-Standard definitions and types, Bob Jenkins
-
-Modified for inclusion with flam3 by Erik Reckase
-
-------------------------------------------------------------------------------
-
-*/
-
-#ifndef _ISAACS_H_
-
-#define _ISAACS_H_
-
-typedef  unsigned long  int  ub4;   /* unsigned 4-byte quantities */
-
-#define UB4MAXVAL 0xffffffff
-
-typedef    signed long  int  sb4;
-
-#define SB4MAXVAL 0x7fffffff
-
-typedef  unsigned short int  ub2;
-
-#define UB2MAXVAL 0xffff
-
-typedef    signed short int  sb2;
-
-#define SB2MAXVAL 0x7fff
-
-typedef  unsigned       char ub1;
-
-#define UB1MAXVAL 0xff
-
-typedef    signed       char sb1;   /* signed 1-byte quantities */
-
-#define SB1MAXVAL 0x7f
-
-typedef                 int  word;  /* fastest type available */
-
-
-
-#endif
-
diff --git a/rect.c b/rect.c
index d08247d..0b18993 100644
--- a/rect.c
+++ b/rect.c
@@ -859,10 +859,7 @@ static int render_rectangle(flam3_frame *spec, void *out,
 
             int rk;
             /* Create a new isaac state for this thread */
-            for (rk = 0; rk < RANDSIZ; rk++)
-               fth[thi].rc.randrsl[rk] = irand(&spec->rc);
-
-            irandinit(&(fth[thi].rc),1);
+			xorshift_seed (&(fth[thi].rc));
 
             if (0==thi) {
 
diff --git a/xorshift.c b/xorshift.c
new file mode 100644
index 0000000..b0c6e0a
--- /dev/null
+++ b/xorshift.c
@@ -0,0 +1,44 @@
+/*	“An experimental exploration of Marsaglia’s xorshift generators,
+ *	scrambled”, Sebastiano Vigna */
+
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <assert.h>
+#include <unistd.h>
+
+#include "xorshift.h"
+
+uint64_t xorshift_step (randctx * const st) {
+	uint64_t s0 = st->s[ st->p ];
+	uint64_t s1 = st->s[ st->p = ( st->p + 1 ) & (XORSHIFT_S-1) ];
+	s1 ^= s1 << 31; // a
+	s1 ^= s1 >> 11; // b
+	s0 ^= s0 >> 30; // c
+	return ( st->s[ st->p ] = s0 ^ s1 ) * 1181783497276652981LL;
+}
+
+static uint64_t rand64 () {
+	uint64_t rand;
+	while (!__builtin_ia32_rdrand64_step (&rand));
+	return rand;
+}
+
+void xorshift_seed (randctx * const st) {
+	/* seed with high-quality randomness */
+	for (unsigned char i = 0; i < XORSHIFT_S; i++) {
+		st->s[i] = rand64 ();
+	}
+}
+
+#if 0
+uint64_t xorshift_step (randctx * const st) {
+	uint64_t x = st->s[0];
+	x ^= x >> 12; // a
+	x ^= x << 25; // b
+	x ^= x >> 27; // c
+	st->s[0] = x;
+	return x * 2685821657736338717LL;
+}
+#endif
+
diff --git a/xorshift.h b/xorshift.h
new file mode 100644
index 0000000..64c3ea2
--- /dev/null
+++ b/xorshift.h
@@ -0,0 +1,14 @@
+#pragma once
+
+#include <stdint.h>
+
+#define XORSHIFT_S 16
+
+typedef struct {
+	uint64_t s[XORSHIFT_S];
+	int p;
+} randctx;
+
+uint64_t xorshift_step (randctx * const st);
+void xorshift_seed (randctx * const st);
+