Thanks for the insight. It looks like the sort of problem custom-made for Pari GP. I noticed that a fair amount of the code in the solutions to this is subsumed in libraries - many invoke an external primality checker.
On Mon, Jan 6, 2014 at 2:01 PM, Mike Day <[email protected]> wrote: > 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 > -- Devon McCormick, CFA ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
