You make a good point that the modulus reduction brings everything into a
mundane range but my brute force rendition of the problem starts like this:
perrinNext=: ],[:+/_2 _3{]
perrinSeq=: 3 : 'perrinNext^:y]3 0 2'
6!:2 'pp=. perrinSeq 271442'
17721.4
This should get me to the first violation where 0=N|N{pp is not prime but
the elements of "pp" are getting quite large:
#":_1{pp
33150
The prime "N" to index into "pp" and take the modulus:
$nn=. p:i.p:^:_1]#pp
23760
_3{.nn
271393 271409 271429
This brute force check takes a very long time:
6!:2 'pres=. nn|nn{pp'
32251.6
$pres
23760
pres +/ . = 0
23760
and I still don't seem to have found this first exception.
I think you're implying that I need to be more clever about this and avoid
calculating P[N] but directly calculate N|P[N] for each new N? I'll sleep
on it now, too.
On Sun, Jan 5, 2014 at 7:37 PM, Mike Day <[email protected]> wrote:
> 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
>
--
Devon McCormick, CFA
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
