my program is a bit ugly because I am using the decfloat routines I posted recently in posts https://www.purebasic.fr/english/viewto ... 66#p585566 and https://www.purebasic.fr/english/viewto ... 77#p585577
you need to combine the top post with the second and save as decfloat.pbi
the precision needs to be the length in decimal digits of N + the number of digits to show
Code: Select all
#number_of_digits = 640
XIncludeFile "decfloat.pbi"
OpenConsole()
CompilerIf 1
; fibonacci-number-modulo-m
Procedure fib(*result.DecFloat, *k.DecFloat, *m.DecFloat)
Define.DecFloat b, x, y, xx, temp, p2, one, tmp, tmp2, z5
Define.i i, bit_length, c
Define.s a, s
si2fp(one, 1)
str2fp(z5, ".05")
copydec(*k, b)
If fpcmp(*k, one) <= 0
copydec(*k, *result)
ProcedureReturn
EndIf
si2fp(x, 1) : si2fp(y, 0)
fplog(tmp, *k)
si2fp(tmp2, 2)
fplog(tmp2, tmp2)
fpdiv(tmp, tmp, tmp2)
fpadd(tmp, tmp, z5)
bit_length=fp2dbl(tmp)
si2fp(p2, 2)
fpipow(p2, p2, bit_length-1)
For i = bit_length-1 To 0 Step -1
;xx = x*x Mod m
fpmul(xx, x, x) : fpdiv(xx, xx, *m)
fpfrac(xx, xx)
fpmul(xx, xx, *m)
fpadd(xx, xx, z5)
fptrunc(xx, xx)
fpmul(x, x, y)
fpmul_si(x, x, 2)
fpadd(x, x, xx)
fpdiv(x, x, *m)
fpfrac(x, x)
fpmul(x, x, *m)
fpadd(x, x, z5)
fptrunc(x, x)
;y = (xx + y*y) Mod m
fpmul(y, y, y)
fpadd(y, y, xx)
fpdiv(y, y, *m)
fpfrac(y, y)
fpmul(y, y, *m)
fpadd(y, y, z5)
fptrunc(y, y)
fpdiv(temp, b, p2)
;If Bit(k,i) Then
If fptrunc_is_odd(temp)
temp = x
;x = (x + y) Mod m
fpadd(x, x, y)
fpdiv(x, x, *m)
fpfrac(x, x)
fpmul(x, x, *m)
fpadd(x, x, z5)
fptrunc(x, x)
y = temp
EndIf
fpdiv_si(p2, p2, 2)
fptrunc(p2, p2)
Next
fpadd(temp, x, z5)
copydec(temp, *result)
EndProcedure
Define.DecFloat last9, fn, m
Define.d t
Define.l i, dp, pow_of_two, show_number_of_digits
Define.s a, s, sn, sl
Define.DecFloat x, y, z, tmp, lg10
show_number_of_digits=20
pow_of_two=2048
si2fp(fn, 2)
fpipow(fn, fn, pow_of_two)
;si2fp(fn, 4784969) ; str2fp(fn, "18446744073709551616")
si2fp(m, 10)
fpipow(m, m, show_number_of_digits)
t=ElapsedMilliseconds()
si2fp(x,5)
fpSqr(x, x)
si2fp(y, 1)
fpadd(y, y, x)
fpdiv_si(y, y, 2)
fplog(y, y)
fpmul(y, y, fn)
fplog(tmp, x)
fpsub(y, y, tmp)
si2fp(tmp, 10)
fplog(lg10, tmp)
fpdiv(y, y, lg10)
;y=(Log((1+x)/2)*fn-Log(x))/Log(DecFloat(10))
fpfrac(z, y)
fpmul(y, lg10, z)
fpexp(y, y)
;y=10^z
s=Trim(fp2str(y, show_number_of_digits+2))
dp=FindString(s, ".")
a=Left(s, dp-1) +Mid(s, dp+1)
a=Left(a, show_number_of_digits)
dp=FindString(a,"e")
If dp>0
s=Left(a,dp-1)
Else
s=a
EndIf
fib(last9, fn, m)
sn=fp2str(fn)
dp=FindString(sn, ".")
sn=Trim(Left(sn, dp-1))
sl=fp2str(last9, show_number_of_digits+2)
sl=Trim(sl)
dp=FindString(sl, ".")
a=Left(sl, dp-1)
dp=FindString(a,"e")
If dp
a+Mid(sl, dp+1)
sl=Left(a,dp-1)
dp=Val(Mid(a,dp+1))
sl=Left(sl,dp+1)
Else
sl=a
EndIf
If Len(sl)<show_number_of_digits
For i=0 To show_number_of_digits-Len(sl)
sl="0"+sl
Next
EndIf
PrintN("N = "+sn)
PrintN("the first "+Str(show_number_of_digits)+" digits of Fibonacci N are "+s)
Print("the last "+Str(show_number_of_digits))
PrintN(" digits of Fibonacci N are "+sl)
t=ElapsedMilliseconds()-t
PrintN("time taken is "+StrD(t)+" milliseconds")
CompilerEndIf
Print("press return to end") : Input()
Input()
CloseConsole()
Code: Select all
N = 32317006071311007300714876688669951960444102669715484032130345427524655138867890893197201411522913463688717960921898019494119559150490921095088152386448283120630877367300996091750197750389652106796057638384067568276792218642619756161838094338476170470581645852036305042887575891541065808607552399123930385521914333389668342420684974786564569494856176035326322058077805659331026192708460314150258592864177116725943603718461857357598351152301645904403697613233287231227125684710820209725157101726931323469678542580656697935045997268352998638215525166389437335543602135433229604645318478604952148193555853611059596230656
the first 20 digits of Fibonacci N are 25061437365882579483
the last 20 digits of Fibonacci N are 61690195364736004667
time taken is 368 milliseconds