+- QB64 Phoenix Edition (https://qb64phoenix.com/forum)
+-- Forum: QB64 Rising (https://qb64phoenix.com/forum/forumdisplay.php?fid=1)
+--- Forum: Code and Stuff (https://qb64phoenix.com/forum/forumdisplay.php?fid=3)
+---- Forum: Programs (https://qb64phoenix.com/forum/forumdisplay.php?fid=7)
+---- Thread: decfloat -- again (/showthread.php?tid=886)
decfloat -- again - Jack - 09-13-2022
Code: (Select All)
_Title "defloat by Jack"
$NoPrefix
$Console:Only
Dest Console
Option Explicit
Const NUMBER_OF_DIGITS = 2000
Const NUM_DIGITS = NUMBER_OF_DIGITS
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)
Call initialize_fp
Dim As decfloat n, m, z1, z2, z3, z4
Dim As String s
Dim As Double t
Dim As Long prec, r
t = Timer(.0001)
n = pi_dec
str2fp m, "1e2000"
fpmul n, n, m, NUMBER_OF_DIGITS
fpnroot m, n, 3, NUMBER_OF_DIGITS
t = Timer(.0001) - t
Print "(Pi * 10^2000) ^ (1/3)"
Print fp2str(m, NUMBER_OF_DIGITS)
Print "elapsed time "; t; " seconds"
Print
'================================================
Input "enter the number digits for output (max 2000) "; prec
If prec = 0 Then prec = 2000
If prec > 2000 Then prec = 2000
If prec < 0 Then prec = 80
s = " "
Print "Please enter a number for root extraction, an empty input will exit"
While s <> ""
Input "enter number "; s
If s = "" Then Exit While
If UCase$(s) = "PI" Then
n = pi_dec
Else
str2fp n, s
End If
If n.sign <> 0 Then
Print "can't handle roots of negative numbers, will make the number positive"
n.sign = 0
End If
Input "enter an Integer root# "; r
If r = 0 Then
Print 1
t = 0
ElseIf Abs(r) > 9999999 Then
Print "Absolute value for root# too large"
Print 1
t = 0
Else
t = Timer(.0001)
fpnroot m, n, r, prec + 4
t = Timer(.0001) - t
Print fp2str(m, prec)
End If
Print "elapsed time "; t; " seconds"
Wend
'================================================
'Call str2fp(n, "2")
't = Timer
'Call mp_tan(z, n, NUMBER_OF_DIGITS)
't = Timer - t
'Print fp2str(z, NUMBER_OF_DIGITS)
'Print t
't = Timer
'Call fptan(z, n, NUMBER_OF_DIGITS)
't = Timer - t
'Print fp2str(z, NUMBER_OF_DIGITS)
'Print t
'Print fpfix_is_odd(n)
't = Timer
'Call fpfactorial(z, 10000, NUMBER_OF_DIGITS)
't = Timer - t
'Print "factorial(10000) to "; NUMBER_OF_DIGITS; " digits"
'Print fp2str(z, NUMBER_OF_DIGITS)
'Print
'Print "in "; t; " seconds"
'Print
'Call str2fp(n, "3")
'Call str2fp(m, "9")
't = Timer
'Call fpdiv(z, n, m, NUMBER_OF_DIGITS)
't = Timer - t
'Print "3/9 to "; NUMBER_OF_DIGITS; " digits"
'Print fp2str(z, NUMBER_OF_DIGITS)
'Print
'Print "in "; t; " seconds"
'Print
't = Timer
'Call fpmul(z, z, z, NUMBER_OF_DIGITS)
't = Timer - t
'Print "and square it"
'Print fp2str(z, NUMBER_OF_DIGITS)
'Print
'Call str2fp(n, "2")
't = Timer
'Call fplog(z, n, NUMBER_OF_DIGITS)
't = Timer - t
'Print "ln(2) to "; NUMBER_OF_DIGITS; " digits"
'Print fp2str(z, NUMBER_OF_DIGITS)
'Print
'Print "in "; t; " seconds"
'Print
't = Timer
'Call fpexp(z, z, NUMBER_OF_DIGITS)
't = Timer - t
'Print "and the inverse"
'Print fp2str(z, NUMBER_OF_DIGITS)
'Print
'Print "in "; t; " seconds"
'Call si2fp(n, 3, NUMBER_OF_DIGITS)
'Print
'Print "reciprocal of 3"
't = Timer
'Call fpinv(z, n, NUMBER_OF_DIGITS)
't = Timer - t
'Print fp2str(z, NUMBER_OF_DIGITS)
'Print
'Print "in "; t; " seconds"
Sub fptan_half (result As decfloat, x As decfloat, digits_in As Unsigned Long)
Dim As Long digits
digits = digits_in
If digits > NUMBER_OF_DIGITS Then digits = NUMBER_OF_DIGITS
Dim As decfloat accum, xx, thn, thd
Call fpmul(xx, x, x, digits)
'the folowing is a rational polynomial for tan(x/2) derived from the continued fraction
' x^2
'tan(x/2) = 2 - -------------------------
' x^2
' 6 - ----------------------
' x^2
' 10 - ------------------
' x^2
' 14 - --------------
' x^2
' 18 - ---------
' x^2
' 22 - -----
' 26 - ...
' and so on
'
'the nice quality of this expansion is that you can calculate sin, cos and tan very easily
'
'if y = tan(x/2) then
'
' 2*x*y
'sin(x) = -------
' y^2+x^2
'
' y^2-x^2
'cos(x) = ---------
' y^2+x^2
'
' 2*x*y
'tan(x) = -------
' y^2-x^2
Call fpmul(thn, xx, tan_half_num(0), digits)
Call fpsub(thn, thn, tan_half_num(1), digits)
Call fpmul(thn, thn, xx, digits): Call fpadd(thn, tan_half_num(2), thn, digits)
Call fpmul(thn, thn, xx, digits): Call fpadd(thn, tan_half_num(3), thn, digits)
Call fpmul(thn, thn, xx, digits): Call fpadd(thn, tan_half_num(4), thn, digits)
Call fpmul(thn, thn, xx, digits): Call fpadd(thn, tan_half_num(5), thn, digits)
Call fpmul(thn, thn, xx, digits): Call fpadd(thn, tan_half_num(6), thn, digits)
Call fpmul(thn, thn, xx, digits): Call fpadd(thn, tan_half_num(7), thn, digits)
Call fpmul(thn, thn, xx, digits): Call fpadd(thn, tan_half_num(8), thn, digits)
Call fpmul(thn, thn, xx, digits): Call fpadd(thn, tan_half_num(9), thn, digits)
Call fpmul(thn, thn, xx, digits): Call fpadd(thn, tan_half_num(10), thn, digits)
Call fpmul(thn, thn, xx, digits): Call fpadd(thn, tan_half_num(11), thn, digits)
Call fpmul(thn, thn, xx, digits): Call fpadd(thn, tan_half_num(12), thn, digits)
Call fpmul(thn, thn, xx, digits): Call fpadd(thn, tan_half_num(13), thn, digits)
Call fpmul(thn, thn, xx, digits): Call fpadd(thn, tan_half_num(14), thn, digits)
Call fpmul(thn, thn, xx, digits): Call fpadd(thn, tan_half_num(15), thn, digits)
Call fpsub(thd, xx, tan_half_den(0), digits)
Call fpmul(thd, thd, xx, digits): Call fpadd(thd, tan_half_den(1), thd, digits)
Call fpmul(thd, thd, xx, digits): Call fpadd(thd, tan_half_den(2), thd, digits)
Call fpmul(thd, thd, xx, digits): Call fpadd(thd, tan_half_den(3), thd, digits)
Call fpmul(thd, thd, xx, digits): Call fpadd(thd, tan_half_den(4), thd, digits)
Call fpmul(thd, thd, xx, digits): Call fpadd(thd, tan_half_den(5), thd, digits)
Call fpmul(thd, thd, xx, digits): Call fpadd(thd, tan_half_den(6), thd, digits)
Call fpmul(thd, thd, xx, digits): Call fpadd(thd, tan_half_den(7), thd, digits)
Call fpmul(thd, thd, xx, digits): Call fpadd(thd, tan_half_den(8), thd, digits)
Call fpmul(thd, thd, xx, digits): Call fpadd(thd, tan_half_den(9), thd, digits)
Call fpmul(thd, thd, xx, digits): Call fpadd(thd, tan_half_den(10), thd, digits)
Call fpmul(thd, thd, xx, digits): Call fpadd(thd, tan_half_den(11), thd, digits)
Call fpmul(thd, thd, xx, digits): Call fpadd(thd, tan_half_den(12), thd, digits)
Call fpmul(thd, thd, xx, digits): Call fpadd(thd, tan_half_den(13), thd, digits)
Call fpmul(thd, thd, xx, digits): Call fpadd(thd, tan_half_den(14), thd, digits)
Call fpdiv(accum, thn, thd, digits) 'tan(x/2)
result = accum
End Sub
Sub mp_tan (result As decfloat, x As decfloat, digits_in As Unsigned Long)
Dim As Long digits
digits = digits_in
If digits > NUMBER_OF_DIGITS Then digits = NUMBER_OF_DIGITS
Dim As decfloat ab, accum, factor, x_2, xx, pi1, pi2, circ, tmp, tmp2, tmp3
Dim As decfloat half
Dim As Long limit2, i, sign
pi1 = pi_dec
Call fpdiv_si(pi2, pi_dec, 2, digits)
x_2 = x
For i = 1 To NUM_DWORDS
Call set(half.mantissa, i, 0)
Next
Call set(half.mantissa, 0, 50000000)
half.exponent = BIAS
half.sign = 0
'calculate the sign of tan(x)
'if integer part of x/(pi/2) is odd then the sign is negative
'unless x is negative then the sign is positive
Call fpdiv(tmp, x, pi2, digits)
sign = fpfix_is_odd(tmp)
If sign = 1 Then
sign = &H8000
Else
sign = 0
End If
Call fpmul_si(circ, pi1, 2, digits)
Call fpabs(ab, x)
If fpcmp(ab, circ, digits) = 1 Then
'======== centralize ==============
'floor/ceil to centralize
Call fpadd(pi1, pi1, pi1, digits) 'got 2*pi
Call fpdiv(tmp, x_2, pi1, digits)
tmp2 = tmp
Call fpfix(tmp, tmp) 'int part
Call fpsub(tmp, tmp2, tmp, digits) 'frac part
Call fpmul(tmp, tmp, pi1, digits)
x_2 = tmp
End If
'==================================
'lm = digits of precision, here's a polynomial fit to get an exstimate for limit2
'works well for precision from 60 to 10000 digits
limit2 = 1 + Int(-0.45344993886092585968# + 0.022333002852398072433# * NUM_DIGITS + 5.0461814408333079844D-7 * NUM_DIGITS * NUM_DIGITS - 4.2338453039804235772D-11 * NUM_DIGITS * NUM_DIGITS * NUM_DIGITS)
Call si2fp(factor, 5, digits)
For i = 1 To limit2
Call fpmul_si(factor, factor, 5, digits)
Next
'factor=factor^limit2
Call fpdiv(x_2, x_2, factor, digits) 'x_2=x_2/5^limit2
Call fptan_half(accum, x_2, digits)
Call fpmul(xx, x_2, x_2, digits)
'now convert to sin(x)
Call fpmul_si(tmp, x_2, 2, digits)
Call fpmul(tmp, tmp, accum, digits)
Call fpmul(tmp2, accum, accum, digits)
Call fpadd(tmp2, tmp2, xx, digits)
Call fpdiv(accum, tmp, tmp2, digits)
'multiply the result by 5^limit2
For i = 1 To limit2 + 1
Call fpmul(tmp, accum, accum, digits)
Call fpmul(tmp2, accum, tmp, digits)
'sin(5*x) = 5 * sin(x) - 20 * sin(x)^3 + 16 * sin(x)^5
Call fpmul_si(accum, accum, 5, digits)
Call fpmul_si(tmp3, tmp2, 20, digits)
Call fpsub(accum, accum, tmp3, digits)
Call fpmul(tmp, tmp2, tmp, digits)
Call fpmul_si(tmp3, tmp, 16, digits)
Call fpadd(accum, accum, tmp3, digits)
Next i
Call fpmul(tmp, accum, accum, digits)
Call si2fp(tmp3, 1, digits)
Call fpsub(tmp, tmp3, tmp, digits)
Call fpsqr(tmp, tmp, digits)
Call fpdiv(tmp, accum, tmp, digits)
tmp.sign = sign Xor x.sign
result = tmp
End Sub
Sub set (s$, i&, v~&)
Mid$(s$, 4 * i& + 1, 4) = MKL$(v~&)
End Sub
Function git~& (s$, i&)
git = CVL(Mid$(s$, 4 * i& + 1, 4))
End Function
Function fp2str_exp$ (n As decfloat, places As Long)
Dim As Long i, ex
Dim As String v, f, ts
If n.exponent > 0 Then
ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
Else
ex = 0
End If
If n.sign Then v = "-" Else v = " "
ts = Str$(git(n.mantissa, 0))
ts = Trim$(ts)
If Len(ts) < 8 Then
ts = ts + String$(8 - Len(ts), "0")
End If
v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
For i = 1 To NUM_DWORDS - 1
ts = Str$(git(n.mantissa, i))
ts = Trim$(ts)
If Len(ts) < 8 Then
ts = String$(8 - Len(ts), "0") + ts
End If
v = v + ts
Next
v = Left$(v, places + 3)
f = Trim$(Str$(Abs(ex)))
f = String$(5 - Len(f), "0") + f
If ex < 0 Then v = v + "E-" Else v = v + "E+"
v = v + f
fp2str_exp$ = v
End Function
'long part of num
Sub fpfix (result As decfloat, num As decfloat)
Dim As decfloat ip
Dim As Long ex, ex2, j, k
ex = (num.exponent And &H7FFFFFFF) - BIAS
If ex < 1 Then
result = ip: Exit Sub
End If
If ex >= (NUM_DIGITS) Then
result = num: Exit Sub
End If
ex2 = ex \ 8
k = ex2
j = ex Mod 8
While ex2 > 0
ex2 = ex2 - 1
Call set(ip.mantissa, ex2, git(num.mantissa, ex2))
Wend
If j = 1 Then
Call set(ip.mantissa, k, 10000000 * (git(num.mantissa, k) \ 10000000))
ElseIf j = 2 Then
Call set(ip.mantissa, k, 1000000 * (git(num.mantissa, k) \ 1000000))
ElseIf j = 3 Then
Call set(ip.mantissa, k, 100000 * (git(num.mantissa, k) \ 100000))
ElseIf j = 4 Then
Call set(ip.mantissa, k, 10000 * (git(num.mantissa, k) \ 10000))
ElseIf j = 5 Then
Call set(ip.mantissa, k, 1000 * (git(num.mantissa, k) \ 1000))
ElseIf j = 6 Then
Call set(ip.mantissa, k, 100 * (git(num.mantissa, k) \ 100))
ElseIf j = 7 Then
Call set(ip.mantissa, k, 10 * (git(num.mantissa, k) \ 10))
ElseIf j = 8 Then
Call set(ip.mantissa, k, git(num.mantissa, k))
End If
ip.exponent = ex + BIAS
ip.sign = num.sign
result = ip
End Sub
Function fpfix_is_odd& (num As decfloat)
Dim As Long ex, j, k
ex = (num.exponent And &H7FFFFFFF) - BIAS
If ex < 1 Then
fpfix_is_odd = 0: Exit Function
End If
If ex >= (NUM_DIGITS) Then
Print "error in function fpfix_is_odd"
fpfix_is_odd = 99999999: Exit Function
End If
k = ex \ 8
j = ex Mod 8
If j = 1 Then
fpfix_is_odd = (git(num.mantissa, k) \ 10000000) And 1: Exit Function
ElseIf j = 2 Then
fpfix_is_odd = (git(num.mantissa, k) \ 1000000) And 1: Exit Function
ElseIf j = 3 Then
fpfix_is_odd = (git(num.mantissa, k) \ 100000) And 1: Exit Function
ElseIf j = 4 Then
fpfix_is_odd = (git(num.mantissa, k) \ 10000) And 1: Exit Function
ElseIf j = 5 Then
fpfix_is_odd = (git(num.mantissa, k) \ 1000) And 1: Exit Function
ElseIf j = 6 Then
fpfix_is_odd = (git(num.mantissa, k) \ 100) And 1: Exit Function
ElseIf j = 7 Then
fpfix_is_odd = (git(num.mantissa, k) \ 10) And 1: Exit Function
ElseIf j = 8 Then
fpfix_is_odd = git(num.mantissa, k) And 1: Exit Function
End If
fpfix_is_odd = 0
End Function
Function fp2dbl# (n As decfloat)
Dim As Long ex
Dim As String v, f, ts
If n.exponent > 0 Then
ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
Else
ex = 0
End If
If n.sign Then v = "-" Else v = " "
ts = Trim$(Str$(git(n.mantissa, 0)))
If Len(ts) < 8 Then
ts = ts + String$(8 - Len(ts), "0")
End If
v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
ts = Trim$(Str$(git(n.mantissa, 1)))
If Len(ts) < 8 Then
ts = String$(8 - Len(ts), "0") + ts
End If
v = v + ts
f = Str$(Abs(ex))
f = String$(5 - Len(f), "0") + f
If ex < 0 Then v = v + "E-" Else v = v + "E+"
v = v + f
fp2dbl# = Val(v)
End Function
Function fp2str_fix$ (n As decfloat, places As Long)
Dim As Long i, ex
Dim As String v, ts, s
If n.exponent > 0 Then
ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
Else
ex = 0
End If
If n.sign Then s = "-" Else s = " "
ts = Trim$(Str$(git(n.mantissa, 0)))
If Len(ts) < 8 Then
ts = ts + String$(8 - Len(ts), "0")
End If
v = ts
For i = 1 To NUM_DWORDS - 1
ts = Trim$(Str$(git(n.mantissa, i)))
If Len(ts) < 8 Then
ts = String$(8 - Len(ts), "0") + ts
End If
v = v + ts
Next
If places < NUM_DIGITS Then
v = Left$(v, places)
End If
If ex = 0 Then
v = Left$(v, 1) + "." + Mid$(v, 2)
ElseIf ex < 0 Then
v = "0." + String$(Abs(ex) - 1, "0") + v
ElseIf ex > 0 Then
v = Left$(v, ex + 1) + "." + Mid$(v, ex + 2)
End If
fp2str_fix$ = s + v
End Function
Function fp2str$ (n As decfloat, places As Long)
Dim As Long ex
If n.exponent > 0 Then
ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
Else
ex = 0
End If
If Abs(ex) < places Then
fp2str = fp2str_fix(n, places)
Else
fp2str = fp2str_exp(n, places)
End If
End Function
Sub str2fp (result As decfloat, value As String)
Dim As Long j, s, d, e, ep, ex, es, i, f, fp, fln
Dim As String c, f1, f2, f3, ts
Dim As Unsigned Long ulng
Dim n As decfloat
j = 1
s = 1
d = 0
e = 0
ep = 0
ex = 0
es = 1
i = 0
f = 0
fp = 0
f1 = ""
f2 = ""
f3 = ""
value = UCase$(value)
fln = Len(value)
While j <= fln
c = Mid$(value, j, 1)
If ep = 1 Then
If c = " " Then
j = j + 1
GoTo skip_while
End If
If c = "-" Then
es = -es
c = ""
End If
If c = "+" Then
j = j + 1
GoTo skip_while
End If
If (c = "0") And (f3 = "") Then
j = j + 1
GoTo skip_while
End If
If (c > "/") And (c < ":") Then 'c is digit between 0 and 9
f3 = f3 + c
ex = 10 * ex + (Asc(c) - 48)
j = j + 1
GoTo skip_while
End If
End If
If c = " " Then
j = j + 1
GoTo skip_while
End If
If c = "-" Then
s = -s
j = j + 1
GoTo skip_while
End If
If c = "+" Then
j = j + 1
GoTo skip_while
End If
If c = "." Then
If d = 1 Then
j = j + 1
GoTo skip_while
End If
d = 1
End If
If (c > "/") And (c < ":") Then 'c is digit between 0 and 9
If ((c = "0") And (i = 0)) Then
If d = 0 Then
j = j + 1
GoTo skip_while
End If
If (d = 1) And (f = 0) Then
e = e - 1
j = j + 1
GoTo skip_while
End If
End If
If d = 0 Then
f1 = f1 + c
i = i + 1
Else
If (c > "0") Then
fp = 1
End If
f2 = f2 + c
f = f + 1
End If
End If
If c = "E" Or c = "D" Then
ep = 1
End If
j = j + 1
skip_while:
Wend
If fp = 0 Then
f = 0
f2 = ""
End If
If s = -1 Then s = &H8000 Else s = 0
n.sign = s
ex = es * ex - 1 + i + e
f1 = f1 + f2
f1 = Mid$(f1, 1, 1) + Right$(f1, Len(f1) - 1)
fln = Len(f1)
If Len(f1) > (NUM_DIGITS + 1 + 8) Then
f1 = Mid$(f1, 1, (NUM_DIGITS + 1 + 8))
End If
While Len(f1) < (NUM_DIGITS + 1 + 8)
f1 = f1 + "0"
Wend
j = 1
For i = 0 To NUM_DWORDS
ts = Mid$(f1, j, 8)
ulng = Val(ts)
Call set(n.mantissa, i, ulng)
If ulng <> 0 Then fp = 1
j = j + 8
Next
If fp Then n.exponent = (ex + BIAS + 1) Else n.exponent = 0
result = n
End Sub
Sub si2fp (result As decfloat, m As Integer64, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As decfloat fac1
Dim As Long i
Dim As Integer64 n
n = Abs(m)
If n > 9999999999999999 Then
Call str2fp(fac1, Str$(m))
result = fac1: Exit Sub
End If
For i = 1 To digits
Call set(fac1.mantissa, i, 0)
Next
If m = 0 Then
fac1.exponent = 0
fac1.sign = 0
result = fac1: Exit Sub
End If
fac1.exponent = BIAS
If n < 100000000 Then
If n < 10 Then
Call set(fac1.mantissa, 0, n * 10000000)
fac1.exponent = fac1.exponent + 1
ElseIf n < 100 Then
Call set(fac1.mantissa, 0, n * 1000000)
fac1.exponent = fac1.exponent + 2
ElseIf n < 1000 Then
Call set(fac1.mantissa, 0, n * 100000)
fac1.exponent = fac1.exponent + 3
ElseIf n < 10000 Then
Call set(fac1.mantissa, 0, n * 10000)
fac1.exponent = fac1.exponent + 4
ElseIf n < 100000 Then
Call set(fac1.mantissa, 0, n * 1000)
fac1.exponent = fac1.exponent + 5
ElseIf n < 1000000 Then
Call set(fac1.mantissa, 0, n * 100)
fac1.exponent = fac1.exponent + 6
ElseIf n < 10000000 Then
Call set(fac1.mantissa, 0, n * 10)
fac1.exponent = fac1.exponent + 7
ElseIf n < 100000000 Then
Call set(fac1.mantissa, 0, n)
fac1.exponent = fac1.exponent + 8
End If
End If
If n > 99999999 Then
fac1.exponent = fac1.exponent + 8
If n < 1000000000 Then
Call set(fac1.mantissa, 0, n \ 10)
Call set(fac1.mantissa, 1, (n Mod 10) * 10000000)
fac1.exponent = fac1.exponent + 1
ElseIf n < 100000000000 Then
Call set(fac1.mantissa, 0, n \ 100)
Call set(fac1.mantissa, 1, (n Mod 100) * 1000000)
fac1.exponent = fac1.exponent + 2
ElseIf n < 1000000000000 Then
Call set(fac1.mantissa, 0, n \ 1000)
Call set(fac1.mantissa, 1, (n Mod 1000) * 100000)
fac1.exponent = fac1.exponent + 3
ElseIf n < 10000000000000 Then
Call set(fac1.mantissa, 0, n \ 10000)
Call set(fac1.mantissa, 1, (n Mod 10000) * 10000)
fac1.exponent = fac1.exponent + 4
ElseIf n < 100000000000000 Then
Call set(fac1.mantissa, 0, n \ 100000)
Call set(fac1.mantissa, 1, (n Mod 100000) * 1000)
fac1.exponent = fac1.exponent + 5
ElseIf n < 1000000000000000 Then
Call set(fac1.mantissa, 0, n \ 1000000)
Call set(fac1.mantissa, 1, (n Mod 1000000) * 100)
fac1.exponent = fac1.exponent + 6
ElseIf n < 10000000000000000 Then
Call set(fac1.mantissa, 0, n \ 10000000)
Call set(fac1.mantissa, 1, (n Mod 10000000) * 10)
fac1.exponent = fac1.exponent + 7
ElseIf n < 100000000000000000 Then
Call set(fac1.mantissa, 0, n \ 100000000)
Call set(fac1.mantissa, 1, n Mod 100000000)
fac1.exponent = fac1.exponent + 8
End If
End If
If m < 0 Then
fac1.sign = &H8000
Else
fac1.sign = 0
End If
result = fac1
End Sub
Sub RSHIFT_1 (mantissa As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Unsigned Long v1, v2
Dim As Long i
For i = digits To 1 Step -1
v1 = git(mantissa.mantissa, i) \ 10
v2 = git(mantissa.mantissa, i - 1) Mod 10
v2 = v2 * 10000000 + v1
Call set(mantissa.mantissa, i, v2)
Next
Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 10)
End Sub
Sub LSHIFT_1 (mantissa As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Unsigned Long v1, v2
Dim As Long i
For i = 0 To digits - 1
v1 = git(mantissa.mantissa, i) Mod 10000000
v2 = git(mantissa.mantissa, i + 1) \ 10000000
Call set(mantissa.mantissa, i, v1 * 10 + v2)
Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 10000000)
Next
Call set(mantissa.mantissa, digits, 10 * (git(mantissa.mantissa, digits) Mod 10000000))
End Sub
Sub RSHIFT_2 (mantissa As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Unsigned Long v1, v2
Dim As Long i
For i = digits To 1 Step -1
v1 = git(mantissa.mantissa, i) \ 100
v2 = git(mantissa.mantissa, i - 1) Mod 100
v2 = v2 * 1000000 + v1
Call set(mantissa.mantissa, i, v2)
Next
Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 100)
End Sub
Sub LSHIFT_2 (mantissa As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Unsigned Long v1, v2
Dim As Long i
For i = 0 To digits - 1
v1 = git(mantissa.mantissa, i) Mod 1000000
v2 = git(mantissa.mantissa, i + 1) \ 1000000
Call set(mantissa.mantissa, i, v1 * 100 + v2)
Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 1000000)
Next
Call set(mantissa.mantissa, digits, 100 * (git(mantissa.mantissa, digits) Mod 1000000))
End Sub
Sub RSHIFT_3 (mantissa As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Unsigned Long v1, v2
Dim As Long i
For i = digits To 1 Step -1
v1 = git(mantissa.mantissa, i) \ 1000
v2 = git(mantissa.mantissa, i - 1) Mod 1000
v2 = v2 * 100000 + v1
Call set(mantissa.mantissa, i, v2)
Next
Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 1000)
End Sub
Sub LSHIFT_3 (mantissa As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Unsigned Long v1, v2
Dim As Long i
For i = 0 To digits - 1
v1 = git(mantissa.mantissa, i) Mod 100000
v2 = git(mantissa.mantissa, i + 1) \ 100000
Call set(mantissa.mantissa, i, v1 * 1000 + v2)
Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 100000)
Next
Call set(mantissa.mantissa, digits, 1000 * (git(mantissa.mantissa, digits) Mod 100000))
End Sub
Sub RSHIFT_4 (mantissa As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Unsigned Long v1, v2
Dim As Long i
For i = digits To 1 Step -1
v1 = git(mantissa.mantissa, i) \ 10000
v2 = git(mantissa.mantissa, i - 1) Mod 10000
v2 = v2 * 10000 + v1
Call set(mantissa.mantissa, i, v2)
Next
Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 10000)
End Sub
Sub LSHIFT_4 (mantissa As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Unsigned Long v1, v2
Dim As Long i
For i = 0 To digits - 1
v1 = git(mantissa.mantissa, i) Mod 10000
v2 = git(mantissa.mantissa, i + 1) \ 10000
Call set(mantissa.mantissa, i, v1 * 10000 + v2)
Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 10000)
Next
Call set(mantissa.mantissa, digits, 10000 * (git(mantissa.mantissa, digits) Mod 10000))
End Sub
Sub RSHIFT_5 (mantissa As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Unsigned Long v1, v2
Dim As Long i
For i = digits To 1 Step -1
v1 = git(mantissa.mantissa, i) \ 100000
v2 = git(mantissa.mantissa, i - 1) Mod 100000
v2 = v2 * 1000 + v1
Call set(mantissa.mantissa, i, v2)
Next
Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 100000)
End Sub
Sub LSHIFT_5 (mantissa As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Unsigned Long v1, v2
Dim As Long i
For i = 0 To digits - 1
v1 = git(mantissa.mantissa, i) Mod 1000
v2 = git(mantissa.mantissa, i + 1) \ 1000
Call set(mantissa.mantissa, i, v1 * 100000 + v2)
Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 1000)
Next
Call set(mantissa.mantissa, digits, 100000 * (git(mantissa.mantissa, digits) Mod 1000))
End Sub
Sub RSHIFT_6 (mantissa As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Unsigned Long v1, v2
Dim As Long i
For i = digits To 1 Step -1
v1 = git(mantissa.mantissa, i) \ 1000000
v2 = git(mantissa.mantissa, i - 1) Mod 1000000
v2 = v2 * 100 + v1
Call set(mantissa.mantissa, i, v2)
Next
Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 1000000)
End Sub
Sub LSHIFT_6 (mantissa As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Unsigned Long v1, v2
Dim As Long i
For i = 0 To digits - 1
v1 = git(mantissa.mantissa, i) Mod 100
v2 = git(mantissa.mantissa, i + 1) \ 100
Call set(mantissa.mantissa, i, v1 * 1000000 + v2)
Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 100)
Next
Call set(mantissa.mantissa, digits, 1000000 * (git(mantissa.mantissa, digits) Mod 100))
End Sub
Sub RSHIFT_7 (mantissa As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Unsigned Long v1, v2
Dim As Long i
For i = digits To 1 Step -1
v1 = git(mantissa.mantissa, i) \ 10000000
v2 = git(mantissa.mantissa, i - 1) Mod 10000000
v2 = v2 * 10 + v1
Call set(mantissa.mantissa, i, v2)
Next
Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 10000000)
End Sub
Sub LSHIFT_7 (mantissa As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Unsigned Long v1, v2
Dim As Long i
For i = 0 To digits - 1
v1 = git(mantissa.mantissa, i) Mod 10
v2 = git(mantissa.mantissa, i + 1) \ 10
Call set(mantissa.mantissa, i, v1 * 10000000 + v2)
Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 10)
Next
Call set(mantissa.mantissa, digits, 10000000 * (git(mantissa.mantissa, digits) Mod 10))
End Sub
Sub RSHIFT_8 (mantissa As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Long i
For i = digits To 1 Step -1
Call set(mantissa.mantissa, i, git(mantissa.mantissa, i - 1))
Next
Call set(mantissa.mantissa, 0, 0)
End Sub
Sub LSHIFT_8 (mantissa As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Long i
For i = 0 To digits - 1
Call set(mantissa.mantissa, i, git(mantissa.mantissa, i + 1))
Next
Call set(mantissa.mantissa, digits, 0)
End Sub
Function fpcmp& (x As decfloat, y As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Long c, i
If x.sign < y.sign Then fpcmp = -1: Exit Function
If x.sign > y.sign Then fpcmp = 1: Exit Function
If x.sign = y.sign Then
If x.exponent = y.exponent Then
For i = 0 To digits
c = git(x.mantissa, i) - git(y.mantissa, i)
If c <> 0 Then Exit For
Next
If c < 0 Then fpcmp = -1: Exit Function
If c = 0 Then fpcmp = 0: Exit Function
If c > 0 Then fpcmp = 1: Exit Function
End If
If x.exponent < y.exponent Then fpcmp = -1: Exit Function
If x.exponent > y.exponent Then fpcmp = 1: Exit Function
End If
' if we reach this point it means that the signs are different
' and if the sign of x is set meaning that x is negative then x < y
End Function
Sub NORM_FAC1 (fac1 As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
' normalize the number in fac1
' all routines exit through this one.
'see if the mantissa is all zeros.
'if so, set the exponent and sign equal to 0.
Dim As Long i, er, f
er = 0: f = 0
For i = 0 To digits
If git(fac1.mantissa, i) > 0 Then f = 1
Next
If f = 0 Then
fac1.exponent = 0
fac1.sign = 0
Exit Sub
'if the highmost Digit in fac1_man is nonzero,
'shift the mantissa right 1 Digit and
'increment the exponent
ElseIf git(fac1.mantissa, 0) > 99999999 Then
Call RSHIFT_1(fac1, digits)
fac1.exponent = fac1.exponent + 1
Else
'now shift fac1_man 1 to the left until a
'nonzero digit appears in the next-to-highest
'Digit of fac1_man. decrement exponent for
'each shift.
While git(fac1.mantissa, 0) = 0
Call LSHIFT_8(fac1, digits)
fac1.exponent = fac1.exponent - 8
If fac1.exponent = 0 Then
Print "NORM_FAC1=EXPU_ERR"
Exit Sub
End If
Wend
If git(fac1.mantissa, 0) < 10 Then
Call LSHIFT_7(fac1, digits)
fac1.exponent = fac1.exponent - 7
ElseIf git(fac1.mantissa, 0) < 100 Then
Call LSHIFT_6(fac1, digits)
fac1.exponent = fac1.exponent - 6
ElseIf git(fac1.mantissa, 0) < 1000 Then
Call LSHIFT_5(fac1, digits)
fac1.exponent = fac1.exponent - 5
ElseIf git(fac1.mantissa, 0) < 10000 Then
Call LSHIFT_4(fac1, digits)
fac1.exponent = fac1.exponent - 4
ElseIf git(fac1.mantissa, 0) < 100000 Then
Call LSHIFT_3(fac1, digits)
fac1.exponent = fac1.exponent - 3
ElseIf git(fac1.mantissa, 0) < 1000000 Then
Call LSHIFT_2(fac1, digits)
fac1.exponent = fac1.exponent - 2
ElseIf git(fac1.mantissa, 0) < 10000000 Then
Call LSHIFT_1(fac1, digits)
fac1.exponent = fac1.exponent - 1
End If
End If
'check for overflow/underflow
If fac1.exponent < 0 Then
Print "NORM_FAC1=EXPO_ERR"
End If
End Sub
Sub fpadd_aux (fac1 As decfloat, fac2 As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Unsigned Long v, c, i
c = 0
For i = digits To 1 Step -1
v = git(fac2.mantissa, i) + git(fac1.mantissa, i) + c
If v > 99999999 Then
v = v - 100000000
c = 1
Else
c = 0
End If
Call set(fac1.mantissa, i, v)
Next
v = git(fac1.mantissa, 0) + git(fac2.mantissa, 0) + c
Call set(fac1.mantissa, 0, v)
Call NORM_FAC1(fac1, digits)
End Sub
Sub fpsub_aux (fac1 As decfloat, fac2 As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Long v, c, i
c = 0
For i = digits To 1 Step -1
v = git(fac1.mantissa, i) - git(fac2.mantissa, i) - c
If v < 0 Then
v = v + 100000000
c = 1
Else
c = 0
End If
Call set(fac1.mantissa, i, v)
Next
v = git(fac1.mantissa, 0) - git(fac2.mantissa, 0) - c
Call set(fac1.mantissa, 0, v)
Call NORM_FAC1(fac1, digits)
End Sub
Sub fpadd (result As decfloat, x As decfloat, y As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As decfloat fac1, fac2
Dim As Long i, t, c, xsign, ysign
xsign = x.sign: x.sign = 0
ysign = y.sign: y.sign = 0
c = fpcmp(x, y, digits)
x.sign = xsign
y.sign = ysign
If c < 0 Then
fac1 = y
fac2 = x
Else
fac1 = x
fac2 = y
End If
t = fac1.exponent - fac2.exponent
t = ((fac1.exponent And &H7FFFFFFF) - BIAS - 1) - ((fac2.exponent And &H7FFFFFFF) - BIAS - 1)
If t < (NUM_DIGITS + 8) Then
'The difference between the two
'exponents indicate how many times
'we have to multiply the mantissa
'of FAC2 by 10 (i.e., shift it right 1 place).
'If we have to shift more times than
'we have digits, the result is already in FAC1.
t = fac1.exponent - fac2.exponent
If t > 0 And t < (NUM_DIGITS + 8) Then 'shift
i = t \ 8
While i > 0
Call RSHIFT_8(fac2, digits)
t = t - 8
i = i - 1
Wend
If t = 7 Then
Call RSHIFT_7(fac2, digits)
ElseIf t = 6 Then
Call RSHIFT_6(fac2, digits)
ElseIf t = 5 Then
Call RSHIFT_5(fac2, digits)
ElseIf t = 4 Then
Call RSHIFT_4(fac2, digits)
ElseIf t = 3 Then
Call RSHIFT_3(fac2, digits)
ElseIf t = 2 Then
Call RSHIFT_2(fac2, digits)
ElseIf t = 1 Then
Call RSHIFT_1(fac2, digits)
End If
End If
'See if the signs of the two numbers
'are the same. If so, add; if not, subtract.
If fac1.sign = fac2.sign Then 'add
Call fpadd_aux(fac1, fac2, digits)
Else
Call fpsub_aux(fac1, fac2, digits)
End If
End If
result = fac1
End Sub
Sub fpsub (result As decfloat, x As decfloat, y As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As decfloat fac1, fac2
fac1 = x
fac2 = y
fac2.sign = fac2.sign Xor &H8000
Call fpadd(fac1, fac1, fac2, digits)
result = fac1
End Sub
Sub fpmul (result As decfloat, x As decfloat, y As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As decfloat fac1, fac2
Dim As Long i, j, ex, er, den, num
Dim As Integer64 digit, carry
Dim As Unsigned Integer64 fac3(0 To digits)
fac1 = x
fac2 = y
'check exponents. if either is zero,
'the result is zero
If fac1.exponent = 0 Or fac2.exponent = 0 Then 'result is zero...clear fac1.
fac1.sign = 0
fac1.exponent = 0
For i = 0 To digits
Call set(fac1.mantissa, i, 0)
Next
'NORM_FAC1(fac1)
result = fac1: Exit Sub
Else
If ex < 0 Then
er = EXPO_ERR
result = fac1: Exit Sub
End If
'clear fac3 mantissa
For i = 0 To digits
fac3(i) = 0
Next
den = digits
While git(fac2.mantissa, den) = 0
den = den - 1
Wend
num = digits
While git(fac1.mantissa, num) = 0
num = num - 1
Wend
If num < den Then
Swap fac1, fac2
'fac1=y
'fac2=x
Swap den, num
End If
For j = den To 0 Step -1
carry = 0
digit = git(fac2.mantissa, j)
For i = num To 0 Step -1
fac3(i) = fac3(i) + digit * git(fac1.mantissa, i)
Next
For i = num To 0 Step -1
fac3(i) = fac3(i) + carry
carry = fac3(i) \ 100000000
fac3(i) = (fac3(i) Mod 100000000)
Next
For i = digits To 1 Step -1
fac3(i) = fac3(i - 1)
Next
fac3(0) = carry
Next
For i = 0 To digits
Call set(fac1.mantissa, i, fac3(i))
Next
End If
'now determine exponent of result.
'as you do...watch for overflow.
ex = fac2.exponent - BIAS + fac1.exponent
fac1.exponent = ex
'determine the sign of the product
fac1.sign = fac1.sign Xor fac2.sign
Call NORM_FAC1(fac1, digits)
result = fac1
End Sub
Sub fpmul_si (result As decfloat, x As decfloat, y As Integer64, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As decfloat fac1, fac2
Dim As Long count, ex, er, i
Dim As Integer64 carry, digit, prod, value
fac1 = x
digit = Abs(y)
If digit > 99999999 Then
Call si2fp(fac2, y, digits)
Call fpmul(fac1, fac1, fac2, digits)
result = fac1: Exit Sub
End If
'check exponents. if either is zero,
'the result is zero
If fac1.exponent = 0 Or y = 0 Then 'result is zero...clear fac1.
fac1.sign = 0
fac1.exponent = 0
For count = 0 To digits
Call set(fac1.mantissa, count, 0)
Next
Call NORM_FAC1(fac1, digits)
result = fac1: Exit Sub
Else
If digit = 1 Then
If y < 0 Then
fac1.sign = fac1.sign Xor &H8000
End If
result = fac1: Exit Sub
End If
'now determine exponent of result.
'as you do...watch for overflow.
If ex < 0 Then
er = EXPO_ERR
result = fac1: Exit Sub
End If
carry = 0
For i = digits To 0 Step -1
prod = digit * git(fac1.mantissa, i) + carry
value = (prod Mod 100000000)
Call set(fac1.mantissa, i, value)
carry = prod \ 100000000
Next
If carry < 10 Then
Call RSHIFT_1(fac1, digits)
fac1.exponent = fac1.exponent + 1
Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10000000)
ElseIf carry < 100 Then
Call RSHIFT_2(fac1, digits)
fac1.exponent = fac1.exponent + 2
Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 1000000)
ElseIf carry < 1000 Then
Call RSHIFT_3(fac1, digits)
fac1.exponent = fac1.exponent + 3
Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 100000)
ElseIf carry < 10000 Then
Call RSHIFT_4(fac1, digits)
fac1.exponent = fac1.exponent + 4
Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10000)
ElseIf carry < 100000 Then
Call RSHIFT_5(fac1, digits)
fac1.exponent = fac1.exponent + 5
Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 1000)
ElseIf carry < 1000000 Then
Call RSHIFT_6(fac1, digits)
fac1.exponent = fac1.exponent + 6
Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 100)
ElseIf carry < 10000000 Then
Call RSHIFT_7(fac1, digits)
fac1.exponent = fac1.exponent + 7
Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10)
ElseIf carry < 100000000 Then
Call RSHIFT_8(fac1, digits)
fac1.exponent = fac1.exponent + 8
Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry)
End If
End If
Call NORM_FAC1(fac1, digits)
If y < 0 Then
fac1.sign = fac1.sign Xor &H8000
End If
result = fac1
End Sub
Sub fpmul_ll (result As decfloat, x As decfloat, y As Integer64, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As decfloat fac1, fac2
Dim As Long count, ex, er, i
Dim As Integer64 carry, digit, prod, value, n0, n1
fac1 = x
digit = Abs(y)
If digit > 99999999 And digit < 10000000000000000 Then
n0 = digit \ 100000000
n1 = digit Mod 100000000
Call fpmul_si(fac2, fac1, n1, digits)
fac1.exponent = fac1.exponent + 8
Call fpmul_si(fac1, fac1, n0, digits)
Call fpadd(fac1, fac1, fac2, digits)
If y < 0 Then fac1.sign = fac1.sign Xor &H8000
result = fac1: Exit Sub
End If
If digit > 9999999999999999 Then
Call str2fp(fac2, Str$(y))
Call fpmul(fac1, fac1, fac2, digits)
result = fac1: Exit Sub
End If
Call si2fp(fac2, y, digits)
'check exponents. if either is zero,
'the result is zero
If fac1.exponent = 0 Or y = 0 Then 'result is zero...clear fac1.
fac1.sign = 0
fac1.exponent = 0
For count = 0 To digits
Call set(fac1.mantissa, count, 0)
Next
Call NORM_FAC1(fac1, digits)
result = fac1: Exit Sub
Else
If digit = 1 Then
If y < 0 Then
fac1.sign = fac1.sign Xor &H8000
End If
result = fac1: Exit Sub
End If
'now determine exponent of result.
'as you do...watch for overflow.
If ex < 0 Then
er = EXPO_ERR
result = fac1: Exit Sub
End If
carry = 0
For i = digits To 0 Step -1
prod = digit * git(fac1.mantissa, i) + carry
value = (prod Mod 100000000)
Call set(fac1.mantissa, i, value)
carry = prod \ 100000000
Next
If carry < 10 Then
Call RSHIFT_1(fac1, digits)
fac1.exponent = fac1.exponent + 1
Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10000000)
ElseIf carry < 100 Then
Call RSHIFT_2(fac1, digits)
fac1.exponent = fac1.exponent + 2
Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 1000000)
ElseIf carry < 1000 Then
Call RSHIFT_3(fac1, digits)
fac1.exponent = fac1.exponent + 3
Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 100000)
ElseIf carry < 10000 Then
Call RSHIFT_4(fac1, digits)
fac1.exponent = fac1.exponent + 4
Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10000)
ElseIf carry < 100000 Then
Call RSHIFT_5(fac1, digits)
fac1.exponent = fac1.exponent + 5
Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 1000)
ElseIf carry < 1000000 Then
Call RSHIFT_6(fac1, digits)
fac1.exponent = fac1.exponent + 6
Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 100)
ElseIf carry < 10000000 Then
Call RSHIFT_7(fac1, digits)
fac1.exponent = fac1.exponent + 7
Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10)
ElseIf carry < 100000000 Then
Call RSHIFT_8(fac1, digits)
fac1.exponent = fac1.exponent + 8
Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry)
End If
End If
Call NORM_FAC1(fac1, digits)
If y < 0 Then
fac1.sign = fac1.sign Xor &H8000
End If
result = fac1
End Sub
Sub fpinv (result As decfloat, m As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As Double x
Dim As Long k, l, ex
Dim As Long prec
Dim As decfloat r, r2, one, n
n = m
l = Log((NUM_DIGITS + 8) * 0.0625) * 1.5
Dim As String v, ts
If n.exponent > 0 Then
ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
Else
ex = 0
End If
If n.sign Then v = "-" Else v = " "
ts = Str$(git(n.mantissa, 0))
If Len(ts) < 8 Then
ts = ts + String$(8 - Len(ts), "0")
End If
v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
ts = Str$(git(n.mantissa, 1))
If Len(ts) < 8 Then
ts = String$(8 - Len(ts), "0") + ts
End If
v = v + ts
x = Val(v)
If x = 0 Then Print "Div 0": result = r: Exit Sub
If x = 1 And ex = 0 Then
Call str2fp(r, "1")
result = r: Exit Sub
ElseIf x = 1 Then
x = 10
ex = ex - 1
End If
ex = (-1) - ex
x = 1 / x
Call str2fp(r, Str$(x))
r.exponent = ex + BIAS + 1
Call set(one.mantissa, 0, 10000000)
one.exponent = BIAS + 1
r2 = r
prec = 3
For k = 1 To l
prec = 2 * prec - 1
Call fpmul(r2, n, r, prec)
Call fpsub(r2, one, r2, prec)
Call fpmul(r2, r, r2, prec)
Call fpadd(r, r, r2, prec)
Next
result = r
End Sub
Function min& (a As Long, b As Long)
If a < b Then min& = a Else min& = b
End Function
Function RealW# (w() As Double, j As Long)
Dim wx As Double
wx = ((w(j - 1) * 10000 + w(j)) * 10000 + w(j + 1)) * 10000
If UBound(w) >= (j + 2) Then wx = wx + w(j + 2)
RealW# = wx
End Function
Sub subtract (w() As Double, q As Long, d() As Double, ka As Long, kb As Long)
Dim As Long j
For j = ka To kb
w(j) = w(j) - q * d(j - ka + 2)
Next
End Sub
Sub normalize (w() As Double, ka As Long, q As Long)
w(ka) = w(ka) + w(ka - 1) * 10000
w(ka - 1) = q
End Sub
Sub finalnorm (w() As Double, kb As Long)
Dim As Long carry, j
For j = kb To 3 Step -1
If w(j) < 0 Then
carry = ((-w(j) - 1) \ 10000) + 1
Else
If w(j) >= 10000 Then
carry = -(w(j) \ 10000)
Else
carry = 0
End If
End If
w(j) = w(j) + carry * 10000
w(j - 1) = w(j - 1) - carry
Next
End Sub
Sub fpdiv (result As decfloat, x As decfloat, y As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As decfloat fac1, fac2
Dim As Long i, er, is_power_of_ten
fac1 = x
fac2 = y
If fac2.exponent = 0 Then ' if fac2 = 0, return
' a divide-by-zero error and
' bail out.
For i = 0 To digits
Call set(fac1.mantissa, i, 99999999)
Next
fac1.exponent = 99999 + BIAS + 1
er = DIVZ_ERR
result = fac1
Exit Sub
ElseIf fac1.exponent = 0 Then 'fact1=0, just return
er = 0
result = fac1
Exit Sub
Else
'check to see if fac2 is a power of ten
is_power_of_ten = 0
If git(fac2.mantissa, 0) = 10000000 Then
is_power_of_ten = 1
For i = 1 To digits
If git(fac2.mantissa, i) <> 0 Then
is_power_of_ten = 0
Exit For
End If
Next
End If
'if fac2 is a power of ten then all we need to do is to adjust the sign and exponent and we are finished
If is_power_of_ten = 1 Then
fac1.sign = fac1.sign Xor fac2.sign
fac1.exponent = fac1.exponent - fac2.exponent + BIAS + 1
result = fac1
Exit Sub
End If
Dim As Double result(1 To 2 * digits + 3), n(1 To 2 * digits + 3), d(1 To 2 * digits + 3)
Const b = 10000
Dim As Long j, last, laststep, q, t
Dim As Long stp
Dim As Double xd, xn, rund
Dim As Double w(1 To UBound(n) + 4)
For j = 0 To digits
n(2 * j + 2) = git(fac1.mantissa, j) \ 10000
n(2 * j + 3) = git(fac1.mantissa, j) Mod 10000
d(2 * j + 2) = git(fac2.mantissa, j) \ 10000
d(2 * j + 3) = git(fac2.mantissa, j) Mod 10000
Next
n(1) = (fac1.exponent And &H7FFFFFFF) - BIAS - 1
d(1) = (fac2.exponent And &H7FFFFFFF) - BIAS - 1
For j = UBound(n) To UBound(w)
w(j) = 0
Next
t = UBound(n) - 1
w(1) = n(1) - d(1) + 1
w(2) = 0
For j = 2 To UBound(n)
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, UBound(w))
Call subtract(w(), q, d(), (stp + 2), last)
Call normalize(w(), (stp + 2), q)
Next
Call finalnorm(w(), (laststep + 1))
If w(2) <> 0 Then laststep = laststep - 1
rund = w(laststep + 1) / b
If rund >= 0.5 Then w(laststep) = w(laststep) + 1
If w(2) = 0 Then
For j = 1 To t + 1
result(j) = w(j + 1)
Next
Else
For j = 1 To t + 1
result(j) = w(j)
Next
End If
If w(2) = 0 Then result(1) = w(1) - 1 Else result(1) = w(1)
For j = 0 To digits
Call set(fac1.mantissa, j, result(2 * j + 2) * 10000 + result(2 * j + 3))
Next
Call NORM_FAC1(fac1, digits)
fac1.exponent = (result(1) + BIAS)
End If
fac1.sign = fac1.sign Xor fac2.sign
result = fac1
End Sub
' sqrt(num)
Sub fpsqr (result As decfloat, num As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As decfloat r, r2, tmp, n, half
Dim As Long ex, k, l, prec
Dim As String ts, v
Dim As Double x
l = Log((NUM_DIGITS + 8) * 0.0625) * 1.5
'l=estimated number of iterations needed
'first estimate is accurate to about 16 digits
'l is approximatly = to log2((NUM_DIGITS+9)/16)
'NUM_DIGITS+9 because decfloat has an extra 9 guard digits
n = num
Call si2fp(tmp, 0, digits)
If fpcmp(n, tmp, digits) = 0 Then
Call si2fp(r, 0, digits)
result = r
Exit Sub
End If
Call si2fp(tmp, 1, digits)
If fpcmp(n, tmp, digits) = 0 Then
Call si2fp(r, 1, digits)
result = r
Exit Sub
End If
Call si2fp(tmp, 0, digits)
If fpcmp(n, tmp, digits) < 0 Then
Call si2fp(r, 0, digits)
result = r
Exit Sub
End If
'=====================================================================
'hack to bypass the limitation of double exponent range
'in case the number is larger than what a double can handle
'for example, if the number is 2e500
'we separate the exponent and mantissa in this case 2
'if the exponent is odd then multiply the mantissa by 10
'take the square root and assign it to decfloat
'divide the exponent in half for square root
'in this case 1.414213562373095e250
If n.exponent > 0 Then
ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
Else
ex = 0
End If
ts = Trim$(Str$(git(n.mantissa, 0)))
If Len(ts) < 8 Then
ts = ts + String$(8 - Len(ts), "0")
End If
v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
ts = Trim$(Str$(git(n.mantissa, 1)))
If Len(ts) < 8 Then
ts = String$(8 - Len(ts), "0") + ts
End If
v = v + ts
x = Val(v)
If x = 0 Then Print "Div 0": result = r: Exit Sub
If x = 1 And ex = 0 Then
Call si2fp(r, 1, digits)
result = r
Exit Sub
End If
If Abs(ex) And 1 Then
x = x * 10
ex = ex - 1
End If
x = Sqr(x) 'approximation
v = Trim$(Str$(x))
k = InStr(v, ".")
Call str2fp(r, v)
r.exponent = ex \ 2 + BIAS + 1
If Len(v) > 1 And k = 0 Then r.exponent = r.exponent + 1
For k = 1 To digits
Call set(half.mantissa, k, 0)
Next
Call set(half.mantissa, 0, 50000000)
half.exponent = BIAS
half.sign = 0
'=====================================================================
'Newton-Raphson method
prec = 3
For k = 1 To l + 1
prec = 2 * prec - 1
Call fpdiv(tmp, n, r, prec)
Call fpadd(r2, r, tmp, prec)
Call fpmul(r, r2, half, prec)
Next
result = r
End Sub
Sub fpdiv_si (result As decfloat, num As decfloat, den As Long, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As decfloat fac1
Dim As Unsigned Integer64 carry, remder
Dim As Integer64 i, divisor
Dim As Integer64 quotient
remder = 0
divisor = Abs(den)
fac1 = num
If divisor = 0 Then
Print "error: divisor = 0"
result = fac1: Exit Sub
End If
If divisor > 99999999 Then
Print "error: divisor too large"
result = fac1: Exit Sub
End If
For i = 0 To digits
quotient = git(fac1.mantissa, i) + remder * 100000000
remder = quotient Mod divisor
Call set(fac1.mantissa, i, quotient \ divisor)
Next
quotient = remder * 100000000
quotient = quotient \ divisor
carry = git(fac1.mantissa, 0)
If carry = 0 Then
Call LSHIFT_8(fac1, digits)
fac1.exponent = fac1.exponent - 8
Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient)
ElseIf carry < 10 Then
Call LSHIFT_7(fac1, digits)
fac1.exponent = fac1.exponent - 7
Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 10)
ElseIf carry < 100 Then
Call LSHIFT_6(fac1, digits)
fac1.exponent = fac1.exponent - 6
Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 100)
ElseIf carry < 1000 Then
Call LSHIFT_5(fac1, digits)
fac1.exponent = fac1.exponent - 5
Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 1000)
ElseIf carry < 10000 Then
Call LSHIFT_4(fac1, digits)
fac1.exponent = fac1.exponent - 4
Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 10000)
ElseIf carry < 100000 Then
Call LSHIFT_3(fac1, digits)
fac1.exponent = fac1.exponent - 3
Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 100000)
ElseIf carry < 1000000 Then
Call LSHIFT_2(fac1, digits)
fac1.exponent = fac1.exponent - 2
Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 1000000)
ElseIf carry < 10000000 Then
Call LSHIFT_1(fac1, digits)
fac1.exponent = fac1.exponent - 1
Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 10000000)
End If
'NORM_FAC1(fac1)
If den < 0 Then
fac1.sign = fac1.sign Xor &H8000
End If
result = fac1
End Sub
'fractional part of num
Sub fpfrac (result As decfloat, num As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As decfloat n
Call fpfix(n, num)
Call fpsub(n, num, n, digits)
result = n
End Sub
'returns the positive of n
Sub fpabs (result As decfloat, n As decfloat)
Dim As decfloat x
x = n
x.sign = 0
result = x
End Sub
'changes the sign of n, if n is positive then n will be negative & vice versa
Sub fpneg (result As decfloat, n As decfloat)
Dim As decfloat x
x = n
x.sign = x.sign Xor &H8000
result = x
End Sub
'returns the negative of n regardless of the sign of n
Sub fpnegative (result As decfloat, n As decfloat)
Dim As decfloat x
x = n
x.sign = &H8000
result = x
End Sub
Sub fpfmod (quotient As decfloat, f_mod As decfloat, num As decfloat, denom As decfloat, digits As Long)
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As decfloat q, fm 'make copy in case the destination and source are the same
Call fpdiv(fm, num, denom, digits): Call fpfix(q, fm)
Call fpsub(fm, fm, q, digits): Call fpmul(fm, fm, denom, digits)
quotient = q
f_mod = fm
End Sub
Sub fpeps (result As decfloat, digits As Long)
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As decfloat ep
Call si2fp(ep, 1, digits)
ep.exponent = (-(NUM_DIGITS) + BIAS + 1)
result = ep
End Sub
Sub fpipow (result As decfloat, x As decfloat, e As Integer64, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
'take x to an Long power
Dim As decfloat y
Dim As decfloat z
Dim As Integer64 n, c
Dim As Long i
c = 0
y = x
n = Abs(e)
z.sign = 0
z.exponent = (BIAS + 1)
Call set(z.mantissa, 0, 10000000)
For i = 1 To NUM_DWORDS
Call set(z.mantissa, i, 0)
Next
While n > 0
While (n And 1) = 0
n = n \ 2
Call fpmul(y, y, y, digits)
c = c + 1
Wend
n = n - 1
Call fpmul(z, y, z, digits)
c = c + 1
Wend
If e < 0 Then
Call fpinv(z, z, digits)
End If
result = z
End Sub
Sub fpfactorial (result As decfloat, n As Long, digits_in As Long)
Dim As Long digits
digits = digits_in
Dim As Long i
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As decfloat f
If n < 0 Then Print "inf": result = f: Exit Sub
If n = 0 Or n = 1 Then
Call si2fp(f, 1, digits)
result = f: Exit Sub
End If
Call si2fp(f, 2, digits)
If n = 2 Then result = f: Exit Sub
For i = 3 To n
Call fpmul_si(f, f, i, digits)
Next
result = f
End Sub
Sub fplogTaylor (result As decfloat, x As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
'taylor series
'====================Log Guard==================
Dim As decfloat g, zero
Dim As Long i
Call fpabs(g, x)
zero.sign = 0
zero.exponent = 0
For i = 0 To NUM_DWORDS
Call set(zero.mantissa, i, 0)
Next
If fpcmp(g, x, digits) <> 0 Then result = zero: Exit Sub
If fpcmp(x, zero, digits) = 0 Then result = zero: Exit Sub
'=============================================
Dim As Long invflag
Dim As decfloat XX, Term, Accum, x9, tmp, tmp2
Dim As decfloat T, B, one, Q, two
one.sign = 0
one.exponent = (BIAS + 1)
Call set(one.mantissa, 0, 10000000)
two.sign = 0
two.exponent = (BIAS + 1)
Call set(two.mantissa, 0, 20000000)
For i = 1 To NUM_DWORDS
Call set(one.mantissa, i, 0)
Call set(two.mantissa, i, 0)
Next
x9 = x
If fpcmp(x, one, digits) < 0 Then
invflag = 1
Call fpdiv(x9, one, x9, digits)
End If
Call fpsub(T, x9, one, digits)
Call fpadd(B, x9, one, digits)
Call fpdiv(Accum, T, B, digits)
Call fpdiv(Q, T, B, digits)
tmp = Q
Call fpmul(XX, Q, Q, digits)
Dim As Long c
c = 1
Do
c = c + 2
tmp2 = tmp
Call fpmul(Q, Q, XX, digits)
Call fpdiv_si(Term, Q, c, digits)
Call fpadd(Accum, tmp, Term, digits)
Swap tmp, Accum
Loop Until fpcmp(tmp, tmp2, digits) = 0
Call fpmul_si(Accum, Accum, 2, digits)
If invflag Then
Call fpneg(Accum, Accum)
result = Accum: Exit Sub
End If
result = Accum
End Sub
Sub fplog (result As decfloat, x As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
'====================Log Guard==================
Dim As decfloat g, one, zero
Dim As Long i, factor
zero.sign = 0
zero.exponent = 0
For i = 0 To NUM_DWORDS
Call set(zero.mantissa, i, 0)
Next
Call si2fp(one, 1, digits)
Call fpabs(g, x)
If fpcmp(g, x, digits) <> 0 Then result = zero: Exit Sub
If fpcmp(x, zero, digits) = 0 Then result = zero: Exit Sub
If fpcmp(x, one, digits) = 0 Then result = zero: Exit Sub
'=============================================
Dim As decfloat approx, ans, logx
approx = x
factor = 4096
Call fpsqr(approx, approx, digits)
Call fpsqr(approx, approx, digits)
Call fpsqr(approx, approx, digits)
Call fpsqr(approx, approx, digits)
Call fpsqr(approx, approx, digits)
Call fpsqr(approx, approx, digits)
Call fpsqr(approx, approx, digits)
Call fpsqr(approx, approx, digits)
Call fpsqr(approx, approx, digits)
Call fpsqr(approx, approx, digits)
Call fpsqr(approx, approx, digits)
Call fpsqr(approx, approx, digits)
Call fplogTaylor(logx, approx, digits)
Call fpmul_si(ans, logx, factor, digits)
result = ans
End Sub
Sub fpexp (result As decfloat, x As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
'taylor series
Dim As decfloat fac, x9, temp, accum, p, term
Dim As Long i, c
Call si2fp(temp, 0, digits)
fac.sign = 0
fac.exponent = (BIAS + 1)
Call set(fac.mantissa, 0, 10000000)
For i = 1 To NUM_DWORDS
Call set(fac.mantissa, i, 0)
Next
If fpcmp(x, temp, digits) = 0 Then result = fac: Exit Sub
c = 1
Call fpdiv_si(x9, x, 8192, digits) 'fpdiv_si(x, 67108864) '
p = x9
Call fpadd(accum, fac, x9, digits) '1 + x
Do
c = c + 1
temp = accum
Call fpdiv_si(fac, fac, c, digits)
Call fpmul(p, p, x9, digits)
Call fpmul(term, p, fac, digits)
Call fpadd(accum, temp, term, digits)
Loop Until fpcmp(accum, temp, digits) = 0
For i = 1 To 13
Call fpmul(accum, accum, accum, digits)
Next
result = accum
End Sub
Sub fpexp_cf (result As decfloat, x As decfloat, n As Long, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
' x
' cf = ------------------------------
' x^2
' 2 + -------------------------
' x^2
' 6 + ----------------------
' x^2
' 10 + ------------------
' x^2
' 14 + --------------
' x^2
' 18 + ---------
' x^2
' 22 + -----
' 26 + ...
' (1 + cf)
' exp(x) = ----------
' (1 - cf)
Dim i As Long
Dim As decfloat s, x1, x2, tmp, tmp2
Call fpdiv_si(x1, x, 8192, digits) ' 2^13
Call fpmul(x2, x1, x1, digits)
Call si2fp(s, 4 * n + 6, digits)
For i = n To 0 Step -1
Call si2fp(tmp, 4 * i + 2, digits)
Call fpdiv(s, x2, s, digits)
Call fpadd(s, s, tmp, digits)
Next
Call fpdiv(s, x1, s, digits)
Call si2fp(tmp, 1, digits)
tmp2 = tmp
Call fpadd(tmp, tmp, s, digits)
Call fpsub(tmp2, tmp2, s, digits)
Call fpdiv(s, tmp, tmp2, digits)
For i = 1 To 13 ' s^13
Call fpmul(s, s, s, digits)
Next
result = s
End Sub
Sub fplog_r (result As decfloat, x As decfloat, digits_in As Long)
Dim As Long digits
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
'====================Log Guard==================
Dim As decfloat x2, y, tmp
Dim As Long i
Dim As Long ex, terms, prec, n
Dim As Double t, t2
tmp.sign = 0
tmp.exponent = (BIAS + 1)
Call set(tmp.mantissa, 0, 10000000)
For i = 1 To NUM_DWORDS
Call set(tmp.mantissa, i, 0)
Next
If fpcmp(x, tmp, digits) = 0 Then result = x2: Exit Sub
ex = (x.exponent And &H7FFFFFFF) - BIAS - 1
t = git(x.mantissa, 0) + git(x.mantissa, 1) / 100000000
t = t / 10000000
t2 = Log(t) + Log(10) * ex
n = Log((digits + 1) * 0.5) * 1.5
Call str2fp(x2, Str$(t2))
prec = 3
terms = 2
prec = 2 * prec - 1
Call fpexp_cf(y, x2, terms, prec)
Call fpsub(tmp, y, x, prec)
Call fpdiv(tmp, tmp, y, prec)
Call fpsub(x2, x2, tmp, prec)
If digits > 4 And digits <= 9 Then
prec = 2 * prec - 1
terms = 2 * terms + 1
Call fpexp_cf(y, x2, terms, prec)
Call fpsub(tmp, y, x, prec)
Call fpdiv(tmp, tmp, y, prec)
Call fpsub(x2, x2, tmp, prec)
ElseIf digits > 9 And digits < 18 Then
For i = 1 To 3
prec = 2 * prec - 1
terms = 2 * terms + 1
Call fpexp_cf(y, x2, terms, prec)
Call fpsub(tmp, y, x, prec)
Call fpdiv(tmp, tmp, y, prec)
Call fpsub(x2, x2, tmp, prec)
Next
ElseIf digits > 17 And digits < 34 Then
For i = 1 To 4
prec = 2 * prec - 1
terms = 2 * terms + 1
Call fpexp_cf(y, x2, terms, prec)
Call fpsub(tmp, y, x, prec)
Call fpdiv(tmp, tmp, y, prec)
Call fpsub(x2, x2, tmp, prec)
Next
ElseIf digits > 33 And digits < 65 Then
For i = 1 To 5
prec = 2 * prec - 1
terms = 2 * terms + 1
Call fpexp_cf(y, x2, terms, prec)
Call fpsub(tmp, y, x, prec)
Call fpdiv(tmp, tmp, y, prec)
Call fpsub(x2, x2, tmp, prec)
Next
ElseIf digits > 64 Then
For i = 1 To 6
prec = 2 * prec - 1
terms = 2 * terms + 1
Call fpexp_cf(y, x2, terms, prec)
Call fpsub(tmp, y, x, prec)
Call fpdiv(tmp, tmp, y, prec)
Call fpsub(x2, x2, tmp, prec)
Next
End If
result = x2
End Sub
Sub fppow (result As decfloat, lhs As decfloat, rhs As decfloat)
Dim As decfloat lhs2
Call fplog(lhs2, lhs, NUM_DWORDS)
Call fpmul(lhs2, lhs2, rhs, NUM_DWORDS)
Call fpexp(lhs2, lhs2, NUM_DWORDS)
result = lhs2
End Sub
Sub fpnroot (result As decfloat, x As decfloat, p_in As Long, digits_in As Long)
Dim As Long digits, p, psign
digits = digits_in
If digits > NUM_DWORDS Then digits = NUM_DWORDS
Dim As decfloat ry, tmp, tmp2
Dim As Double t, t2
Dim As Long i, ex, l, prec
x.sign = 0
psign = Sgn(p_in)
p = Abs(p_in)
l = Log((NUM_DIGITS + 8) * 0.0625) * 1.5 'calculate the number of iterations needed
ex = (x.exponent And &H7FFFFFFF) - BIAS - 1 'extract the exponent
t = git(x.mantissa, 0) + git(x.mantissa, 1) / 100000000 'get 16 digits from x.mantissa
'for example, if x = 3.1415926535897932384626433832795028842 then the above would give
't = 31415926.53589793 because the mantissa doesn't have a decimal point, it's an integer
'each element of the mantissa holds 8 digits
'in this example ex = 0
t = t / 10000000 'now t = 3.141592653589793
t2 = Log(t) / p '+ Log(10) * ex / p 'log(t ^ (1/p)) = log(t)/p + Log(10) * ex / p
'in this example since ex = 0, it becomes: log(t ^ (1/p)) = log(t)/p
t2 = Exp(t2) 't2=t ^ (1/p)
Call str2fp(ry, Str$(t2)) 'convert the double t2 to decfloat in ry
t = Log(10) * ex / p
t2 = Exp(t - Fix(t))
Call str2fp(tmp, Str$(t2)) 'convert the double t2 to decfloat in tmp
Call fpmul(ry, ry, tmp, prec) 'ry = ry * Log(10) * ex / p
str2fp tmp, "2.7182818284590452353602874713527"
fpipow tmp, tmp, Fix(t), 24
Call fpmul(ry, ry, tmp, 24)
prec = 3 '3 * 8 = 24 digits, prec here means number of 8 digit elements
Call fpipow(tmp, ry, p - 1, prec) 'tmp = ry ^ (p-1)
Call fpdiv(tmp, x, tmp, prec) 'tmp = x * tmp
Call fpmul_si(tmp2, ry, p - 1, prec) 'tmp2 = ry * (p-1)
Call fpadd(tmp2, tmp2, tmp, prec) 'tmp2 = tmp2 + tmp
Call fpdiv_si(ry, tmp2, p, prec) 'ry = tmp2 / p
For i = 1 To l + 1
prec = 2 * prec - 1
Call fpipow(tmp, ry, p - 1, prec) 'tmp = ry^(p-1)
Call fpdiv(tmp, x, tmp, prec) 'tmp = x/tmp
Call fpmul_si(tmp2, ry, p - 1, prec) 'tmp2 = ry*(p-1)
Call fpadd(tmp2, tmp2, tmp, prec) 'tmp2 = tmp2+tmp
Call fpdiv_si(ry, tmp2, p, prec) 'ry = tmp2/p
Next
If psign < 0 Then
si2fp tmp, 1, digits
fpdiv ry, tmp, ry, digits
End If
result = ry
End Sub
Sub pi_brent_salamin (pi_bs As decfloat, digits_in As Unsigned Long)
Dim As Unsigned Long digits
digits = digits_in
If digits > NUMBER_OF_DIGITS Then digits = NUMBER_OF_DIGITS
Dim As Long limit2
Dim As decfloat c0, c1, c2, c05
Dim As decfloat a, b, sum
Dim As decfloat ak, bk, ck
Dim As decfloat ab, asq
Dim As decfloat pow2, tmp
limit2 = -digits + BIAS + 1
Call si2fp(c0, 0, NUMBER_OF_DIGITS): ak = c0: bk = c0: ab = c0: asq = c0
Call si2fp(c1, 1, NUMBER_OF_DIGITS): a = c1: ck = c1: pow2 = c1
Call si2fp(c2, 2, NUMBER_OF_DIGITS): b = c2
Call str2fp(c05, ".5"): sum = c05
Call si2fp(pi_bs, 3, NUMBER_OF_DIGITS)
Call fpsqr(b, b, NUMBER_OF_DIGITS)
Call fpdiv(b, c1, b, NUMBER_OF_DIGITS)
While fpcmp(ck, c0, NUMBER_OF_DIGITS) <> 0 And ck.exponent > limit2
Call fpadd(ak, a, b, NUMBER_OF_DIGITS)
Call fpmul(ak, c05, ak, NUMBER_OF_DIGITS)
Call fpmul(ab, a, b, NUMBER_OF_DIGITS)
Call fpsqr(bk, ab, NUMBER_OF_DIGITS)
Call fpmul(asq, ak, ak, NUMBER_OF_DIGITS)
Call fpsub(ck, asq, ab, NUMBER_OF_DIGITS)
Call fpmul(pow2, pow2, c2, NUMBER_OF_DIGITS)
Call fpmul(tmp, pow2, ck, NUMBER_OF_DIGITS)
Call fpsub(sum, sum, tmp, NUMBER_OF_DIGITS)
a = ak: b = bk
Wend
Call fpdiv(tmp, asq, sum, NUMBER_OF_DIGITS)
Call fpmul(pi_bs, c2, tmp, NUMBER_OF_DIGITS)
End Sub
Sub fpsin (result As decfloat, x As decfloat, digits_in As Unsigned Long)
Dim As decfloat XX, Term, Accum, p, temp2, fac, x_2
Dim As decfloat pi2, circ, Ab
Dim As Long precision
precision = digits_in
x_2 = x
pi2 = pi_dec
Call fpmul_si(circ, pi2, 2, precision)
Call fpabs(Ab, x)
If fpcmp(Ab, circ, precision) > 0 Then
'======== CENTRALIZE ==============
'floor/ceil to centralize
Dim As decfloat tmp, tmp2
If precision > 20 Then
pi2 = pi_dec
End If
Call fpmul_si(pi2, pi2, 2, precision) 'got 2*pi
Call fpdiv(tmp, x_2, pi2, precision)
tmp2 = tmp
Call fpfix(tmp, tmp) 'int part
Call fpsub(tmp, tmp2, tmp, precision) 'frac part
Call fpmul(tmp, tmp, pi2, precision)
x_2 = tmp
End If
Dim As Long lm, limit2, i
Dim As decfloat factor
lm = precision
limit2 = Int(-0.45344993886092585968 + 0.022333002852398072433 * lm + 5.0461814408333079844E-7 * lm * lm - 4.2338453039804235772E-11 * lm * lm * lm)
Call si2fp(factor, 5, precision)
Call fpipow(factor, factor, limit2, precision)
Call fpdiv(x_2, x_2, factor, precision) 'x_=x_/5^limit2
'==================================
Dim As Long sign(3): sign(3) = 1
Dim As Long c: c = 1
Accum = x_2
Call si2fp(fac, 1, precision)
p = x_2
Call fpmul(XX, x_2, x_2, precision)
Do
c = c + 2
temp2 = Accum
Call fpmul_si(fac, fac, c * (c - 1), precision)
Call fpmul(p, p, XX, precision)
Call fpdiv(Term, p, fac, precision)
If sign(c And 3) Then
Call fpsub(Accum, temp2, Term, precision)
Else
Call fpadd(Accum, temp2, Term, precision)
End If
Loop Until fpcmp(Accum, temp2, precision) = 0
'multiply the result by 5^limit2
For i = 1 To limit2
Call fpmul(p, Accum, Accum, precision)
Call fpmul(temp2, Accum, p, precision)
'*** sin(5*x) = 5 * sin(x) - 20 * sin(x)^3 + 16 * sin(x)^5
Call fpmul_si(Accum, Accum, 5, precision)
Call fpmul_si(Term, temp2, 20, precision)
Call fpmul_si(XX, temp2, 16, precision)
Call fpmul(XX, XX, p, precision)
Call fpsub(Accum, Accum, Term, precision)
Call fpadd(Accum, Accum, XX, precision)
Next i
result = Accum
End Sub
Sub fpcos (result As decfloat, z As decfloat, digits_in As Unsigned Long)
Dim As decfloat x_2, pi2
Dim As Long precision
precision = digits_in
Call fpdiv_si(pi2, pi_dec, 2, precision)
Call fpsub(x_2, pi2, z, precision)
Call fpsin(result, x_2, precision)
End Sub
Sub fptan (result As decfloat, z As decfloat, digits_in As Unsigned Long)
Dim As decfloat x_2, s, c
Dim As Long precision
precision = digits_in
x_2 = z
Call fpsin(s, x_2, precision)
x_2 = z
Call fpcos(c, x_2, precision)
Call fpdiv(result, s, c, precision)
End Sub
Sub fpatn (result As decfloat, x As decfloat, digits_in As Unsigned Long)
Dim As Long precision, z
precision = digits_in
Dim As Long sign(3): sign(3) = 1
Dim As Unsigned Long c: c = 1
Dim As decfloat XX, Term, Accum, strC, x_2, mt, mt2, p
Dim As decfloat decnum, one, decnum2, factor
decnum2 = x
decnum2.sign = 0
Call si2fp(one, 1, precision)
If fpcmp(decnum2, one, precision) = 0 Then
Call fpdiv_si(result, pi_dec, 4, precision)
result.sign = x.sign
Exit Sub
End If
decnum2.sign = x.sign
Dim As Long limit2: limit2 = 16
Call si2fp(factor, ShL(2, limit2 - 1), precision)
For z = 1 To limit2
Call fpmul(decnum, decnum2, decnum2, precision)
Call fpadd(decnum, decnum, one, precision)
Call fpsqr(decnum, decnum, precision)
Call fpadd(decnum, decnum, one, precision)
Call fpdiv(decnum, decnum2, decnum, precision)
decnum2 = decnum
Next z
mt = decnum
x_2 = decnum
p = decnum
Call fpmul(XX, x_2, x_2, precision)
Do
c = c + 2
mt2 = mt
Call si2fp(strC, c, precision)
Call fpmul(p, p, XX, precision)
Call fpdiv(Term, p, strC, precision)
If sign(c And 3) Then
Call fpsub(Accum, mt, Term, precision)
Else
Call fpadd(Accum, mt, Term, precision)
End If
Swap mt, Accum
Loop Until fpcmp(mt, mt2, precision) = 0
Call fpmul(result, factor, mt, precision)
End Sub
Sub fpasin (result As decfloat, x As decfloat, digits_in As Unsigned Long)
Dim As Long precision
precision = digits_in
Dim As Double num
' ARCSIN = ATN(x / SQR(-x * x + 1))
'============= ARCSIN GUARD =========
num = fp2dbl(x)
If num > 1 Then Exit Sub
If num < -1 Then Exit Sub
'========================
Dim As decfloat one, T, B, term1, minusone
Call si2fp(one, 1, precision)
Call si2fp(minusone, -1, precision)
T = x
Call fpmul(B, x, x, precision) 'x*x
'for 1 and -1
If fpcmp(B, one, precision) = 0 Then
Dim As decfloat two, atn1
Call si2fp(two, 2, precision)
Call fpatn(atn1, one, precision)
If fpcmp(x, minusone, precision) = 0 Then
Call fpmul(two, two, atn1, precision)
Call fpmul(result, two, minusone, precision)
Exit Sub
Else
Call fpmul(result, two, atn1, precision)
Exit Sub
End If
End If
Call fpsub(B, one, B, precision) '1-x*x
Call fpsqr(B, B, precision) 'sqr(1-x*x)
Call fpdiv(term1, T, B, precision)
Call fpatn(result, term1, precision)
End Sub
Sub fpacos (result As decfloat, x As decfloat, digits_in As Unsigned Long)
Dim As Long precision
precision = digits_in
Dim As decfloat one, minusone, two, atn1, tail, T, B, term1, atnterm1 ',_x,temp
Dim As Double num
'ARCCOS = ATN(-x / SQR(-x * x + 1)) + 2 * ATN(1)
'============= ARCCOS GUARD =========
num = fp2dbl(x)
If num > 1 Then Exit Sub
If num < -1 Then Exit Sub
'========================
Call si2fp(one, 1, precision)
Call si2fp(minusone, -1, precision)
Call si2fp(two, 2, precision)
Call fpatn(atn1, one, precision)
Call fpmul(tail, two, atn1, precision) '2*atn(1)
Call fpmul(T, minusone, x, precision) '-x
Call fpmul(B, x, x, precision) 'x*x
If fpcmp(B, one, precision) = 0 Then
'for 1 and -1
If fpcmp(x, minusone, precision) = 0 Then
Dim As decfloat four
Call si2fp(four, 4, precision)
Call fpmul(result, four, atn1, precision)
Exit Sub
Else
Call si2fp(result, 0, precision)
Exit Sub
End If
End If
Call fpsub(B, one, B, precision) '1-x*x
Call fpsqr(B, B, precision) 'sqr(1-x*x)
Call fpdiv(term1, T, B, precision)
Call fpatn(atnterm1, term1, precision)
Call fpadd(result, atnterm1, tail, precision)
End Sub
Sub initialize_fp
Print "initializing constants: please wait"
Call pi_brent_salamin(pi_dec, NUMBER_OF_DIGITS)
Call si2fp(tan_half_num(0), 992, NUMBER_OF_DIGITS)
Call str2fp(tan_half_num(1), "161388480")
Call str2fp(tan_half_num(2), "7686610525440")
Call str2fp(tan_half_num(3), "-167256984742848000")
Call str2fp(tan_half_num(4), "2000393537524462080000")
Call str2fp(tan_half_num(5), "-14467646220791919547392000")
Call str2fp(tan_half_num(6), "66677300813465118447390720000")
Call str2fp(tan_half_num(7), "-201117789910786072985458237440000")
Call str2fp(tan_half_num(8), "400342706504764747636935691468800000")
Call str2fp(tan_half_num(9), "-521967288977995909272835160309760000000")
Call str2fp(tan_half_num(10), "435052278687602761865918494187323392000000")
Call str2fp(tan_half_num(11), "-221463964316902607512694240578598338560000000")
Call str2fp(tan_half_num(12), "63663507608965602906315837691661069058048000000")
Call str2fp(tan_half_num(13), "-8994510946805140308046160658488525397688320000000")
Call str2fp(tan_half_num(14), "470550277118574335327341015729793791741132800000000")
Call str2fp(tan_half_num(15), "-3827142253897737927329040261268989506161213440000000")
Call si2fp(tan_half_den(0), 491040, NUMBER_OF_DIGITS)
Call str2fp(tan_half_den(1), "39540177600")
Call str2fp(tan_half_den(2), "-1232419887578880")
Call str2fp(tan_half_den(3), "19569067214913216000")
Call str2fp(tan_half_den(4), "-180435497084706479616000")
Call str2fp(tan_half_den(5), "1036847979156754234229760000")
Call str2fp(tan_half_den(6), "-3857758118493338995884748800000")
Call str2fp(tan_half_den(7), "9452536125806945430316537159680000")
Call str2fp(tan_half_den(8), "-15257505370126034271052104685977600000")
Call str2fp(tan_half_den(9), "15972199042726674823748755905478656000000")
Call str2fp(tan_half_den(10), "-10480804895655884717678945541785518080000000")
Call str2fp(tan_half_den(11), "4060172679143214471066061077274302873600000000")
Call str2fp(tan_half_den(12), "-837419984702547545921539095790310985302016000000")
Call str2fp(tan_half_den(13), "75810877980214754024960496978688999780515840000000")
Call str2fp(tan_half_den(14), "-1913571126948868963664520130634494753080606720000000")
Cls
End SubRE: decfloat -- again - Jack - 09-13-2022
sample run
Code: (Select All)
(Pi * 10^2000) ^ (1/3)
6798033351105428972796748712399168896288062921435309418105683598078476864738298960827006729808515426265943084793580899570645182344095301929614027410723460129225988067511365977568705516838730926288492233363900259587955916964663852003454674744146788305116431034206795837640371505764617850080034461533129170778315679586851291320339973886105059270306455700926462667888048241388351931217312091552532072744334169572736176567701342465982087666609795187459304748316788025109585867953771470483306863071529154237239012481213716924794320232073073207197866771876365643607933611217695182503432360129741097941230465215353461448280987779091219043128991891321941453822525251856646896.5162217573545876925479621003234982405877425169332532134776380519701920896262086090634232587956995974174358710237357682879063357708089640366285593961817396810578301066266996880451474299616453161076725390163822554703142811524839508435133033396737457622596770653810396626157075556180227381227490588758212715009553628957698357573793013178292902984345684802531630351754142944175206160908513237350232013216980072591690824614855337946089140129092639878670786691860592315904254319887961195947977404455314864350648172442459221386785069930645535320152617455911747093215902405426698189236047769996088564520865324838768138514373321804697588101054138539792705783350114868020073573890641634263186524746778983139506101446912304993725079604388965143093210115906759913920317434802641280484937166169955427962151488638366363774632645661667087745350133370837260829277370291220032608830931847940817480098205283930562277564189420954300058096820688252682197752401285069691156372366569517089444322697778977462467249743716611187150659816884398200059645223334944496859169139011306257908882950326502168199674231375682727946092815765848694903837220377568464422525690223793221005170218652708523128614184339214197304809820866749104314399394267004331659545190261752743714384285714595286263226031941293481840586569284991380553181517304082907209088957837406337916458
elapsed time .2900000000008731 seconds
enter the number digits for output (max 2000) ? 200
Please enter a number for root extraction, an empty input will exit
enter number ? 2
enter an Integer root# ? 100
1.0069555500567188088326982141132397854535407405341259051168656964887105974087030025255827827631333722908158160857148534344012233193291411228058249641443268336500649821316291594690183759920154612534944
elapsed time .150999999998021 secondsRE: decfloat -- again - Jack - 09-14-2022
the other day I noticed in the QB64 wiki that a _Float takes 32 bytes and that made me think that a rewrite of decfloat might be a good idea, I was thinking of using a 32 byte structure of 8 longs, 1 for the exponent and 7 for the mantissa at 8 digits per long it would give 56 digits, any thoughts ?
the structure would look like this
Code: (Select All)
Type decfloat
As Unsigned Long exponent
As Long M0
As Long M1
As Long M2
As Long M3
As Long M4
As Long M5
As Long M6
End TypeRE: decfloat -- again - Pete - 09-14-2022
I understood everything in the first post up to and including the line written as:
Code: (Select All)
_Title "defloat by Jack"Actually, I think I could learn a lot from that post, if I can take the time to pick it apart. Very interested in the trig stuff.
Jack, did you ever wonder why none of the earlier developers did not include a way to connect to the C math lib? I mean QB64 is a c/c++ translator, after all. Maybe I'm missing something, but when you ask about building more precise math functions at a C forum, the first reply you get is to include the math lib.
Anyway, are you suggesting to make _FLOAT from a 32 bit to a 64 bit type variable designation? That was my take on your last reply.
Pete
RE: decfloat -- again - SpriggsySpriggs - 09-14-2022
So, what would be the use case for a number with this many decimal points? Would you even be able to reliably perform math on this number?
RE: decfloat -- again - Jack - 09-14-2022
I think lack of sleep led me to make a stupid post, please ignore this.
RE: decfloat -- again - Pete - 09-14-2022
There are no such things as stupid posts... only stupid threads. Oops, that didn't help at all.
I hear you on the lack of sleep issue, but my wife says that's when I make my best decisions.... which means she waits until I fall asleep, and then goes and buys crap online. Oh well, who knew all 256 shades of pink lipstick are all considered essential?
Pete
RE: decfloat -- again - SpriggsySpriggs - 09-14-2022
(09-14-2022, 06:26 PM)Jack Wrote: I think lack of sleep led me to make a stupid post, please ignore this.
@Jack
I don't think it's stupid. I just don't know how to use it or when.
RE: decfloat -- again - BSpinoza - 09-15-2022
It could be a preparation for the Pi-Day: https://www.piday.org/
RE: decfloat -- again - Pete - 09-15-2022
Beats "Pride" month. Rainbows are just half-ascii circles.
Pete