was searching the web for Fibonacci modulo and came across https://www.uni-math.gwdg.de/tschinkel/gauss/Fibon.pdf
in page 14 it list a program in Aribas https://www.mathematik.uni-muenchen.de/~...ribas.html
the pdf is not there anymore, but I got the Aribas program here
but I found this https://math.uni-paderborn.de/fileadmin-...ibams3.pdf the Aribas program is on page 16
so then I tried fib(2**128, 100000000000000000000). and instantly it printed 17229_32409_58826_54267
so the first 20 digits of F(2**128) = 1.52627_28879_74047_15656 and last 20 digits = 17229_32409_58826_54267
the first 20+ digits were calculated to 1024 bit precision by 10**frac(((log((1+sqrt(5))/2)*n)-log(sqrt(5)))/log(10)). where n = 2**128
well, I kludged/translated the program to QB64pe using my decfloat routines
here's the relevant code, you will have to stich/replace/paste this code with the code posted above
output
in page 14 it list a program in Aribas https://www.mathematik.uni-muenchen.de/~...ribas.html
the pdf is not there anymore, but I got the Aribas program here
but I found this https://math.uni-paderborn.de/fileadmin-...ibams3.pdf the Aribas program is on page 16
Code: (Select All)
(*------------------------------------------------------------------*)
(*
** Schnelle Berechnung der Fibonacci-Zahlen mittels der Formeln
** fib(2*k-1) = fib(k)**2 + fib(k-1)**2
** fib(2*k) = fib(k)**2 + 2*fib(k)*fib(k-1)
**
** Dabei werden alle Berechnungen mod m durchgef¨uhrt
*)
function fib(k,m : integer): integer;
var
b, x, y, xx, temp: integer;
begin
if k <= 1 then return k end;
x := 1; y := 0;
for b := bit_length(k)-2 to 0 by -1 do
xx := x*x mod m;
x := (xx + 2*x*y) mod m;
y := (xx + y*y) mod m;
if bit_test(k,b) then
temp := x;
x := (x + y) mod m;
y := temp;
end;
end;
return x;
end.
(*------------------------------------------------------------------*)
set_floatprec(1024);
n:= 2**128;
fib(n, 10**40).
f:=10**frac(((log((1+sqrt(5))/2)*n)-log(sqrt(5)))/log(10));
f:=trunc(f*10**40).
so then I tried fib(2**128, 100000000000000000000). and instantly it printed 17229_32409_58826_54267
so the first 20 digits of F(2**128) = 1.52627_28879_74047_15656 and last 20 digits = 17229_32409_58826_54267
the first 20+ digits were calculated to 1024 bit precision by 10**frac(((log((1+sqrt(5))/2)*n)-log(sqrt(5)))/log(10)). where n = 2**128
well, I kludged/translated the program to QB64pe using my decfloat routines
here's the relevant code, you will have to stich/replace/paste this code with the code posted above
Code: (Select All)
_Title "decfloat-offset-v2-fib by Jack"
$Console:Only
_Dest _Console
Option _Explicit
Declare CustomType Library
Sub memcpy Alias "memmove" (ByVal dest As _Offset, ByVal source As _Offset, ByVal bytes As Long)
End Declare
Const NUMBER_OF_DIGITS = 100
Const NUM_DIGITS = 8 * (1 + NUMBER_OF_DIGITS \ 8)
Const NUM_DWORDS = NUM_DIGITS \ 8
Const BIAS = 1073741824 '2 ^ 30
Const MANTISSA_BYTES = (NUM_DWORDS + 1) * 4
Type decfloat
As Long sign
As _Unsigned Long exponent
As String * Mantissa_bytes mantissa
End Type
' Error definitions
Const DIVZ_ERR = 1 'Divide by zero
Const EXPO_ERR = 2 'Exponent overflow error
Const EXPU_ERR = 3 'Exponent underflow error
Dim Shared As decfloat pi_dec, tan_half_num(15), tan_half_den(14)
Dim As decfloat x, y, z, last9, f, m, lt, lx, one
Dim As Double t
Dim As Long dp, pow_of_two, show_number_of_digits
Dim As String s, sn, sl
'initialize_fp
show_number_of_digits = 20
pow_of_two = 1024
'si2fp f, 2, NUMBER_OF_DIGITS
'fpipow f, f, pow_of_two, NUMBER_OF_DIGITS
si2fp f, 1000000, NUMBER_OF_DIGITS
'fn=4784969 '1000 '"10000000" '"18446744073709551616" '
si2fp m, 10, NUMBER_OF_DIGITS
fpipow m, m, show_number_of_digits, NUMBER_OF_DIGITS
t = Timer(.0001)
si2fp one, 1, NUMBER_OF_DIGITS
si2fp x, 10, NUMBER_OF_DIGITS
fplog lt, x, NUMBER_OF_DIGITS
si2fp x, 5, NUMBER_OF_DIGITS
fpsqr x, x, NUMBER_OF_DIGITS
fplog lx, x, NUMBER_OF_DIGITS
fpadd y, x, one, NUMBER_OF_DIGITS
fpdiv_si y, y, 2, NUMBER_OF_DIGITS
fplog y, y, NUMBER_OF_DIGITS
fpmul y, y, f, NUMBER_OF_DIGITS
fpsub y, y, lx, NUMBER_OF_DIGITS
fpdiv y, y, lt, NUMBER_OF_DIGITS
fpfrac z, y, NUMBER_OF_DIGITS
si2fp x, 10, NUMBER_OF_DIGITS
fpmul y, lt, z, NUMBER_OF_DIGITS
fpexp y, y, NUMBER_OF_DIGITS
s = _Trim$(fp2str_fix$(y, show_number_of_digits))
dp = InStr(s, ".")
s = Left$(s, dp - 1) + Mid$(s, dp + 1)
s = Left$(s, 20)
fib last9, f, m
sn = fp2str_fix(f, 310)
dp = InStr(sn, ".")
sn = _Trim$(Left$(sn, dp - 1))
sl = fp2str_fix(last9, show_number_of_digits + 2)
dp = InStr(sl, ".")
sl = _Trim$(Left$(sl, dp - 1))
Print "the first"; Len(s); " digits of Fibonacci "; sn; " are "; s
Print "the last "; show_number_of_digits;
Print " digits of Fibonacci "; sn; " are "; sl
t = Timer(.0001) - t
Print "time taken is"; t; " seconds"
Sub fib (result As decfloat, k As decfloat, m As decfloat)
Dim As decfloat b, x, y, xx, temp, p2, rn, one, two, l2
Dim As Integer i, bit_length
si2fp one, 1, NUMBER_OF_DIGITS
b = k
If fpcmp(k, one, NUMBER_OF_DIGITS) <= 0 Then
result = k
Exit Sub
End If
si2fp l2, 2, NUMBER_OF_DIGITS
fplog l2, l2, NUMBER_OF_DIGITS
str2fp x, ".05"
fplog y, k, NUMBER_OF_DIGITS
fpdiv y, y, l2, NUMBER_OF_DIGITS
fpadd y, y, x, NUMBER_OF_DIGITS
fpfix y, y
bit_length = fp2dbl(y)
si2fp y, 0, NUMBER_OF_DIGITS
x = one
si2fp p2, 2, NUMBER_OF_DIGITS
si2fp two, 2, NUMBER_OF_DIGITS
fpipow p2, p2, bit_length - 1, NUMBER_OF_DIGITS
str2fp rn, ".5"
For i = bit_length - 1 To 0 Step -1
'xx = x*x Mod m
fpmul xx, x, x, NUMBER_OF_DIGITS
fpdiv xx, xx, m, NUMBER_OF_DIGITS
fpfrac xx, xx, NUMBER_OF_DIGITS
fpmul xx, xx, m, NUMBER_OF_DIGITS
fpadd xx, xx, rn, NUMBER_OF_DIGITS
fpfix xx, xx
fpmul_si x, x, 2, NUMBER_OF_DIGITS
fpmul x, x, y, NUMBER_OF_DIGITS
fpadd x, x, xx, NUMBER_OF_DIGITS
fpdiv x, x, m, NUMBER_OF_DIGITS
fpfrac x, x, NUMBER_OF_DIGITS
fpmul x, x, m, NUMBER_OF_DIGITS
fpadd x, x, rn, NUMBER_OF_DIGITS
fpfix x, x
'y = (xx + y*y) Mod m
fpmul y, y, y, NUMBER_OF_DIGITS
fpadd y, y, xx, NUMBER_OF_DIGITS
fpdiv y, y, m, NUMBER_OF_DIGITS
fpfrac y, y, NUMBER_OF_DIGITS
fpmul y, y, m, NUMBER_OF_DIGITS
fpadd y, y, rn, NUMBER_OF_DIGITS
fpfix y, y
fpdiv temp, b, p2, NUMBER_OF_DIGITS
'If Bit(k,i) Then
If fpfix_is_odd(temp) Then
temp = x
'x = (x + y) Mod m
fpadd x, x, y, NUMBER_OF_DIGITS
fpdiv x, x, m, NUMBER_OF_DIGITS
fpfrac x, x, NUMBER_OF_DIGITS
fpmul x, x, m, NUMBER_OF_DIGITS
fpadd x, x, rn, NUMBER_OF_DIGITS
fpfix x, x
y = temp
End If
fpdiv_si p2, p2, 2, NUMBER_OF_DIGITS
fpfix p2, p2
Next
result = x
End Sub
output
Quote:the first 20 digits of Fibonacci 1000000 are 19532821287077577316
the last 20 digits of Fibonacci 1000000 are 68996526838242546875
time taken is 4.000000000814907D-03 seconds