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()
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