18

On-the-fly Linear Congruential Generator Using Emacs Calc

 4 years ago
source link: https://nullprogram.com/blog/2019/11/19/
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.

I regularly make throwaway “projects” and do a surprising amount of programming in /tmp . For Emacs Lisp, the equivalent is the *scratch* buffer. These are places where I can make a mess, and the mess usually gets cleaned up before it becomes a problem. A lot of my established projects (ex.) start out in volatile storage and only graduate to more permanent storage once the concept has proven itself.

Throughout my whole career, this sort of throwaway experimentation has been an important part of my personal growth, and I try to encourage it in others . Even if the idea I’m trying doesn’t pan out, I usually learn something new, and occasionally it translates into an article here.

I also enjoy small programming challenges. One of the most abused tools in my mental toolbox is the Monte Carlo method, and I readily apply it to solve toy problems. Even beyond this, random number generators are frequently a useful tool (1,2), so I find myself reaching for one all the time.

Virtually every programming language comes with a pseudo-random number generation function or library. Unfortunately the language standard PRNG often not ideal for my purposes. It’s probably not high quality, slower than it needs to be ( also ), lacks reliable semantics or behavior between implementations , or is missing some other property I want. So I’ve long been a fan of BYOPRNG: Bring Your Own Pseudo-random Number Generator. Just embed a generator with the desired properties directly into the program. The best non-cryptographic PRNGs today are tiny and exceptionally friendly to embedding. Though, depending on what you’re doing, you might need to be creative about seeding .

Crafting a PRNG

On occasion I don’t have an established, embeddable PRNG in reach, and I have yet to commit xoshiro256** to memory. Or maybe I want to use a totally unique PRNG for a particular project. In these cases I make one up. With just a bit of know-how it’s not too difficult.

Probably the easiest decent PRNG to code from scratch is the venerable Linear Congruential Generator (LCG). It’s a simple recurrence relation:

x[1] = (x[0] * A + C) % M

That’s trivial to remember once you know the details. You only need to choose appropriate values for A , C , and M . Done correctly, it will be a full-period generator — a generator that visits a permutation of each of the numbers between 0 and M - 1 . The seed — the value of x[0] — is chooses a starting position in this (looping) permutation.

M has a natural, obvious choice: a power of two matching the range of operands, such as 2^32 or 2^64. With this the modulo operation is free as a natural side effect of the computer architecture.

Choosing C also isn’t difficult. It must be co-prime with M , and since M is a power of two, any odd number is valid. Even 1. In theory choosing a small value like 1 is faster since the compiler won’t need to embed a large integer in the code, but this difference doesn’t show up in any micro-benchmarks I tried. If you want a cool, unique generator, then choose a large random integer. More on that below.

The tricky value is A , and getting it right is the linchpin of the whole LCG. It must be coprime with M (i.e. not even), and, for a full-period generator, A-1 must be divisible by four. For better results, A-1 should not be divisible by 8. A good choice is a prime number that satisfies these properties.

If your operands are 64-bit integers, or larger, how are you going to generate a prime number?

Primes from Emacs Calc

Emacs Calc can solve this problem. I’venoted before how featureful it is. It has arbitrary precision, random number generation, and primality testing. It’s everything we need to choose A . (In fact, this is nearly identical to the process I used to implement RSA .) For this example I’m going to generate a 64-bit LCG for the C programming language, but it’s easy to use whatever width you like and mostly whatever language you like. If you wanted a minimal standard 128-bit LCG , this will still work.

Start by opening up Calc with M-x calc , then:

  1. Push 2 on the stack
  2. Push 64 on the stack
  3. Press ^ , computing 2^64 and pushing it on the stack
  4. Press k r to generate a random number in this range
  5. Press d r 16 to switch to hexadecimal display
  6. Press k n to find the next prime following the random value
  7. Repeat step 6 until you get a number that ends with 5 or D
  8. Press k p a few times to avoid false positives.

What’s left on the stack is your A ! If you want a random value for C , you can follow a similar process. Heck, make it prime, too!

The reason for using hexadecimal (step 5) and looking for 5 or D (step 7) is that such numbers satisfy both of the important properties for A-1 .

Calc doesn’t try to factor your random integer. Instead it uses the Miller–Rabin primality test , a probabilistic test that, itself, requires random numbers. It has false positives but no false negatives. The false positives can be mitigated by repeating the test multiple times, hence step 8.

Trying this all out right now, I got this implementation (in C):

uint64_t lcg1(void)
{
    static uint64_t s = 0;
    s = s*UINT64_C(0x7c3c3267d015ceb5) + UINT64_C(0x24bd2d95276253a9);
    return s;
}

However, we can still do a little better. Outputting the entire state doesn’t have great results, so instead it’s better to create a truncated LCG and only return some portion of the most significant bits.

uint32_t lcg2(void)
{
    static uint64_t s = 0;
    s = s*UINT64_C(0x7c3c3267d015ceb5) + UINT64_C(0x24bd2d95276253a9);
    return s >> 32;
}

This won’t quite pass BigCrush in 64-bit form, but the results are pretty reasonable for most purposes.

But we can still do better without needing to remember much more than this.

Appending permutation

A Permuted Congruential Generator (PCG) is really just a truncated LCG with a permutation applied to its output. Like LCGs themselves, there are arbitrarily many variations. The “official” implementation has adata-dependent shift, for which I can never remember the details. Fortunately a couple of simple, easy to remember transformations is sufficient. Basically anything I used while prospecting for hash functions . I love xorshifts, so lets add one of those:

uint32_t pcg1(void)
{
    static uint64_t s = 0;
    s = s*UINT64_C(0x7c3c3267d015ceb5) + UINT64_C(0x24bd2d95276253a9);
    uint32_t r = s >> 32;
    r ^= r >> 16;
    return r;
}

This is a big improvement, but it still fails one BigCrush test. As they say, when xorshift isn’t enough, use xorshift-multiply! Below I generated a 32-bit prime for the multiply, but any odd integer is a valid permutation.

uint32_t pcg2(void)
{
    static uint64_t s = 0;
    s = s*UINT64_C(0x7c3c3267d015ceb5) + UINT64_C(0x24bd2d95276253a9);
    uint32_t r = s >> 32;
    r ^= r >> 16;
    r *= UINT32_C(0x60857ba9);
    return r;
}

This passes BigCrush, and I can reliably build a new one entirely from scratch using Calc any time I need it.

Bonus: Adapting to other languages

Sometimes it’s not so straightforward to adapt this technique to other languages. For example, JavaScript has limited support for 32-bit integer operations (enough for a poor 32-bit LCG) and no 64-bit integer operations. Though BigInt is now a thing, and should make a great 96- or 128-bit LCG easy to build.

Java doesn’t have unsigned integers, so how could you build the above PCG in Java? Easy! First, remember is that Java has two’s complement semantics, including wrap around, and that two’s complement doesn’t care about unsigned or signed for multiplication (or addition, or subtraction). The result is identical. Second, the oft-forgotten >>> operator does an unsigned right shift. With these two tips:

long s = 0;

int pcg2() {
    s = s*0x7c3c3267d015ceb5L + 0x24bd2d95276253a9L;
    int r = (int)(s >>> 32);
    r ^= r >>> 16;
    r *= 0x60857ba9;
    return r;
}

So, in addition to the Calc step list above, you may need to know some of the finer details of your target language.


About Joyk


Aggregate valuable and interesting links.
Joyk means Joy of geeK