Loading ...
Sorry, an error occurred while loading the content.

Re: Lucas super-pseudoprime puzzle

Expand Messages
  • mikeoakes2
    ... That s a beautiful function tn(). I suspected Mod()s would be involved, but could never have come up with that polcoeff() construct. It can be speeded up
    Message 1 of 33 , May 24, 2010
    • 0 Attachment
      --- In primenumbers@yahoogroups.com, "djbroadhurst" <d.broadhurst@...> wrote:
      >
      > Using Pari-GP, I found a candidate solution as follows:
      >
      > tn(n,a)=2*polcoeff(lift(Mod((a+x)/Mod(2,n),x^2+4-a^2)^n),0)==Mod(a,n);
      > ts(n)=tn(n,3)&&tn(n,4)&&tn(n,5)&&tn(n,6)&&tn(n,7)&&tn(n,8);
      > {forstep(n=3,10^8,2,if(n%3&&!isprime(n)&&ts(n),
      > print(n);print(round(gettime/10^3)" seconds");break()));}
      >
      > 7056721
      > 84 seconds

      That's a beautiful function tn().
      I suspected Mod()s would be involved, but could never have come up with that polcoeff() construct.

      It can be speeded up by a factor of 6/5 by omitting the tn(n,7) test, at the expense of a very few spurions, which can easily be eliminated later by hand. Just 4 up to n=10^8, in fact:
      8646121, 17639999, 39528721, 90336961.

      Your code scales amazingly well with n: 90 secs for first 10^7, increasing by a mere 27 secs for the slab of 10^7 starting at n=10^8, on my 3.6Ghz box.

      > http://adm.lnpu.edu.ua/downloads/issues/2008/N2/adm-n2-4.pdf
      > which proves that n is a solution if and only if it is an odd
      > square-free composite integer such that for each prime p|n
      > n = +/- 1 mod p-1 ... [1]
      > n = +/- 1 mod p+1 ... [2]
      > This paper claims that 7056721 is the unique solution with n < 10^10.

      It might be interesting to program this test as it might be significantly faster for large n, not so?

      Mike
    • mikeoakes2
      ... You did very much what I did. I nearly fell off the chair when I averaged everything out and saw 0.57... :-) ... Not very. Would you buy an appeal to
      Message 33 of 33 , May 27, 2010
      • 0 Attachment
        --- In primenumbers@yahoogroups.com, "djbroadhurst" <d.broadhurst@...> wrote:
        >
        > I tried 1/n^c:
        >
        > v=[1237.1, 328.7, 105.4, 28.01, 6.22 , 1.510, 0.439, 0.0939];
        > print(vector(8,k,-log(v[k]/10^6)/(k+4)/log(10)));
        >
        > [0.5815, 0.5805, 0.5682, 0.5691, 0.5785, 0.5821, 0.5780, 0.5856]
        >
        > and then A/n^c, using the first datum to remove A:
        >
        > print(vector(7,k,-log(v[k+1]/v[1])/k/log(10)));
        >
        > [0.5756, 0.5348, 0.5484, 0.5747, 0.5827, 0.5750, 0.5885]
        >
        > In both cases c =~ Euler looks rather convincing,
        > given the statistics. Well spotted, Sir!

        You did very much what I did.
        I nearly fell off the chair when I averaged everything out and saw 0.57... :-)

        > How strongly are you committed to A = 1, for the average,
        > given the variability of the overall factor with a?

        Not very.
        Would you buy an appeal to Occam's razor, mon vieux?

        Mike
      Your message has been successfully submitted and would be delivered to recipients shortly.