I don't think you need worry about testing primality of 33000 digit
numbers, only numbers less than 100 million.
The overflow problem concerns the size not of N but of P(N) for N <=
1e8; for what non-prime values of N does 0 = N|P(N) ?
I wonder if there's some finite arithmetic trick to keeping these
numbers small. It would be easy for a fixed modulus, but here
the modulus increases by one at each step!?.... Time for bed over here
(UK). Perhaps I'll dream the answer.
Mike
On 05/01/2014 21:19, Devon McCormick wrote:
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...
---
This email is free from viruses and malware because avast! Antivirus protection
is active.
http://www.avast.com
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
