I gave up on idea "one".
Idea "two" is so-so - it works but is slow. I can provide the code if
you wish.
Idea "three" is "borrowed" (plagiarised) from the Pari code quoted at
the codegolf site.
It's obvious when you see it, but it's still slow, though not as slow
as my method two.
Here's my J version:
mission =: 3 : 0 NB. check numbers from x to y
2 mission y
:
i =. x
I =. =i.3
n =. y
l =. '' NB. required list of incorrect pseudo-primes such as 271441,
904631 ...
mat =. 0 1 0, 0 0 1,:1 1 0x
while. i < y do.
if. 0=100000| i do. NB. reassuring occasional output
smoutput i
wd'msgs'
end.
if. -. 1 p: i do. NB. no point testing true primes!
mati =. mat (i mMpower) i
trace=. +/,I*mati
if. 0 = i|trace do.
l =. l, i
end.
end.
i =. >: i
end.
l
)
NB. this is from my "finite.ijs" script of useful finite arithmetic
functions
NB. Matrix x to high power y modulo m, using binary rep of y with
repeated squaring of the matrix
mMpower =: 1 : 0
:
M =. = i.# M0 =. x
b =. #:y
NB. mm =. (m&|@+)/ . (m&|@*)
NB. mm =. (m&|@+)/ . (m&|@*) NB. seems too slow AND inaccurate!!!
mm =. <.@:(m| +/ . *)
while.#b do.
if. {:b do. M =. M0 mm M end.
M0 =. mm~ M0 NB. square it
b =. }:b
end.
M
)
timer' mission 2715' NB. best to use fixed font on these .........
+-------++
|1.96802||
+-------++
timer' mission 27150' NB. idea "two" is about 3 times slower here....
+------++
|35.374||
+------++
NB. The Pari GP version takes about 25 seconds for the first 27150 numbers.
timer'271430 mission 271500' NB. ok if you know an answer!
+--------+------+
|0.136963|271441|
+--------+------+
It would be good to be able to get directly from i | mat^ i to (i+1) |
mat ^(i+1) rather than have to construct each from scratch. I don't see how.
All rather long-winded for golfers I should think.
Mike
On 06/01/2014 10:11, Mike Day wrote:
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
