2021.10.06 MATH 3600 Pisano

603 days ago by calkin

def F(n,m): return(mod(fibonacci(n),m)) 
       
mod(fibonacci(10),2)==1 
        
True
True
F(10,2) 
        
1
1
for i in srange(5): print(i,F(i,2)) 
        
(0, 0)
(1, 1)
(2, 1)
(3, 0)
(4, 1)
(0, 0)
(1, 1)
(2, 1)
(3, 0)
(4, 1)
for i in srange(10): print(i,F(i,3)) 
        
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 0)
(5, 2)
(6, 2)
(7, 1)
(8, 0)
(9, 1)
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 0)
(5, 2)
(6, 2)
(7, 1)
(8, 0)
(9, 1)
for i in srange(10): print(i,F(i,4)) 
        
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 3)
(5, 1)
(6, 0)
(7, 1)
(8, 1)
(9, 2)
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 3)
(5, 1)
(6, 0)
(7, 1)
(8, 1)
(9, 2)
for i in srange(25): print(i,F(i,5)) 
        
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 3)
(5, 0)
(6, 3)
(7, 3)
(8, 1)
(9, 4)
(10, 0)
(11, 4)
(12, 4)
(13, 3)
(14, 2)
(15, 0)
(16, 2)
(17, 2)
(18, 4)
(19, 1)
(20, 0)
(21, 1)
(22, 1)
(23, 2)
(24, 3)
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 3)
(5, 0)
(6, 3)
(7, 3)
(8, 1)
(9, 4)
(10, 0)
(11, 4)
(12, 4)
(13, 3)
(14, 2)
(15, 0)
(16, 2)
(17, 2)
(18, 4)
(19, 1)
(20, 0)
(21, 1)
(22, 1)
(23, 2)
(24, 3)
for i in srange(30): print(i,F(i,6)) 
        
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 3)
(5, 5)
(6, 2)
(7, 1)
(8, 3)
(9, 4)
(10, 1)
(11, 5)
(12, 0)
(13, 5)
(14, 5)
(15, 4)
(16, 3)
(17, 1)
(18, 4)
(19, 5)
(20, 3)
(21, 2)
(22, 5)
(23, 1)
(24, 0)
(25, 1)
(26, 1)
(27, 2)
(28, 3)
(29, 5)
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 3)
(5, 5)
(6, 2)
(7, 1)
(8, 3)
(9, 4)
(10, 1)
(11, 5)
(12, 0)
(13, 5)
(14, 5)
(15, 4)
(16, 3)
(17, 1)
(18, 4)
(19, 5)
(20, 3)
(21, 2)
(22, 5)
(23, 1)
(24, 0)
(25, 1)
(26, 1)
(27, 2)
(28, 3)
(29, 5)
for i in srange(30): print(i,F(i,7)) 
        
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 3)
(5, 5)
(6, 1)
(7, 6)
(8, 0)
(9, 6)
(10, 6)
(11, 5)
(12, 4)
(13, 2)
(14, 6)
(15, 1)
(16, 0)
(17, 1)
(18, 1)
(19, 2)
(20, 3)
(21, 5)
(22, 1)
(23, 6)
(24, 0)
(25, 6)
(26, 6)
(27, 5)
(28, 4)
(29, 2)
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 3)
(5, 5)
(6, 1)
(7, 6)
(8, 0)
(9, 6)
(10, 6)
(11, 5)
(12, 4)
(13, 2)
(14, 6)
(15, 1)
(16, 0)
(17, 1)
(18, 1)
(19, 2)
(20, 3)
(21, 5)
(22, 1)
(23, 6)
(24, 0)
(25, 6)
(26, 6)
(27, 5)
(28, 4)
(29, 2)
for i in srange(15): print(i,F(i,8)) 
        
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 3)
(5, 5)
(6, 0)
(7, 5)
(8, 5)
(9, 2)
(10, 7)
(11, 1)
(12, 0)
(13, 1)
(14, 1)
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 3)
(5, 5)
(6, 0)
(7, 5)
(8, 5)
(9, 2)
(10, 7)
(11, 1)
(12, 0)
(13, 1)
(14, 1)
for i in srange(80): print(i,F(i,9)) 
        
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 3)
(5, 5)
(6, 8)
(7, 4)
(8, 3)
(9, 7)
(10, 1)
(11, 8)
(12, 0)
(13, 8)
(14, 8)
(15, 7)
(16, 6)
(17, 4)
(18, 1)
(19, 5)
(20, 6)
(21, 2)
(22, 8)
(23, 1)
(24, 0)
(25, 1)
(26, 1)
(27, 2)
(28, 3)
(29, 5)
(30, 8)
(31, 4)
(32, 3)
(33, 7)
(34, 1)
(35, 8)
(36, 0)
(37, 8)
(38, 8)
(39, 7)
(40, 6)
(41, 4)
(42, 1)
(43, 5)
(44, 6)
(45, 2)
(46, 8)
(47, 1)
(48, 0)
(49, 1)
(50, 1)
(51, 2)
(52, 3)
(53, 5)
(54, 8)
(55, 4)
(56, 3)
(57, 7)
(58, 1)
(59, 8)
(60, 0)
(61, 8)
(62, 8)
(63, 7)
(64, 6)
(65, 4)
(66, 1)
(67, 5)
(68, 6)
(69, 2)
(70, 8)
(71, 1)
(72, 0)
(73, 1)
(74, 1)
(75, 2)
(76, 3)
(77, 5)
(78, 8)
(79, 4)
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 3)
(5, 5)
(6, 8)
(7, 4)
(8, 3)
(9, 7)
(10, 1)
(11, 8)
(12, 0)
(13, 8)
(14, 8)
(15, 7)
(16, 6)
(17, 4)
(18, 1)
(19, 5)
(20, 6)
(21, 2)
(22, 8)
(23, 1)
(24, 0)
(25, 1)
(26, 1)
(27, 2)
(28, 3)
(29, 5)
(30, 8)
(31, 4)
(32, 3)
(33, 7)
(34, 1)
(35, 8)
(36, 0)
(37, 8)
(38, 8)
(39, 7)
(40, 6)
(41, 4)
(42, 1)
(43, 5)
(44, 6)
(45, 2)
(46, 8)
(47, 1)
(48, 0)
(49, 1)
(50, 1)
(51, 2)
(52, 3)
(53, 5)
(54, 8)
(55, 4)
(56, 3)
(57, 7)
(58, 1)
(59, 8)
(60, 0)
(61, 8)
(62, 8)
(63, 7)
(64, 6)
(65, 4)
(66, 1)
(67, 5)
(68, 6)
(69, 2)
(70, 8)
(71, 1)
(72, 0)
(73, 1)
(74, 1)
(75, 2)
(76, 3)
(77, 5)
(78, 8)
(79, 4)
for i in srange(50): print(i,F(i,14)) 
        
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 3)
(5, 5)
(6, 8)
(7, 13)
(8, 7)
(9, 6)
(10, 13)
(11, 5)
(12, 4)
(13, 9)
(14, 13)
(15, 8)
(16, 7)
(17, 1)
(18, 8)
(19, 9)
(20, 3)
(21, 12)
(22, 1)
(23, 13)
(24, 0)
(25, 13)
(26, 13)
(27, 12)
(28, 11)
(29, 9)
(30, 6)
(31, 1)
(32, 7)
(33, 8)
(34, 1)
(35, 9)
(36, 10)
(37, 5)
(38, 1)
(39, 6)
(40, 7)
(41, 13)
(42, 6)
(43, 5)
(44, 11)
(45, 2)
(46, 13)
(47, 1)
(48, 0)
(49, 1)
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 3)
(5, 5)
(6, 8)
(7, 13)
(8, 7)
(9, 6)
(10, 13)
(11, 5)
(12, 4)
(13, 9)
(14, 13)
(15, 8)
(16, 7)
(17, 1)
(18, 8)
(19, 9)
(20, 3)
(21, 12)
(22, 1)
(23, 13)
(24, 0)
(25, 13)
(26, 13)
(27, 12)
(28, 11)
(29, 9)
(30, 6)
(31, 1)
(32, 7)
(33, 8)
(34, 1)
(35, 9)
(36, 10)
(37, 5)
(38, 1)
(39, 6)
(40, 7)
(41, 13)
(42, 6)
(43, 5)
(44, 11)
(45, 2)
(46, 13)
(47, 1)
(48, 0)
(49, 1)
for i in srange(30): print(i,F(i,16)) 
        
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 3)
(5, 5)
(6, 8)
(7, 13)
(8, 5)
(9, 2)
(10, 7)
(11, 9)
(12, 0)
(13, 9)
(14, 9)
(15, 2)
(16, 11)
(17, 13)
(18, 8)
(19, 5)
(20, 13)
(21, 2)
(22, 15)
(23, 1)
(24, 0)
(25, 1)
(26, 1)
(27, 2)
(28, 3)
(29, 5)
(0, 0)
(1, 1)
(2, 1)
(3, 2)
(4, 3)
(5, 5)
(6, 8)
(7, 13)
(8, 5)
(9, 2)
(10, 7)
(11, 9)
(12, 0)
(13, 9)
(14, 9)
(15, 2)
(16, 11)
(17, 13)
(18, 8)
(19, 5)
(20, 13)
(21, 2)
(22, 15)
(23, 1)
(24, 0)
(25, 1)
(26, 1)
(27, 2)
(28, 3)
(29, 5)
def pisano(m): for i in srange(1,10000): if F(i,m)==0 and F(i+1,m)==1: return(i) 
       
pisano(2) 
        
3
3
pisano(3) 
        
8
8
def pisano(m): i=1 while (F(i,m)==0 and F(i+1,m)==1)==False: i+=1 return(i) 
       
pisano(8) 
        
12
12
for i in srange(2,10): print(i,pisano(i)) 
        
(2, 3)
(3, 8)
(4, 6)
(5, 20)
(6, 24)
(7, 16)
(8, 12)
(9, 24)
(2, 3)
(3, 8)
(4, 6)
(5, 20)
(6, 24)
(7, 16)
(8, 12)
(9, 24)
pisano_list=[] for i in srange(2,1000): pisano_list.append([i,pisano(i)]) 
       
print(pisano_list) 
        
[[2, 3], [3, 8], [4, 6], [5, 20], [6, 24], [7, 16], [8, 12], [9, 24],
[10, 60], [11, 10], [12, 24], [13, 28], [14, 48], [15, 40], [16, 24],
[17, 36], [18, 24], [19, 18], [20, 60], [21, 16], [22, 30], [23, 48],
[24, 24], [25, 100], [26, 84], [27, 72], [28, 48], [29, 14]]
[[2, 3], [3, 8], [4, 6], [5, 20], [6, 24], [7, 16], [8, 12], [9, 24], [10, 60], [11, 10], [12, 24], [13, 28], [14, 48], [15, 40], [16, 24], [17, 36], [18, 24], [19, 18], [20, 60], [21, 16], [22, 30], [23, 48], [24, 24], [25, 100], [26, 84], [27, 72], [28, 48], [29, 14]]
list_plot(pisano_list) 
       
pisano_list=[] for i in srange(2,10000): if is_prime(i): pisano_list.append([i,pisano(i)]) 
       
list_plot(pisano_list) 
       
pisano_list_1=[] pisano_list_2=[] for i in srange(2,1000): if is_prime(i): p_value=pisano(i) if p_value>3/2*i: pisano_list_2.append([i,pisano(i)]) elif p_value>9/10*i: pisano_list_1.append([i,pisano(i)]) 
       
list_plot(pisano_list_1) 
       
list_plot(pisano_list_2) 
       
print(pisano_list_1[:20]) 
        
[[2, 3], [11, 10], [19, 18], [31, 30], [41, 40], [59, 58], [61, 60],
[71, 70], [79, 78], [109, 108], [131, 130], [149, 148], [179, 178],
[191, 190], [239, 238], [241, 240], [251, 250], [269, 268], [271, 270],
[311, 310]]
[[2, 3], [11, 10], [19, 18], [31, 30], [41, 40], [59, 58], [61, 60], [71, 70], [79, 78], [109, 108], [131, 130], [149, 148], [179, 178], [191, 190], [239, 238], [241, 240], [251, 250], [269, 268], [271, 270], [311, 310]]
print(pisano_list_2[:20]) 
        
[[3, 8], [5, 20], [7, 16], [13, 28], [17, 36], [23, 48], [37, 76], [43,
88], [53, 108], [67, 136], [73, 148], [83, 168], [97, 196], [103, 208],
[127, 256], [137, 276], [157, 316], [163, 328], [167, 336], [173, 348]]
[[3, 8], [5, 20], [7, 16], [13, 28], [17, 36], [23, 48], [37, 76], [43, 88], [53, 108], [67, 136], [73, 148], [83, 168], [97, 196], [103, 208], [127, 256], [137, 276], [157, 316], [163, 328], [167, 336], [173, 348]]
def pisano(m): i=1 while (F(i,m)==0 and F(i+1,m)==1)==False: i+=1 return(i) 
       
[F(i, 11) for i in srange(20)] 
        
[0, 1, 1, 2, 3, 5, 8, 2, 10, 1, 0, 1, 1, 2, 3, 5, 8, 2, 10, 1]
[0, 1, 1, 2, 3, 5, 8, 2, 10, 1, 0, 1, 1, 2, 3, 5, 8, 2, 10, 1]
print([F(i, 13) for i in srange(30)]) 
        
[0, 1, 1, 2, 3, 5, 8, 0, 8, 8, 3, 11, 1, 12, 0, 12, 12, 11, 10, 8, 5, 0,
5, 5, 10, 2, 12, 1, 0, 1]
[0, 1, 1, 2, 3, 5, 8, 0, 8, 8, 3, 11, 1, 12, 0, 12, 12, 11, 10, 8, 5, 0, 5, 5, 10, 2, 12, 1, 0, 1]
 
       

Let's explore when 5 is a square mod($p$).  
def is_5_square(p): for i in srange(3,(p+1)/2): if mod(i^2,p)==5: return(True) return(False) 
       
is_5_square(19) 
        
True
True
for p in srange(10,100): if is_prime(p): if is_5_square(p): print(p) 
        
11
19
29
31
41
59
61
71
79
89
11
19
29
31
41
59
61
71
79
89
for p in srange(10,100): if is_prime(p): if is_5_square(p)==False: print(p) 
        
13
17
23
37
43
47
53
67
73
83
97
13
17
23
37
43
47
53
67
73
83
97

This is consistent with the conjecture that 

5 is a square (mod $p$) if and only if $p\equiv \pm 1 $(mod 10).
def root_5(p): for i in srange(3,(p+1)/2): if mod(i^2,p)==5: return(i) 
       
for p in srange(10,100): if is_prime(p): if is_5_square(p): root5=root_5(p) print(p, root5) 
        
(11, 4)
(19, 9)
(29, 11)
(31, 6)
(41, 13)
(59, 8)
(61, 26)
(71, 17)
(79, 20)
(89, 19)
(11, 4)
(19, 9)
(29, 11)
(31, 6)
(41, 13)
(59, 8)
(61, 26)
(71, 17)
(79, 20)
(89, 19)

We wish to explore $\pi(p)$ when $p$ is $\pm 3$ (mod 10).
pisano_list_3=[] for i in srange(1000): if is_prime(10*i+3): pisano_list_3.append([10*i+3,pisano(10*i+3)]) if is_prime(10*i+7): pisano_list_3.append([10*i+7, pisano(10*i+7)]) 
       
list_plot(pisano_list_3) 
       

It appears to be the case that the majority of the values of $\pi(p)$ in this category are $2(p-1)$.  Are the remaining terms divisors of $2(p+1)$?

Creat a ratio list, containing pairs $[p, 2(p+1)/\pi(p)]$
ratio_list=[] for pisano_pair in pisano_list_3: rat= 2*(pisano_pair[0]+1)/pisano_pair[1] ratio_list.append([pisano_pair[0],rat]) 
       
list_plot(ratio_list) 
       
print(ratio_list[60:80]) 
        
[[647, 1], [653, 1], [673, 1], [677, 3], [683, 1], [727, 1], [733, 1],
[743, 3], [757, 1], [773, 1], [787, 1], [797, 7], [823, 1], [827, 1],
[853, 1], [857, 1], [863, 1], [877, 1], [883, 1], [887, 1]]
[[647, 1], [653, 1], [673, 1], [677, 3], [683, 1], [727, 1], [733, 1], [743, 3], [757, 1], [773, 1], [787, 1], [797, 7], [823, 1], [827, 1], [853, 1], [857, 1], [863, 1], [877, 1], [883, 1], [887, 1]]

How slow is our Pisano code?  It takes about $p$ computations to compute all Fibonacci numbers up to $p$ (need up to 2$p$+3 for some primes).

There are about $N/\log(N)$ primes up to $N$, so we have to do about $N^2/\log(N)$ computations of Fibonacci numbers to get all periods for primes less than $N$.  How many is this if $N=10000$?
N=10000. print(N^2/log(N)) 
        
1.08573620475813e7
1.08573620475813e7
fibonacci(10000) 
        
336447648764317832666216120051075433103021484606800639065647699746800814\
421666623681555955136337340255820653326808361593737347904838652682630408\
924630564318873545443695598274916066020998841839338646527313000888302692\
356736131351175792974378544137521305205043477016022647583189065278908551\
543661595829872796829875106312005754287834532155151038708182989697916131\
278562650331954871402142875326981879620469360978799003509623022910263681\
314931952756302278376284415403605844025721143349611800230912082870460889\
239623288354615057765832712525460935911282039252853934346209042452489294\
039017062338889910858410651831733604374707379085526317643257339937128719\
375877468974799263058370657428301616374089691784263786242128352581128205\
163702980893320999057079200643674262023897831114700540749984592503606335\
609338838319233867830561364353518921332797329081337326426526339897639227\
234078829281779535805709936910491754708089318410561463223382174656373212\
482263830921032977016480547262438423748624114530938122065649140327510866\
433945175121615265453613331113140424368548051067658434935238369596534280\
717687753283482343455573667197313927462736291082106792807847180353291311\
767789246590899386354593278945237776744061922403376386740040213303432974\
969020283281459334188268176838930720036347956231171031012919531697946076\
327375892535307725523759437884345040677155557790564504430166401194625809\
722167297586150269684431469520346149322911059706762432685159928347098912\
847067408620085871350162603120719031720860940812983215810772820763531866\
246112782455372085323653057759564300725177443150515396009051686032203491\
632226408852488524331580515348496224348482993809050704834824493274537326\
245677558790891871908036620580095947431500524025327097469953187707243768\
259074199396322659841474981936092852239450397071654431564213281576889080\
587831834049174345562705202235648464951961124602683139709750693826487066\
132645076650746115126775227486215986425307112984411826226610571635150692\
600298617049454250474913781151541399415506712562711971332527636319396069\
028956502882686083622410820505624307017949761711212330660733100599473668\
75
33644764876431783266621612005107543310302148460680063906564769974680081442166662368155595513633734025582065332680836159373734790483865268263040892463056431887354544369559827491606602099884183933864652731300088830269235673613135117579297437854413752130520504347701602264758318906527890855154366159582987279682987510631200575428783453215515103870818298969791613127856265033195487140214287532698187962046936097879900350962302291026368131493195275630227837628441540360584402572114334961180023091208287046088923962328835461505776583271252546093591128203925285393434620904245248929403901706233888991085841065183173360437470737908552631764325733993712871937587746897479926305837065742830161637408969178426378624212835258112820516370298089332099905707920064367426202389783111470054074998459250360633560933883831923386783056136435351892133279732908133732642652633989763922723407882928177953580570993691049175470808931841056146322338217465637321248226383092103297701648054726243842374862411453093812206564914032751086643394517512161526545361333111314042436854805106765843493523836959653428071768775328348234345557366719731392746273629108210679280784718035329131176778924659089938635459327894523777674406192240337638674004021330343297496902028328145933418826817683893072003634795623117103101291953169794607632737589253530772552375943788434504067715555779056450443016640119462580972216729758615026968443146952034614932291105970676243268515992834709891284706740862008587135016260312071903172086094081298321581077282076353186624611278245537208532365305775956430072517744315051539600905168603220349163222640885248852433158051534849622434848299380905070483482449327453732624567755879089187190803662058009594743150052402532709746995318770724376825907419939632265984147498193609285223945039707165443156421328157688908058783183404917434556270520223564846495196112460268313970975069382648706613264507665074611512677522748621598642530711298441182622661057163515069260029861704945425047491378115154139941550671256271197133252763631939606902895650288268608362241082050562430701794976171121233066073310059947366875

Let's create a more efficient pisano(p) function.  We'll assume that the input will always be a prime $p$.

Do we need a more efficient Fibonacci (mod p) algorithm?

We conjecture that the period is always a divisor of $g(p)$, where $g(p) =p-1$ when $p$ is $\pm 1 $(mod 10) and $g(p)=2(p+1)$ when $p$ is $\pm 3$ (mod 10).
g=10 
       
g.divisors() 
        
[1, 2, 5, 10]
[1, 2, 5, 10]
def faster_pisano(p): if p.is_prime(): if mod(p,10)==1 or mod(p,10)==9: g=p-1 elif mod(p,10)==3 or mod(p,10)==7: g=2*(p+1) d_list=g.divisors() for d in d_list: # We'll check if F_d, F_(d+1) are 0, 1 mod p if F(d,p)==0 and F(d+1,p)==1: return(d) return('Ooops') 
       
pisano(11) 
        
10
10
faster_pisano(11) 
        
10
10
pisano(13) 
        
28
28
faster_pisano(13) 
        
28
28
next_prime(10^5) 
        
100003
100003
pisano(100003) 
        
Traceback (click to the left of this block for traceback)
...
__SAGE__
^CTraceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_32.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("cGlzYW5vKDEwMDAwMyk="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpqLEXEt/___code___.py", line 3, in <module>
    exec compile(u'pisano(_sage_const_100003 )
  File "", line 1, in <module>
    
  File "/tmp/tmp2f5M0d/___code___.py", line 5, in pisano
    while (F(i,m)==_sage_const_0  and F(i+_sage_const_1 ,m)==_sage_const_1 )==False:
  File "/tmp/tmpZ0F4IW/___code___.py", line 3, in F
    return(mod(fibonacci(n),m))
  File "/usr/local/sage-6.10/local/lib/python2.7/site-packages/sage/combinat/combinat.py", line 545, in fibonacci
    return ZZ(pari(n).fibonacci())
  File "sage/libs/pari/gen.pyx", line 5326, in sage.libs.pari.gen.gen.fibonacci (/usr/local/sage-6.10/src/build/cythonized/sage/libs/pari/gen.c:123712)
  File "sage/ext/interrupt/interrupt.pyx", line 88, in sage.ext.interrupt.interrupt.sig_raise_exception (/usr/local/sage-6.10/src/build/cythonized/sage/ext/interrupt/interrupt.c:925)
KeyboardInterrupt
__SAGE__
faster_pisano(100003) 
        
200008
200008
next_prime(10^6) 
        
1000003
1000003
faster_pisano(next_prime(10^6)) 
        
2000008
2000008
faster_pisano(next_prime(10^7)) 
        
10000018
10000018
faster_pisano(next_prime(10^8)) 
        
200000016
200000016

Our faster pisano code is clearly way more efficient than the previous pisano code.  This means that we can get lots more data.
pisano_list_3=[] for i in srange(10000): if is_prime(10*i+3): pisano_list_3.append([10*i+3,faster_pisano(10*i+3)]) if is_prime(10*i+7): pisano_list_3.append([10*i+7, faster_pisano(10*i+7)]) 
       
list_plot(pisano_list_3) 
       

We can now compute out far enough to generate conjectures about how often we see each ratio $g(p)/\pi(p)$, and which ratios we see.

Suggestion: use primes up to 1000 to generate a conjecture, and primes up to 10000 to test the conjecture.

Idea: for each of the two classes of primes, compute the fraction of primes up to N with ratio 1, ratio 2, ratio 3, etc.  Do this for N=1000, and guess at a conjecture.  Compute for N=10000, and see if the results are consistent.
ratio_list=[] for pisano_pair in pisano_list_3: rat= 2*(pisano_pair[0]+1)/pisano_pair[1] ratio_list.append([pisano_pair[0],rat]) 
       
list_plot(ratio_list)