[nvvsbkrb] Primes with primitive root 2

The motivation for this exercise is a screen transition of the "pedal" screensaver, which plots a series of black pixels "randomly" until the screen is blank (actually, it gives up before that). Suppose we wish to have a "random" (this is a graphics application; it need only look random to the eye) permutation of all the pixels on the screen.

We generate pixel coordinates (linearized 1..(m*n)) with the recurrence c <- 2*c mod p, where p is chosen to be slightly larger than m*n, then skipping c if it is larger than m*n. We require primes with primitive root 2. We use 2 because, up to 2^31-1, we can do modular multiplication without overflow in 32 bits.

Essentially, we run through the whole period of a (probably bad) linear congruential random number generator.

Here is the next suitable prime (the delta to the next prime) after certain common screen resolutions.

*Main> next_primitive_2(640*480)
43
*Main> next_primitive_2(800*600)
19
*Main> next_primitive_2(1024*768)
37
*Main> next_primitive_2(1280*1024)
21
*Main> next_primitive_2(1600*1200)
11
*Main> next_primitive_2(1280*768)
29
*Main> next_primitive_2(1366*768)
5
*Main> next_primitive_2(1280*800)
21
*Main> next_primitive_2(1440*900)
11
*Main> next_primitive_2(1680*1050)
53
*Main> next_primitive_2(1920*1200)
11
*Main> next_primitive_2(1920*1080)
13
*Main> next_primitive_2(2560*1600)
13

For example, VGA (640 x 480) generates the sequence [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768, 65536, 131072, 262144, 217045, 126847, 253694, 200145, 93047, 186094] where the modulus is 640*480+43 = 307243 .

Animation of 128x96 (modulus is 128*96+13), magnified by a factor of 2. code

Exercising Pari/GP to see how many there are.

? default(timer,1)
? default(primelimit,2^31+10)
time = 5,120 ms.
%2 = 2147483658
? for(i=1,31,z=0;forprime(p=2^(i-1)+1,2^i,z=z+isprimroot(2,p));print(i," ",z))

1 1
2 1
3 1
4 2
5 2
6 4
7 4
8 10
9 15
10 29
11 51
12 101
13 175
14 316
15 601
16 1131
17 2148
18 3993
19 7635
20 14427
21 27566
22 52564
23 100322
24 192059
25 368570
26 708620
27 1363431
28 2628113
29 5071530
30 9800737
31 18957386               
time = 47mn, 47,672 ms.

Some Haskell code, including factoring by trial division. Note that "a" needs to be at least Int64 in order to handle up to 2^32. Assuming we are interested only in small integers, how much faster can this routine become? So many opportunities for clever algorithmic and number theoretic tricks.

Hardcoding 39301545 (or 39301544 not including 2) integers is not a great option. (Sum of the numbers above.)

next_primitive_2 :: forall a . (Integral (a)) => a -> a; next_primitive_2 x = (let { divide_test_seq :: a -> [](a) -> Maybe([](a)); divide_test_seq n xs = (case xs of { ((:)(x) (xrest))-> (case ((<) n ((*) x x)) of { (True)-> Nothing; (_)-> (case (mod n x) of { (0)-> (Just xs); (_)-> (divide_test_seq n xrest) }) }) }); factor :: a -> [](a); factor n = (factor_with_seq n prime_seq); factor_with_seq :: a -> [](a) -> [](a); factor_with_seq n nums = (case (compare n 1) of { (EQ)-> []; (GT)-> (case (divide_test_seq n nums) of { (Nothing)-> [n]; (Just(xs))-> ((:) (head xs) (factor_with_seq (div n (head xs)) xs)) }) }); is_prime :: a -> Bool; is_prime n = ((&&) ((<) 1 n) (isNothing(trial_division(n)))); is_primitive :: a -> a -> Bool; is_primitive x m = (and((map ((/=) 1))((test_prim_list m x)))); mod_exp :: a -> a -> a -> a; mod_exp m x n = (case (compare n 1) of { (EQ)-> x; (GT)-> (case (divMod n 2) of { ((q), (0))-> (mod_square m (mod_exp m x q)); ((q), (1))-> (mod_mult m x (mod_square m (mod_exp m x q))) }) }); mod_mult :: a -> a -> a -> a; mod_mult modulus x y = (mod ((*) x y) modulus); mod_square :: a -> a -> a; mod_square m x = (mod_mult m x x); prime_seq :: [](a); prime_seq = ((:) 2 (enumFromThen 3 5)); test_prim_list :: a -> a -> [](a); test_prim_list m x = (map (mod_exp m x) (unique_factors_div (pred m))); trial_division :: a -> Maybe([](a)); trial_division n = (divide_test_seq n prime_seq); unique_factors_div :: a -> [](a); unique_factors_div n = (map (div n) (map head (group (factor n)))) } in ((-) (head((filter (is_primitive 2))((filter is_prime)((enumFrom x))))) x))