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

,

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 &lt;stdio.h&gt;
#include &quot;gmp.h&quot;
#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 &lt; 1000000; i++)
{
unsigned long t = power();
if ((i % 100000) == 0)
printf(&quot;%10i,%10lun&quot;,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. aam says:

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!