Well, consider (over breakfast)
(0 1 0, 0 0 1,:1 1 0)&(+/ . *)^:(2 3 4 5 6 - 2)3 0 2x NB. just to
get started
3 0 2
0 2 3
2 3 2
3 2 5
2 5 5
{:"1(0 1 0, 0 0 1,:1 1 0)&(+/ . *)^:( 100 101 103-2) 3 0 2x NB.
P(100 101 103)
1630580875002 2160059765855 3790640640857
103#.inv {:"1(0 1 0, 0 0 1,:1 1 0)&(+/ . *)^:( 100 101 103-2) 3 0
2x NB. P(103)=+/P100 101
1 37 67 51 50 23 59
1 83 33 87 69 90 44
3 17 101 36 17 11 0
+/2{.103#.inv {:"1(0 1 0, 0 0 1,:1 1 0)&(+/ . *)^:( 100 101 103-2)
3 0 2x
2 120 100 138 119 113 103
So my first idea is to consider that each P(i) is (ai0,... aij...) in base j
In principle it can be re-expressed in base j+2 or base j+3.
P(n) = a(n-3) (+) a(n-2) in base n
Relabel P(n-3) = a[n-3] as (b0,b1,...) in base n-3
P(n-3) = (b0-3b1+9b2-27b3..., b1-6bi2+27b3-108b4..., b2-9b3....,...)
in base n
P(n-2) = a(n-2) = (c0,c1,c2,....) in base n-2
= (c0-2c1+4c2-8c3..., c1-6c2+12c3...,c2-6c3+24c3..., ...) in
base n
(apologies for any slip-ups here!)
Second idea is to express P(n) in base M=100000007, ie _4 p: 1e8 .
It's evidently not sufficient unfortunately to just use M as a modulus.
A 33000 decimal digit number has about 4125 coefficients in base M, so
this might still be tricky.
I'll see if I can make one of these ideas work a bit later - got to go
out now.
Mike
On 06/01/2014 07:08, Devon McCormick wrote:
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
---
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
