3 and 5 are twin primes, and the 3rd and 5th Fibonacci numbers are both prime. 5 and 7 are twin primes, and the 5th and 7th Fibonacci numbers are both prime. 11 and 13 are twin primes, and the 11th and 13th Fibonacci numbers are both prime. Are there any other twin prime pairs (p, p+2) such that the pth and (p+2)th Fibonacci numbers are both prime?

(In reply to re: Exploration by Daniel)

The new findings from Mathematica point out a flaw in my program. It relied on the prmdiv function of UBASIC to return the same number as its argument for prime values. However, beyond somewhere near the 72nd Fibonacci number, 498454011879264, UBASIC returns a zero as the prmdiv, instead of crashing, misleadingly making the number seem to be not prime, as diagnosed by:

list
   10   OldFib=1:Fib=1:NewFib=2:FibNo=2
   20   while Ct<47
   30     NxtFib=Fib+NewFib
   40     OldFib=Fib:Fib=NewFib:NewFib=NxtFib
   50     inc FibNo
   55   if prmdiv(NewFib)=0 then print "stopping at F(";FibNo;")",Fib:end
   60     if prmdiv(OldFib)=OldFib and prmdiv(NewFib)=NewFib then
   70        :print FibNo-1;OldFib,FibNo+1;NewFib,:inc Ct
   80        :if prmdiv(FibNo-1)=FibNo-1 and prmdiv(FibNo+1)=FibNo+1 then
   90          :print "*":else print
  100   wend
OK
run
 2  1    4  3
 3  2    5  5   *
 5  5    7  13  *
 11  89          13  233        *
stopping at F( 72 )      498454011879264

A modified program, using a probabilitstic prime test, produces the same results as the Mathematica findings:

   10   OldFib=1:Fib=1:NewFib=2:FibNo=2
   20   while Ct<47
   30     NxtFib=Fib+NewFib
   40     OldFib=Fib:Fib=NewFib:NewFib=NxtFib
   50     inc FibNo
   60     if fnPrime(OldFib)=1 and fnPrime(NewFib)=1 then
   70        :print FibNo-1;OldFib,FibNo+1;NewFib,:inc Ct
   80        :if fnPrime(FibNo-1)=1 and fnPrime(FibNo+1)=1 then
   90          :print "*":else print
  100   wend
  980   end
  990   '
10000   fnOddfact(N)
10010   local K=0,P
10030   while N @ 2=0
10040     N=N
10050     K=K+1
10060   wend
10070   P=pack(N,K)
10080   return(P)
10090   '
10100   fnPrime(N)
10110   local I,X,J,Y,Q,K,T,Ans
10115   if N=2 then Ans=1:goto *EndPrime
10120   if N=1 or N @ 2=0 then Ans=0:goto *EndPrime
10125   O=fnOddfact(N-1)
10130   Q=member(O,1)
10140   K=member(O,2)
10150   I=0
10160   repeat
10170     repeat
10180       X=fnLrand(N)
10190     until X>1
10200     J=0
10210     Y=modpow(X,Q,N)
10220     loop
10230       if or{and{J=0,Y=1},Y=N-1} then goto *ProbPrime
10240       if and{J>0,Y=1} then goto *NotPrime
10250       J=J+1
10260       if J=K then goto *NotPrime
10270       Y=(Y*Y)@N
10280     endloop
10290    *ProbPrime
10300     I=I+1
10310   until I>50
10320   Ans=1
10330   goto *EndPrime
10340   *NotPrime
10350   Ans=0
10360   *EndPrime
10370   return(Ans)
10380   '
10400   fnLrand(N)
10410   local R
10415   N=int(N)
10420   R=(int(rnd*10^(alen(N)+2)))@N
10430   return(R)
10440   '
10500   fnNxprime(X)
10510   if N @ 2=0 then X=X+1
10520   while fnPrime(X)=0
10530     X=X+2
10540   wend
10550   return(X)
10560   '

 3  2    5  5   *
 5  5    7  13  *
 11  89          13  233        *
 431  529892711006095621792039556787784670197112759029534506620905162834769955134424689676262369         
433  1387277127804783827114186103186246392258450358171783690079918032136025225954602593712568353       *

 569  36684474316080978061473613646275630451100586901195229815270242868417768061
193560857904335017879540515228143777781065869    
571  96041200618922553823942883360924865026104917411877067816822264789029014378308478864192589084185254331637646183008074629   *

Edited on May 4, 2013, 11:21 am

Posted by Charlie on 2013-05-04 11:20:35