Did you know that Dr. George Marsaglia (1924-2011), creator of the famed Diehard battery of randomness tests, devised a super-simple PRNG algorithm, just a month prior to his passing? Called the KISS PRNG (or Super-KISS, as there have been multiple KISS versions released previously by Marsaglia), this algorithm boasts a period in excess of 10^40million (10^40,000,000) — an astoundingly large number, many orders of magnitude larger than the famed Mersenne Twister‘s period (only 2^19,937 − 1, or only about 4.3 x 10^6,001 according to gcalctool). Plus, it’s so simple codewise, and also very fast.
Here’s the C implementation (adapted from Marsaglia’s own code):
/* AUTHOR: Shinobu (zuttobenkyou.wordpress.com) */
/* LICENSE: PUBLIC DOMAIN */
#include <stdint.h>
typedef uint64_t u64;
#define QSIZE 0x200000
#define CNG (cng = 6906969069ULL * cng + 13579)
#define XS (xs ^= (xs << 13), xs ^= (xs >> 17), xs ^= (xs << 43))
#define KISS (B64MWC() + CNG + XS)
static u64 QARY[QSIZE];
static int j;
static u64 carry;
static u64 xs;
static u64 cng;
void randk_reset(void)
{
j = QSIZE - 1;
carry = 0;
xs = 362436069362436069ULL;
cng = 123456789987654321ULL; /* use this as the seed */
}
u64 B64MWC(void)
{
u64 t, x;
j = (j + 1) & (QSIZE - 1);
x = QARY[j];
t = (x << 28) + carry;
carry = (x >> 36) - (t < x);
return (QARY[j] = t - x);
}
/* Initialize PRNG with default seed */
void randk_seed(void)
{
u64 i;
/* Seed QARY[] with CNG+XS: */
for (i = 0; i < QSIZE; i++)
QARY[i] = CNG + XS;
}
void randk_seed_manual(u64 seed)
{
cng ^= seed;
xs ^= cng;
randk_seed();
}
void randk_warmup(int rounds)
{
int i;
/* Run through several rounds to warm up the state */
for (i = 0; i < rounds; i++)
randk();
}
/* Generate a pseudorandom 64-bit unsigned integer. */
u64 randk(void)
{
return KISS;
}
Simple, eh? This algorithm is actually 3 PRNG’s in one: the 64-bit Multiply-With-Carry PRNG (B64MWC()), the XOR-Shift PRNG (XS), and the simple Linear Congruential PRNG (CNG). The exorbitant period comes from the fact that this algorithm relies on three different states of the three PRNGs to generate a random number.
Now, where does Haskell come into the picture? Well, I ported the code to Haskell because I wanted a simple PRNG that was of higher quality than the default System.Random RNG. Plus, if you look into the actual source of System.Random, here is an unnerving bit of code:
stdSplit :: StdGen -> (StdGen, StdGen)
stdSplit std@(StdGen s1 s2)
= (left, right)
where
-- no statistical foundation for this!
left = StdGen new_s1 t2
right = StdGen t1 new_s2
new_s1 | s1 == 2147483562 = 1
| otherwise = s1 + 1
new_s2 | s2 == 1 = 2147483398
| otherwise = s2 - 1
StdGen t1 t2 = snd (next std)
See, the RandomGen type class requires the definition of next, split, and genRange functions (see this page). The split function’s purpose is to take one PRNG state, and give two distinct PRNG states, so that you can get multiple unique PRNG’s to work with (this comes up in functional programming in real practice — I speak from experience). The thing is, the statistical robustness of the split function for the StdGen PRNG that comes with Haskell, as can be seen in the source listing, is a bit… annoying/worrying.
Well, when I saw this, I thought: “Hey, why not use KISS? It has 3 PRNGs built into one, so when implementing split, it could just change the state of one of the PRNGs, and you’d get a *completely* different PRNG!” And so that’s exactly what I did:
-- AUTHOR: Shinobu (zuttobenkyou.wordpress.com)
-- LICENSE: PUBLIC DOMAIN
{-# LANGUAGE RecordWildCards #-}
module KISS where
import Data.Array.Unboxed
import Data.Bits
import Data.List
import Data.Word
import System.Random (RandomGen(..))
type U64 = Word64
-- | This is the last KISS-type RNG (2011) that Dr. George Marsaglia (1924-2011)
-- released to the internet before his death. The only difference with this
-- version is that the kissMWCArraySize is 0xfff (4096), instead of 0x200000
-- (2,097,152), for purposes of speed. The period of the original was
-- approximated by Marsaglia as 10^40million, which is practically infinity for
-- most, if not all, needs for everyday programs. The reduced state size for the
-- MWC algorithm from 0x200000 to 0xfff should shorten the period, but it should
-- still be excellent for general usage; because KISS combines not only MWC but
-- also CNG (Congruential) and XSHF (XOR-shift) generators, the period should
-- still be very large.
--
-- TODO: Determine period of this KISS rng.
data KISSRNG = KISSRNG
{ kissMWCArray :: UArray Int U64
, kissMWCArraySize :: U64
, kissMWCIndex :: U64
, kissMWCCarry :: U64
, kissCNG :: U64
, kissXSHF :: U64
}
kissStateSize :: U64
kissStateSize = 0xfff
roundCNG :: U64 -> U64
roundCNG cng = 6906969069 * cng + 13579
roundXSHF :: U64 -> U64
roundXSHF = round3 . round2 . round1
where
round1 b = xor b (unsafeShiftL b 13)
round2 b = xor b (unsafeShiftR b 17)
round3 b = xor b (unsafeShiftL b 43)
roundB64MWC :: KISSRNG -> KISSRNG
roundB64MWC kiss@KISSRNG{..} = kiss
{ kissMWCArray = array'
, kissMWCIndex = index'
, kissMWCCarry = carry'
}
where
index' = (kissMWCIndex + 1) .&. (kissMWCArraySize - 1)
x = kissMWCArray ! (fromIntegral index')
t = unsafeShiftL x 28 + kissMWCCarry
carry' = unsafeShiftR x 36 - (if t < x then 1 else 0)
array' = kissMWCArray // [(fromIntegral index', t - x)]
makeKISSRNG :: U64 -> KISSRNG
makeKISSRNG seed = KISSRNG
{ kissMWCArray = array (0, (fromIntegral kissStateSize - 1)) $ zip [0..] kissArray
, kissMWCArraySize = kissStateSize
, kissMWCIndex = kissStateSize - 1
, kissMWCCarry = 0
, kissCNG = cngWarmed
, kissXSHF = xshfWarmed
}
where
-- seed the MWC array with the Congruential and XOR-Shift RNG's
(kissArray, cngWarmed, xshfWarmed) = foldl' step ([], seedCNG, seedXSHF) [0..(kissStateSize - 1)]
step (ary, cng, xshf) _ = ((cng' + xshf'):ary, cng', xshf')
where
cng' = roundCNG cng
xshf' = roundXSHF xshf
-- default Congruential RNG seed
seedCNG = 123456789987654321
-- default XOR-Shift RNG seed
seedXSHF = 362436069362436069
randKISS :: KISSRNG -> (U64, KISSRNG)
randKISS k = (kissMWC + cng + xshf, k' {kissCNG = cng, kissXSHF = xshf})
where
k' = roundB64MWC k
kissMWC = (kissMWCArray k') ! (fromIntegral $ kissMWCIndex k')
cng = roundCNG $ kissCNG k
xshf = roundXSHF $ kissXSHF k
instance Show KISSRNG where
show KISSRNG{..} = "kissMWC: [kissMWC..]"
++ "\nkissMWCIdx: " ++ show kissMWCIndex ++ "\n"
++ "kissMWCArray!!Idx: " ++ show (kissMWCArray ! (fromIntegral kissMWCIndex)) ++ "\n"
++ "kissMWCArray!!Head: " ++ show (kissMWCArray ! 0) ++ "\n"
++ "kissMWCArray!!Last: " ++ show (kissMWCArray ! (fromIntegral $ kissStateSize - 1)) ++ "\n"
++ "kissMWCCarry: " ++ show kissMWCCarry ++ "\n"
++ "kissCNG: " ++ show kissCNG ++ "\n"
++ "kissXSHF: " ++ show kissXSHF ++ "\n"
instance RandomGen KISSRNG where
next rng = (fromIntegral n, rng')
where
(n, rng') = randKISS rng
split rng = (rng1, rng2)
where
rng1 = warmupXSHF rng
rng2 = warmupMWC rng
genRange _ = (0, 0xffffffffffffffff)
warmupRNG :: RandomGen g => g -> Int -> g
warmupRNG g rounds = foldl' warmup g [1..rounds]
where
warmup g' _ = snd $ next g'
warmupMWC :: KISSRNG -> KISSRNG
warmupMWC rng = roundB64MWC rng
warmupXSHF :: KISSRNG -> KISSRNG
warmupXSHF rng = rng { kissXSHF = roundXSHF $ kissXSHF rng}
In the implementation of split, you can clearly see that we simply warm up one of PRNGs (move on to the next state in the period) to get a new PRNG. Again, since KISS depends on all three PRNGs, simply changing the state of one of the PRNGs will give you a completely different PRNG.
Oh, the only weakness of the Haskell version is that its QSIZE is only 0xfff, not 0x200000 as in the original, for performance reasons. I certainly hope someone could improve the performance of the code and release it on Hackage (my code is hereby released into the PUBLIC DOMAIN, so do what you like with it), restoring the state size to 0x200000 as in Marsaglia’s original. Also, I’m not sure how large the period is, but judging by how the XOR-Shift PRNG has a large period on its own, it should still be very, very large with a 0xfff state size for the MWC algorithm.
I would sincerely appreciate it if someone more familiar with combinatorics/compsci could tell me what the size of the period is with a 0xfff state size for the MWC.
I was also pleasantly surprised by the very good quality of KISS. I used my code to write some random bits into a file, and used the ent program to judge the entroy of it. Here are the results:
Entropy = 7.999829 bits per byte.
Optimum compression would reduce the size
of this 1048576 byte file by 0 percent.
Chi square distribution for 1048576 samples is 248.29, and randomly
would exceed this value 60.65 percent of the times.
Arithmetic mean value of data bytes is 127.5231 (127.5 = random).
Monte Carlo value for Pi is 3.141895835 (error 0.01 percent).
Serial correlation coefficient is 0.001437 (totally uncorrelated = 0.0).
The results show that the KISS RNG has excellent quality random numbers. These figures make it as good (randomness-wise) as, e.g., the one based on AES encryption (AES in counter mode), which has also been analyzed with ent, as stated on the github page:
Using ent, a randomness property maker on one 1Mb sample.
cprng-AES:
Entropy = 7.999837 bits per byte.
Optimum compression would reduce the size of this 1048576 byte file by 0 percent.
Chi square distribution for 1048576 samples is 237.02.
Arithmetic mean value of data bytes is 127.3422 (127.5 = random).
Monte Carlo value for Pi is 3.143589568 (error 0.06 percent).
The rather ugly code I used to generate this file (and believe me, it took forever to generate a 1MiB file because the code is horribly unoptimized…) is below:
-- AUTHOR: Shinobu (zuttobenkyou.wordpress.com)
-- LICENSE: PUBLIC DOMAIN
{-# LANGUAGE RecordWildCards #-}
module Main where
import Data.Array.Unboxed
import Data.Bits
import Data.List
import Data.Word
import System.Random (RandomGen(..))
import qualified Data.ByteString as BS
type U64 = Word64
data KISSRNG = KISSRNG
{ kissMWCArray :: UArray Int U64
, kissMWCArraySize :: U64
, kissMWCIndex :: U64
, kissMWCCarry :: U64
, kissCNG :: U64
, kissXSHF :: U64
}
kissStateSize :: U64
kissStateSize = 0xfff
roundCNG :: U64 -> U64
roundCNG cng = 6906969069 * cng + 13579
roundXSHF :: U64 -> U64
roundXSHF = round3 . round2 . round1
where
round1 b = xor b (unsafeShiftL b 13)
round2 b = xor b (unsafeShiftR b 17)
round3 b = xor b (unsafeShiftL b 43)
roundB64MWC :: KISSRNG -> KISSRNG
roundB64MWC kiss@KISSRNG{..} = kiss
{ kissMWCArray = array'
, kissMWCIndex = index'
, kissMWCCarry = carry'
}
where
index' = (kissMWCIndex + 1) .&. (kissMWCArraySize - 1)
x = kissMWCArray ! (fromIntegral index')
t = unsafeShiftL x 28 + kissMWCCarry
carry' = unsafeShiftR x 36 - (if t < x then 1 else 0)
array' = kissMWCArray // [(fromIntegral index', t - x)]
makeKISSRNG :: U64 -> KISSRNG
makeKISSRNG seed = KISSRNG
{ kissMWCArray = array (0, (fromIntegral kissStateSize - 1)) $ zip [0..] kissArray
, kissMWCArraySize = kissStateSize
, kissMWCIndex = kissStateSize - 1
, kissMWCCarry = 0
, kissCNG = xor cngWarmed seed
, kissXSHF = xor xshfWarmed seed
}
where
-- seed the MWC array with the Congruential and XOR-Shift RNG's
(kissArray, cngWarmed, xshfWarmed) = foldl' step ([], seedCNG, seedXSHF) [0..(kissStateSize - 1)]
step (ary, cng, xshf) _ = ((cng' + xshf'):ary, cng', xshf')
where
cng' = roundCNG cng
xshf' = roundXSHF xshf
-- default Congruential RNG seed
seedCNG = 123456789987654321
-- default XOR-Shift RNG seed
seedXSHF = 362436069362436069
randKISS :: KISSRNG -> (U64, KISSRNG)
randKISS k = (kissMWC + cng + xshf, k' {kissCNG = cng, kissXSHF = xshf})
where
k' = roundB64MWC k
kissMWC = (kissMWCArray k') ! (fromIntegral $ kissMWCIndex k')
cng = roundCNG $ kissCNG k
xshf = roundXSHF $ kissXSHF k
instance Show KISSRNG where
show KISSRNG{..} = "kissMWC: [kissMWC..]"
++ "\nkissMWCIdx: " ++ show kissMWCIndex ++ "\n"
++ "kissMWCArray!!Idx: " ++ show (kissMWCArray ! (fromIntegral kissMWCIndex)) ++ "\n"
++ "kissMWCArray!!Head: " ++ show (kissMWCArray ! 0) ++ "\n"
++ "kissMWCArray!!Last: " ++ show (kissMWCArray ! (fromIntegral $ kissStateSize - 1)) ++ "\n"
++ "kissMWCCarry: " ++ show kissMWCCarry ++ "\n"
++ "kissCNG: " ++ show kissCNG ++ "\n"
++ "kissXSHF: " ++ show kissXSHF ++ "\n"
instance RandomGen KISSRNG where
next rng = (fromIntegral n, rng')
where
(n, rng') = randKISS rng
split rng = (rng1, rng2)
where
rng1 = warmupXSHF rng
rng2 = warmupMWC rng
genRange _ = (0, 0xffffffffffffffff)
warmupRNG :: RandomGen g => g -> Int -> g
warmupRNG g rounds = foldl' warmup g [1..rounds]
where
warmup g' _ = snd $ next g'
warmupMWC :: KISSRNG -> KISSRNG
warmupMWC rng = roundB64MWC rng
warmupXSHF :: KISSRNG -> KISSRNG
warmupXSHF rng = rng { kissXSHF = roundXSHF $ kissXSHF rng}
main :: IO ()
main = do
let
rng = makeKISSRNG 0
(bytes1MiB, _) = genBytesKISS 0x100000 rng
BS.writeFile "data" bytes1MiB
genBytesKISS :: U64 -> KISSRNG -> (BS.ByteString, KISSRNG)
genBytesKISS len kissrng = foldl' step (BS.empty, kissrng) [1..(div len 8)] -- divide by 8, b/c e.g., to generate 8 bytes, we only need 1 U64
where
step (bs, rng) _ = (foldl' BS.snoc bs $ octets u64, rng')
where
(u64, rng') = randKISS rng
smallerChunks :: [U64] -> [Word8]
smallerChunks = concatMap octets
-- | Get a number and split it up into 8 8-bit parts (64 bits total).
octets :: (Bits a, Integral a) => a -> [Word8]
octets w = map (\n -> fromIntegral $ shiftR w n) . reverse $ take 8 [0,8..]
Here are the C and Haskell standalone versions that prove that the Haskell port behaves in the same way as the C version, given the right starting seed values and state size (both 0xfff for the MWC PRNG):
C standalone version (compile with gcc -o ckiss kiss.c):
/* AUTHOR: Shinobu (zuttobenkyou.wordpress.com) */
/* LICENSE: PUBLIC DOMAIN */
#include <stdio.h>
#include <inttypes.h>
#include <stdint.h>
typedef uint64_t u64;
#define QSIZE 0xfff
#define CNG (cng = 6906969069ULL * cng + 13579)
#define XS (xs ^= (xs << 13), xs ^= (xs >> 17), xs ^= (xs << 43))
#define KISS (B64MWC() + CNG + XS)
static u64 QARY[QSIZE];
static int j;
static u64 carry;
static u64 xs;
static u64 cng;
u64 B64MWC(void)
{
u64 t, x;
j = (j + 1) & (QSIZE - 1);
x = QARY[j];
t = (x << 28) + carry;
carry = (x >> 36) - (t < x);
return (QARY[j] = t - x);
}
/* Initialize PRNG with default seed */
void randk_seed(void)
{
u64 i;
j = QSIZE - 1;
carry = 0;
xs = 362436069362436069ULL;
cng = 123456789987654321ULL;
/* Seed QARY[] with CNG+XS: */
for (i = 0; i < QSIZE; i++)
QARY[i] = CNG + XS;
}
/* Generate a pseudorandom 64-bit unsigned integer. */
u64 randk(void)
{
return KISS;
}
int main(void)
{
randk_seed();
printf("randk_seed called!\n");
printf("KISS idx: %"PRIu64"\n", j);
printf("qary[idx] is: %"PRIu64"\n", QARY[j]);
printf("qary[0] is: %"PRIu64"\n", QARY[0]);
printf("qary[QSIZE - 1] is: %"PRIu64"\n", QARY[QSIZE - 1]);
printf("carry: %"PRIu64"\n", carry);
printf("cng: %"PRIu64"\n", cng);
printf("xs: %"PRIu64"\n", xs);
u64 x = KISS;
printf("\nKISS called! KISS num is: %"PRIu64"\n", x);
printf("\nKISS idx: %"PRIu64"\n", j);
printf("qary[idx] is: %"PRIu64"\n", QARY[j]);
printf("qary[0] is: %"PRIu64"\n", QARY[0]);
printf("qary[QSIZE - 1] is: %"PRIu64"\n", QARY[QSIZE - 1]);
printf("carry: %"PRIu64"\n", carry);
printf("cng: %"PRIu64"\n", cng);
printf("xs: %"PRIu64"\n", xs);
printf("x + 18334599312639636657 is: %"PRIu64"\n", x + 18334599312639636657ULL);
}
Haskell standalone version (run with runhaskell kiss.hs):
-- AUTHOR: Shinobu (zuttobenkyou.wordpress.com)
-- LICENSE: PUBLIC DOMAIN
{-# LANGUAGE RecordWildCards #-}
module KISS where
import Data.Array.Unboxed
import Data.Bits
import Data.List
import Data.Word
import System.Random (RandomGen(..))
type U64 = Word64
data KISSRNG = KISSRNG
{ kissMWCArray :: UArray Int U64
, kissMWCArraySize :: U64
, kissMWCIndex :: U64
, kissMWCCarry :: U64
, kissCNG :: U64
, kissXSHF :: U64
}
kissStateSize :: U64
kissStateSize = 0xfff
roundCNG :: U64 -> U64
roundCNG cng = 6906969069 * cng + 13579
roundXSHF :: U64 -> U64
roundXSHF = round3 . round2 . round1
where
round1 b = xor b (unsafeShiftL b 13)
round2 b = xor b (unsafeShiftR b 17)
round3 b = xor b (unsafeShiftL b 43)
roundB64MWC :: KISSRNG -> KISSRNG
roundB64MWC kiss@KISSRNG{..} = kiss
{ kissMWCArray = array'
, kissMWCIndex = index'
, kissMWCCarry = carry'
}
where
index' = (kissMWCIndex + 1) .&. (kissMWCArraySize - 1)
x = kissMWCArray ! (fromIntegral index')
t = unsafeShiftL x 28 + kissMWCCarry
carry' = unsafeShiftR x 36 - (if t < x then 1 else 0)
array' = kissMWCArray // [(fromIntegral index', t - x)]
makeKISSRNG :: U64 -> KISSRNG
makeKISSRNG seed = KISSRNG
{ kissMWCArray = array (0, (fromIntegral kissStateSize - 1)) $ zip [0..] kissArray
, kissMWCArraySize = kissStateSize
, kissMWCIndex = kissStateSize - 1
, kissMWCCarry = 0
, kissCNG = cngWarmed
, kissXSHF = xshfWarmed
}
where
-- seed the MWC array with the Congruential and XOR-Shift RNG's
(kissArray, cngWarmed, xshfWarmed) = foldl' step ([], seedCNG, seedXSHF) [0..(kissStateSize - 1)]
step (ary, cng, xshf) _ = ((cng' + xshf'):ary, cng', xshf')
where
cng' = roundCNG cng
xshf' = roundXSHF xshf
-- default Congruential RNG seed
seedCNG = 123456789987654321
-- default XOR-Shift RNG seed
seedXSHF = 362436069362436069
randKISS :: KISSRNG -> (U64, KISSRNG)
randKISS k = (kissMWC + cng + xshf, k' {kissCNG = cng, kissXSHF = xshf})
where
k' = roundB64MWC k
kissMWC = (kissMWCArray k') ! (fromIntegral $ kissMWCIndex k')
cng = roundCNG $ kissCNG k
xshf = roundXSHF $ kissXSHF k
instance Show KISSRNG where
show KISSRNG{..} = "kissMWC: [kissMWC..]"
++ "\nkissMWCIdx: " ++ show kissMWCIndex ++ "\n"
++ "kissMWCArray!!Idx: " ++ show (kissMWCArray ! (fromIntegral kissMWCIndex)) ++ "\n"
++ "kissMWCArray!!Head: " ++ show (kissMWCArray ! 0) ++ "\n"
++ "kissMWCArray!!Last: " ++ show (kissMWCArray ! (fromIntegral $ kissStateSize - 1)) ++ "\n"
++ "kissMWCCarry: " ++ show kissMWCCarry ++ "\n"
++ "kissCNG: " ++ show kissCNG ++ "\n"
++ "kissXSHF: " ++ show kissXSHF ++ "\n"
instance RandomGen KISSRNG where
next rng = (fromIntegral n, rng')
where
(n, rng') = randKISS rng
split rng = (rng1, rng2)
where
rng1 = warmupXSHF rng
rng2 = warmupMWC rng
genRange _ = (0, 0xffffffffffffffff)
warmupRNG :: RandomGen g => g -> Int -> g
warmupRNG g rounds = foldl' warmup g [1..rounds]
where
warmup g' _ = snd $ next g'
warmupMWC :: KISSRNG -> KISSRNG
warmupMWC rng = roundB64MWC rng
warmupXSHF :: KISSRNG -> KISSRNG
warmupXSHF rng = rng { kissXSHF = roundXSHF $ kissXSHF rng}
main :: IO ()
main = do
let rng = makeKISSRNG 0
putStrLn $ show rng
EDIT May 3, 2012: Sorry about the somewhat redundant-looking code listings, but I’m too lazy/busy to clean them up.
EDIT May 4, 2012: Alternate period number for the Mersenne Twister, for easier comparison of period size with KISS.
Shinobu's Secrets 