In my discussion of efficient computation of Fibonacci's iterator, I exploited a subtle trick we can play to do it faster by expressing the iteration as multiplication by one fixed polynomial modulo another. The initial state is then described by a polynomial and all later states are this one times powers of the iterator polynomial, modulo the invariant one. Since computing powers can be done in logarithmic time, this works nicely for getting sparse values from the iterator, by bypassing the need to compute lots of intermediate states. I naively thought I could apply the same to any other iteratively defined sequence.
My friend Colin wrote about computing Perrin Pseudo-Primes: these use an iterator defined by k(1) = 0, k(2) = 2, k(3) = 3 and k(n+3) = k(n+1) +k(n), from which we can back-track k(0) = 3 if we want (: k |{naturals}). The significance of this sequence is that n divides k(n) whenever n is prime and usually the converse works, too; i.e. if n divides k(n), then n is quite likely prime. An iterative solution wants to keep track of a triple [k(i),k(i+1),k(i+2)] and infer from it [k(i+1),k(i+2),k(i+3)], which requires an iteration step [b,c,a+b] ← [a,b,c]. The successive steps after a [1,0,0] triple are then [0,0,1], [0,1,0], [1,0,1], [0,1,1], [1,1,1], [1,1,2], [1,2,2], [2,2,3] and so on. If we take 3 = k(0) of each plus 2 = k(2) of the next (and 0 = k(1) of the one after that), we get our [k(i),k(i+1),k(i+2)] triples for i = 0, 1, … Notice, however, the peculiarity that k(i).[1,0,0] +k(i+1).[0,1,0] +k(i+2).[0.0,1] has its three unit triples in a different order than the iteration, which maps the first to the third to the second.
I naïvely (and, with hindsight, obviously foolishly) supposed I could represent a triple [a,b,c] by a +b.x +c.x.x (where x stands for the polynomial power(1)) and working modulo x.x.x −x −1 would make multiplication by x be my iterator; this is the simple analogue of what worked for Fibonacci. However, the iterator needs to get b +c.x +(a+b).x.x at its next step and instead we get a.x +b.x.x +c.(x +1) = c +(a+c).x +b.x.x, so it's clearly not that simple ! Perhaps I just need a different multiplier ? Try
in which we want u.a +v.c +w.b = b, so w = 1, u = v = 0, so we need to multiply by x.x rather than x; however, we also want v.a +(u +w).b +(v +w).c = c, so v = 0, u = −w and w = 1 −v = 1 (again), giving x +x.x as factor; and we want w.a +v.b +(u +w).c = a +b, so w = 1 (again) but v = 1 and u = −w = −1. So there's no factor that works modulo x.x.x = x +1, with the given representation of a triple. Clearly I was lucky with the Fibonacci iterator in some way, so let's work out what we actually need.
Maybe I just need to use my triple in some other order as polynomial coefficients. Each entry in the iterated triple has zero coefficient of the corresponding entry in the input triple; this gives us u = 0 from the constant term's coefficient of a and u+w = 0 from the other two terms; so we have u = 0 = w and must have v = 1, so our scaling is indeed the polynomial x and we map a +b.x +c.x.x to c +(a+c).x +b.x.x; so the right way to represent our iteration [b,c,a+b] ←[a,b,c] is to identify the input triple [a,b,c] with a +c.x +b.x.x; scaling this by x gives us a.x +c.x.x +b.(1 +x) = b +(a+b).x +c.x.x which is indeed what we represent as [b,c,a+b]. Remember that I mentioned we get the unit triples by which we multiply the k(i) in a twisted order relative to the iteration ? We're here doing the same thing for the polynomial coefficients.
So our triple [k(i),k(i+1),k(i+2)] gets represented as k(i) +k(i+2).x
+k(i+1).x.x and repeated multiplication by x performs our iteration, starting
with the polynomial 3 +2.x representing [3,0,2]. After reducing modulo x.x.x =
x +1, the coefficient of x.x in (3 +2.x).power(i, x) is k(i+1) for each natural
i. We just need to define multiply modulo x.x.x = x +1
and use it to
compute powers of x times a fixed item. Of course, since a polynomial modulo a
cubic is represented by just three coefficients, it's natural to use our triples
as the actual data by which we represent the cubic and express the
multiplication on polynomials modulo x.x.x = x +1 as a multiplication on
triples; which, aside from having the entries in the wrong order, I already did
above. Correcting the order, we now get [a,b.c].[u.v.w] = [a.u +b.w +c.v, a.v
+c.w +b.u +b.v, a.w +c.u +c.v +b.w +b.v] and can use this
in the usual O(log(n)) way to
compute power(n) to compute k(n) % n for any given natural n. Naturally, I've
implemented this in my pythonic study package
in
module study.maths.Perrin
as a class Perrin, whose .primal(n) method tests a natural
n for whether k(n) is a multiple of n.
We can also wind this iterator backwards, k(i) = k(i+3) −k(i+1), to obtain k(i) for negative whole i. Indeed, in the polynomials, we have x.x.x = 1 +x so 1 = x.x.x −x = (x +1).x.(x −1) so each of x and x ±1 has a multiplicative inverse, the product of the other two. In particular 1/x is x.x −1 and the reverse iterator is generated by powers of [−1, 1, 0] just as the forward one is generated by [0, 0, 1]. We have [-1,0,1].[0,0,1].[1,0,1] = [1,0,0], the multiplicative identity; so any product of powers of [-1,0,1], [0,0,1] and [1,0,1] is necessarily invertible.
These three generators are [a, 0, 1] for a in {-1, 0, 1}. Now, suppose we multiply one of them by [u, v, w] and all three entries in the resulting [a.u +v, a.v +w, a.w +u +v] have some common factor p. Working now modulo p, we thus have v = −a.u, w = −a.v, u +v = −a.w. In the case a = 0, this makes all three of u, v, w multiples of p. In the other cases, we have a.a = 1, so w = −a.v = u, v = −u −a.u = −(1+a).u; combining this with the previous, u = −a.v = a.(1+a).u = (a+1).u; subtracting u we get 0 = a.u, whence u = 0 = w and v = −(1+a).u = 0 so, again, all three of u, v, w have p as a factor. Consequently, any invertible triple – i.e. any product of only our three basic invertible triples – must have highest common factor 1 among its entries (since this is true of [1, 0, 0], the empty product, and multiplication by the basic factors never adds a new common factor). As our starting triple [3, 0, 2] or [0, 2, 3] has highest common factor 1 and our forward iterator is multiplication by one of the basic invertible triples, we can also assert that the highest common factor of any three consecutive values of Perrin's iterator (whether wound forwards or backwards) is 1.
Written by Eddy.