Recently I’ve been experimenting with pseudo random number generators (PRNGs), and I’ve been doing this in C. This is a big step for me, as most of my computational work over the past few decades has been in high level computing environments (Matlab, Sage, Maxima etc). And with a PRNG one wants to generate *lots* of values, so that statistical tests can been performed. I am particularly interested in using the TestU01 library, which is coded in C. I was idly wondering whether the following recurrence:

,

with

would make a reasonable PRNG. (It doesn’t: it’s too slow, and fails some tests). But this led me into some programming for testing. First, here it is in C, using the multiple precision GMP library:

/* A program to compute some values x[i] = a^x[i-1] mod p, with a = 7, p = 2^31-1. Compile with: gcc modpowers.c -o modpowers -lgmp -std=gnu99 -g */ #include <stdio.h> #include "gmp.h" #define pp 2147483647 mpz_t y,b,p; static mpz_t x; unsigned long power(void) { mpz_powm(y,b,x,p); // Set y = b^x mod p mpz_set(x,y); // x = y return (mpz_get_ui(x)); } int main() { mpz_init_set_ui(b,7); // Initialize values for use in GMP functions mpz_init_set_ui(p,pp); mpz_init_set_ui(x,7); mpz_init(y); for (int i = 0; i < 1000000; i++) { unsigned long t = power(); if ((i % 100000) == 0) printf("%10i,%10lun",i,t); } return 0; }

And here is its timing, after compiling and running:

$ time ./modpowers 0, 823543 100000,1919783493 200000, 198289041 300000,1263498529 400000, 887186099 500000, 944134154 600000,1803260766 700000, 555710010 800000, 469238760 900000,1834709953 real 0m0.904s user 0m0.900s sys 0m0.000s

So: a million values in less than a second.

Now, because I’m not very sure of my C skills, I often use Sage to check whether I’ve implemented a function properly in C. So here’s Sage doing the same thing:

def printpowers(n): a,p = 7,2^31-1 for i in range(10^n): a = power_mod(7,a,p) if mod(i,10^(n-1))==0: print i,a

We can check that `printpowers(6)` gives the same results as the output of the C program above. But:

sage: timeit('printpowers(6)',repeat=1) 5 loops, best of 1: 73.8 s per loop

As you see, the time is far greater.

Now, I know I can write code in Cython, and use that within Sage, but sometimes it’s nice – especially given some highly optimized and specialized libraries for C – to be able to work directly with C. Anyway, I’m enjoying the huge speed of compiled code.

You might try this sage code (1000000 iterations, though n=100000, 1 second)

def ppowers(n):

p = 2^31-1

a = 7

x = 7

for j in range(10):

for i in range(n):

x = a.powermod(x, p)

print n * (j + 1), x

timeit(‘ppowers(100000)’, repeat=1, number=1)

And thanks for interesting benches!