Number theoretic computations in C

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:

x_i = b^{x_{i-1}}pmod{p}, u_i=x_i/p

with

b = 7,quad p=2^{31}-1,quad x_0 = 7

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.

One thought on “Number theoretic computations in C

  1. 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!

Leave a Reply

Your email address will not be published. Required fields are marked *