View difference between Paste ID: tbGxRwPB and wyFMP7R6
SHOW: | | - or go back to the newest paste.
1
#include <stdio.h>
2
#include <stdint.h>
3
#include <stdlib.h>
4
5
/* http://marc-b-reynolds.github.io/math/2017/09/18/ModInverse.html */
6
static uint32_t modinv(uint32_t a)
7
{
8
	uint32_t x;
9
	x = a;		//  3 bits
10
	x *= 2 - a*x;	//  6
11
	x *= 2 - a*x;	// 12
12
	x *= 2 - a*x;	// 24
13
	x *= 2 - a*x;	// 48 / retaining bottom 32
14
	return x;
15
}
16
17
/* http://number-none.com/blow/blog/programming/2016/07/08/fabian-on-lcg-fast-forward.html */
18
uint32_t skip(uint32_t x, uint32_t a, uint32_t c, unsigned n)
19
{
20
	while (n) {
21
    		if (n & 1)
22
        		x = a * x + c;
23
		c *= a + 1;
24
		a *= a;
25
		n >>= 1;
26
	}
27
	return x;
28
}
29
30
int main(int argc, char **argv)
31
{
32
	srandom(time(NULL));
33
	unsigned i, s = 1, a = random() | 1UL, c =  random() | 1UL;
34
35
	/* move forward one step at a time */
36
	for (i = 0; i < 100; i++)
37
		s = s * a + c;
38
	printf("%08x\n", s);
39
	/* skip forwaward */
40
	printf("%08x\n", skip(1, a, c, 100));
41
	/* you could loop through the full LCG cycle */
42
	printf("%08x\n", skip(s, a, c, -100));
43
	/* or inverse of LCG is still LCG */
44
	a = modinv(a);
45
	c = -c * a;
46
	/* skip backward */
47
	printf("%08x\n", skip(s, a, c, 100));
48
	/* move backward one step at a time */
49
	for (i = 0; i < 100; i++)
50
		s = s * a + c;
51
	printf("%08x\n", s);
52
}