Fibonacci by reciprocal

Everything else that doesn't fall into one of the other PB categories.
jack
Addict
Addict
Posts: 1336
Joined: Fri Apr 25, 2003 11:10 pm

Fibonacci by reciprocal

Post by jack »

here's an unusual way to generate the Fibonacci numbers https://www.futilitycloset.com/2015/06/ ... o-order-4/ also https://math.stackexchange.com/question ... i-sequence, to illustrate we generate the first 5 Fibonacci numbers using double

Code: Select all

Define.d d, r
Define.l i, j
Define.s s

d=998999
r=1/d
s=StrD(r,16)
s=Right(s,Len(s)-2) ;trim 0.
OpenConsole()
j=1
For i = 1 To 5
    PrintN(Mid(s, j, 3))
    j = j + 3
Next
Input()
CloseConsole()
now we extend that idea dividing 1 by 9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999899999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999 to 48380 digits
the division algorithm is by Dr. David M. Smith http://dmsmith.lmu.build/MComp1996.pdf

Code: Select all

; using the algorithm by Dr. David M. Smith
; "A Multiple-Precision Division Algorithm"
; http://dmsmith.lmu.build/MComp1996.pdf
;
; adapted to PureBASIC by jack
;
; Fibonacci demo

#dimension = 12096 ; 48380 digits
Dim n.d(#dimension)
Dim d.d(#dimension)
Dim result.d(#dimension)
Define.d j, t
Define.s s, digit
Define.i i

Declare divide (Array result.d(1), Array n.d(1), Array d.d(1))
; n(1) holds the exponent of n, n(2) holds the first digit of the numerator
; in sci notation it;s 10^1 * .1 or .1e1
n(1) = 1: n(2) = 1
; d(1) holds the exponent of d, d(2) holds the first digit of the denominator
; d(3) an onward hold the rest of the denominator
d(1) = 202: d(2) = 9

For i = 3 To 50 + 3
    d(i) = 9999
Next
d(24 + 3) = 9998
d(50 + 3) = 9000

t = ElapsedMilliseconds()
    divide(result(), n(), d())
t = ElapsedMilliseconds() - t

s = ""
For i = 2 To #dimension
    digit = Trim(StrD(result(i)))
    While Len(digit) < 4
        digit = "0" + digit
    Wend
    s = s + digit
Next

digit=Space(-result(1)-1)
ReplaceString(digit," ","0",#PB_String_InPlace)

s=digit+s

j = 1

OpenConsole()

For i = 1 To 481
    PrintN(Mid(s, j, 100))
    j = j + 101
Next

PrintN("elapsed time for the division is "+ StrD(t/1000)+ " seconds")

Input()

CloseConsole()

Procedure min (a.l, b.l)
  If a < b
    ProcedureReturn a
  Else
    ProcedureReturn b
  EndIf
EndProcedure

Procedure RealW (Array w.d(1), j.l)
    Define.d wx
    wx = ((w(j - 1) * 10000 + w(j)) * 10000 + w(j + 1)) * 10000
    If #dimension >= (j + 2)
      wx = wx + w(j + 2)
    EndIf
    ProcedureReturn wx
EndProcedure

Procedure subtract (Array w.d(1), q.l, Array d.d(1), ka.l, kb.l)
    Define.l j
    For j = ka To kb
        w(j) = w(j) - q * d(j - ka + 2)
    Next
EndProcedure

Procedure normalize (Array w.d(1), ka.l, q.l)
    w(ka) = w(ka) + w(ka - 1) * 10000
    w(ka - 1) = q
EndProcedure

Procedure finalnorm (Array w.d(1), kb.l)
    Define.l carry, j
    For j = kb To 3 Step -1
        If w(j) < 0
            carry = ((-Int(w(j)) - 1) / 10000) + 1
        Else
            If w(j) >= 10000
                carry = -(Int(w(j)) / 10000)
            Else
                carry = 0
            EndIf
        EndIf
        w(j) = w(j) + carry * 10000
        w(j - 1) = w(j - 1) - carry
    Next
EndProcedure

Procedure divide (Array result.d(1), Array n.d(1), Array d.d(1))
    Define.l b, j, last, laststep, q, t
    Define.l stp
    Define.d xd, xn, round
    Dim w.d(#dimension + 4)
    b = 10000
    For j = #dimension To #dimension+4
        w(j) = 0
    Next
    t = #dimension - 1
    w(1) = n(1) - d(1) + 1
    w(2) = 0
    For j = 2 To #dimension
        w(j + 1) = n(j)
    Next
    xd = (d(2) * b + d(3)) * b + d(4) + d(5) / b
    laststep = t + 2
    For stp = 1 To laststep
        xn = RealW(w(), (stp + 2))
        q = Int(xn / xd)
        last = min(stp + t + 1, #dimension+4)
        subtract(w(), q, d(), stp + 2, last)
        normalize(w(), stp + 2, q)
    Next
    finalnorm(w(), laststep + 1)
    If w(2) <> 0
      laststep = laststep - 1
    EndIf
    round = w(laststep + 1) / b
    If round >= 0.5
      w(laststep) = w(laststep) + 1
    EndIf
    If w(2) = 0
        For j = 1 To t + 1
            result(j) = w(j + 1)
        Next
    Else
        For j = 1 To t + 1
            result(j) = w(j)
        Next
    EndIf
      If w(2) = 0
        result(1) = w(1) - 1
      Else
        result(1) = w(1)
      EndIf
EndProcedure