I was looking at this "code golf" problem - http://codegolf.stackexchange.com/questions/16843/find-all-perrin-pseudoprimes-less-than-100-million- and considering a J solution. Unfortunately, this Fibonacci-like sequence grows quite quickly and my J implementation using extended integers runs for about 5 hours to find the first 271,445 members of the sequence (a 12 GB vector). Also, J's built-in primality tester has difficulty with 33,000 digit numbers.
In any case, given the obvious parallel of Perrin with Fibonacci, I looked into deriving a closed form solution like Binet's formula for Perrin. Following the logic here - http://gozips.uakron.edu/~crm23/fibonacci/fibonacci.htm - I come up with a characteristic equation with a term of (r^3)-r-1. I can solve for the zeros simply enough - p. _1 _1 0 1 +-+--------------------------------------------+ |1|1.32472 _0.662359j0.56228 _0.662359j_0.56228| +-+--------------------------------------------+ but this has limited precision. I can try my old buddy Newton, Newton=: 1 : '] - u % u d. 1' but run into problems (_1 _1 0 1&(] p.~ [: p. [)) Newton 1 |domain error: Newton | ]-u% u d.1 |Newton[0] Presuming a difficulty of the symbolic differentiation, I tried this variant: NewtonEmp=: 1 : '] - u % u D. 1' (_1 _1 0 1&(] p.~ [: p. [)) NewtonEmp 1 1.5j_1.38778e_10 (_1 _1 0 1&(] p.~ [: p. [)) NewtonEmp^:(5+i.5) ] 1 1.32472j_3.95572e_22 1.32472j_4.40342e_29 1.32472j_5.44617e_36 1.32472j_6.73583e_43 1.32472j_8.33089e_50 which gives me answers but fails to respond to extended precision arguments as I'd like: 0j25":,.(_1 _1 0 1x&(] p.~ [: p. [)) NewtonEmp^:(5+i.5) ] 1x 1.3247179572448169000000000 1.3247179572447461000000000 1.3247179572447461000000000 1.3247179572447461000000000 1.3247179572447461000000000 I also tried a variant of Newton where I supply the first differential explicitly: Newton2=: 2 : ']-(u"0)%(v"0)' But this still fails to provide me with further precision: 0j25":,.(_1 _1 0 1x&(] p.~ [: p. [)) Newton2 ((1x - ~ 3x * *:))^:(5+i.5)]1x 1.3247179572447898000000000 1.3247179572447461000000000 1.3247179572447461000000000 1.3247179572447461000000000 1.3247179572447461000000000 Figuring the "p." is what's coercing my extended precision to floating point, I tried a tacit version of "p. _1 _1 0 1" - 0j50":,.(((3x ^~ ]) - 1x + ]) Newton2 (1x -~ 3 * 2x ^~ ]))^:(i.8) ] 13r10 1.30000000000000000000000000000000000000000000000000 1.32530712530712530712530712530712530712530712530713 1.32471828046117299147326745427717351001937908806959 1.32471795724484337903599716224738987122840620163678 1.32471795724474602596090886331016451129247469958035 1.32471795724474602596090885447809734073440405690173 1.32471795724474602596090885447809734073440405690173 1.32471795724474602596090885447809734073440405690173 This looks promising... -- Devon McCormick, CFA ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
