Posts: 1,000
Threads: 50
Joined: May 2022
Reputation:
27
09-20-2022, 09:54 PM
(This post was last modified: 09-20-2022, 09:55 PM by Kernelpanic.)
I think it fits here. I bought this book: Numeric for fans for €13.44. It's better than described. A library copy, outside OK and all pages clean.
A dream for math fans!  Unfortunately, I'm not a mathematician, and the book is too much for me. Well, I understand quite a few things and I'll try it out. I also have all the programs for this, the author sent me after I asked.
So I can only recommend this book to anyone who is familiar with numerics.
Posts: 422
Threads: 27
Joined: Apr 2022
Reputation:
26
proof of concept of a decimal floating point type using 32-bytes
Code: (Select All) _Title "defloat-32 by Jack"
$NoPrefix
$Console:Only
Dest Console
Dim Shared As Long shift_constants(7)
shift_constants(0) = 10
shift_constants(1) = 100
shift_constants(2) = 1000
shift_constants(3) = 10000
shift_constants(4) = 100000
shift_constants(5) = 1000000
shift_constants(6) = 10000000
shift_constants(7) = 100000000
Const BIAS = 16384 '2 ^ 14
Const DWORDS = 7
Type decfloat
As _Unsigned Integer exponent
As Integer sign
As _Unsigned Long M0
As _Unsigned Long M1
As _Unsigned Long M2
As _Unsigned Long M3
As _Unsigned Long M4
As _Unsigned Long M5
As _Unsigned Long M6
End Type
Dim Shared As decfloat pi_dec, pi2_dec, pi_dec_half
str2fp pi_dec, "3.14159265358979323846264338327950288419716939937510582097494459"
str2fp pi2_dec, "6.28318530717958647692528676655900576839433879875021164194988918"
str2fp pi_dec_half, "1.57079632679489661923132169163975144209858469968755291048747230"
' Error definitions
Const DIVZ_ERR = 1 'Divide by zero
Const EXPO_ERR = 2 'Exponent overflow error
Const EXPU_ERR = 3 'Exponent underflow error
Dim As decfloat x, y, z
Dim As Long i
'Dim As Double t
For i = 1 To 20
Si2fp x, i
fpsqr y, x
fpmul z, y, y
fpsub z, z, x
fpdiv z, z, x
Print i, fp2str(y, 58); " rel. error "; fp2str(z, 10)
Next
Sub fpatn (result As decfloat, x As decfloat)
Dim As Long sign(3): sign(3) = 1
Dim As Long z, 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
one.sign = 0
one.exponent = (BIAS + 1)
one.M0 = 100000000
one.M1 = 0
one.M2 = 0
one.M3 = 0
one.M4 = 0
one.M5 = 0
one.M6 = 0
If fpcmp(decnum2, one) = 0 Then
fpdiv_si decnum, pi_dec, 4
decnum.sign = x.sign
result = decnum
Exit Sub
End If
decnum2.sign = x.sign
Dim As Long lm2: lm2 = 16
Si2fp factor, _ShL(2, (lm2 - 1))
For z = 1 To lm2
fpmul decnum, decnum2, decnum2
fpadd decnum, decnum, one
fpsqr decnum, decnum
fpadd decnum, decnum, one
fpdiv decnum, decnum2, decnum
decnum2 = decnum
Next z
mt = decnum
x_2 = decnum
p = decnum
fpmul XX, x_2, x_2
Do
c = c + 2
mt2 = mt
Si2fp strC, c
fpmul p, p, XX
fpdiv Term, p, strC
If sign(c And 3) Then
fpsub Accum, mt, Term
Else
fpadd Accum, mt, Term
End If
Swap mt, Accum
Loop Until fpcmp(mt, mt2) = 0
fpmul result, factor, mt
End Sub
Sub fpasin (result As decfloat, x As decfloat)
Dim As Double num
Dim As decfloat one, T, B, term1, minusone
' ARCSIN = ATN(x / SQR(-x * x + 1))
'============= ARCSIN GUARD =========
num = Val(fp2str(x, 16))
If num > 1 Then
result = one
Exit Sub
End If
If num < -1 Then
result = one
Exit Sub
End If
'========================
one.sign = 0
one.exponent = (BIAS + 1)
one.M0 = 100000000
one.M1 = 0
one.M2 = 0
one.M3 = 0
one.M4 = 0
one.M5 = 0
one.M6 = 0
minusone = one: minusone.sign = &H8000
T = x
fpmul B, x, x 'x*x
'for 1 and -1
If fpcmp(B, one) = 0 Then
Dim As decfloat two, atn1
two.sign = 0
two.exponent = (BIAS + 1)
two.M0 = 200000000
two.M0 = 0
two.M0 = 0
two.M0 = 0
two.M0 = 0
two.M0 = 0
two.M0 = 0
fpdiv_si atn1, pi_dec, 4
If fpcmp(x, minusone) = 0 Then
fpmul two, two, atn1
fpmul result, two, minusone
Exit Sub
Else
result = pi_dec_half
Exit Sub
End If
End If
fpsub B, one, B '1-x*x
fpsqr B, B 'sqr(1-x*x)
fpdiv term1, T, B
fpatn result, term1
End Sub
Sub fpcos (result As decfloat, z As decfloat)
Dim As decfloat x_2
fpsub x_2, pi_dec_half, z
fpsin result, x_2
End Sub
Sub fptan (result As decfloat, z As decfloat)
Dim As decfloat x_2, s, c
x_2 = z
fpsin s, x_2
x_2 = z
fpcos c, x_2
fpdiv result, s, c
End Sub
Sub fpsin (result As decfloat, x As decfloat)
Dim As decfloat XX, Term, Accum, p, temp2, fac, x_2
Dim As decfloat pi2, circ, Ab
x_2 = x
pi2 = pi_dec
circ = pi2_dec
fpabs Ab, x
If fpcmp(Ab, circ) > 0 Then
'======== CENTRALIZE ==============
'floor/ceil to centralize
Dim As decfloat tmp, tmp2
pi2 = pi2_dec 'got 2*pi
fpdiv tmp, x_2, pi2
tmp2 = tmp
fpfix tmp, tmp 'int part
fpsub tmp, tmp2, tmp 'frac part
fpmul tmp, tmp, pi2
x_2 = tmp
End If
Dim As Long lm, lm2, i
Dim As decfloat factor
lm = 63
lm2 = 1 + Int(-0.45344993886092585968 + 0.022333002852398072433 * lm + 5.0461814408333079844E-7 * lm * lm - 4.2338453039804235772E-11 * lm * lm * lm)
If lm2 < 0 Then lm2 = 0
Si2fp factor, 5
fpipow factor, factor, lm2
fpdiv x_2, x_2, factor 'x_=x_/5^lm2
'==================================
Dim As Long sign(3): sign(3) = 1
Dim As Long c: c = 1
Accum = x_2
Si2fp fac, 1
p = x_2
fpmul XX, x_2, x_2
Do
c = c + 2
temp2 = Accum
fpmul_si fac, fac, c * (c - 1)
fpmul p, p, XX
fpdiv Term, p, fac
If sign(c And 3) Then
fpsub Accum, temp2, Term
Else
fpadd Accum, temp2, Term
End If
Loop Until fpcmp(Accum, temp2) = 0
'multiply the result by 5^lm2
For i = 1 To lm2
fpmul p, Accum, Accum
fpmul temp2, Accum, p
'*** sin(5*x) = 5 * sin(x) - 20 * sin(x)^3 + 16 * sin(x)^5
fpmul_si Accum, Accum, 5
fpmul_si Term, temp2, 20
fpmul_si XX, temp2, 16
fpmul XX, XX, p
fpsub Accum, Accum, Term
fpadd Accum, Accum, XX
Next i
result = Accum
End Sub
Sub fppow (result As decfloat, lhs As decfloat, rhs As decfloat)
Dim As decfloat lhs2
fplog lhs2, lhs
fpmul lhs2, lhs2, rhs
fpexp result, lhs2
End Sub
Sub fpexp (result As decfloat, x As decfloat)
'taylor series
Dim As decfloat fac, ux, temp, accum, p, term
Dim As Long i, c
temp.sign = 0
temp.exponent = 0
temp.M0 = 0
temp.M1 = 0
temp.M2 = 0
temp.M3 = 0
temp.M4 = 0
temp.M5 = 0
temp.M6 = 0
fac.sign = 0
fac.exponent = (BIAS + 1)
fac.M0 = 100000000
fac.M1 = 0
fac.M2 = 0
fac.M3 = 0
fac.M4 = 0
fac.M5 = 0
fac.M6 = 0
If fpcmp(x, temp) = 0 Then
result = fac
Exit Sub
End If
c = 1
fpdiv_si ux, x, 512
p = ux
fpadd accum, fac, ux '1 + x
Do
c = c + 1
temp = accum
fpdiv_si fac, fac, c
fpmul p, p, ux
fpmul term, p, fac
fpadd accum, temp, term
Loop Until fpcmp(accum, temp) = 0
For i = 1 To 9
fpmul accum, accum, accum
Next
result = accum
End Sub
Sub fplogTaylor (result As decfloat, x As decfloat)
'taylor series
'====================Log Guard==================
Dim As decfloat g, zero
fpabs g, x
zero.sign = 0
zero.exponent = 0
zero.M0 = 0
zero.M1 = 0
zero.M2 = 0
zero.M3 = 0
zero.M4 = 0
zero.M5 = 0
zero.M6 = 0
If fpcmp(g, x) <> 0 Then
result = zero
Exit Sub
End If
If fpcmp(x, zero) = 0 Then
result = zero
Exit Sub
End If
'=============================================
Dim As Long invflag
Dim As decfloat XX, Term, Accum, ux, tmp, tmp2
Dim As decfloat T, B, one, Q, two
one.sign = 0
one.exponent = (BIAS + 1)
one.M0 = 100000000
two.sign = 0
two.exponent = (BIAS + 1)
two.M0 = 200000000
one.M1 = 0: two.M0 = 0
one.M2 = 0: two.M0 = 0
one.M3 = 0: two.M0 = 0
one.M4 = 0: two.M0 = 0
one.M5 = 0: two.M0 = 0
one.M6 = 0: two.M0 = 0
ux = x
If fpcmp(x, one) < 0 Then
invflag = 1
fpdiv ux, one, ux
End If
fpsub T, ux, one
fpadd B, ux, one
fpdiv Accum, T, B
fpdiv Q, T, B
tmp = Q
fpmul XX, Q, Q
Dim As Long c: c = 1
Do
c = c + 2
tmp2 = tmp
fpmul Q, Q, XX
fpdiv_si Term, Q, c
fpadd Accum, tmp, Term
Swap tmp, Accum
Loop Until fpcmp(tmp, tmp2) = 0
fpmul_si Accum, Accum, 2
If invflag Then
fpneg Accum, Accum
End If
result = Accum
End Sub
Sub fplog (result As decfloat, x As decfloat)
'====================Log Guard==================
Dim As decfloat g, one, zero
Dim As Long factor
zero.sign = 0
zero.exponent = 0
zero.M0 = 0
zero.M1 = 0
zero.M2 = 0
zero.M3 = 0
zero.M4 = 0
zero.M5 = 0
zero.M6 = 0
one.sign = 0
one.exponent = (BIAS + 1)
one.M0 = 100000000
one.M1 = 0
one.M2 = 0
one.M3 = 0
one.M4 = 0
one.M5 = 0
one.M6 = 0
fpabs g, x
If fpcmp(g, x) <> 0 Then
result = zero
Exit Sub
End If
If fpcmp(x, zero) = 0 Then
result = zero
Exit Sub
End If
If fpcmp(x, one) = 0 Then
result = zero
Exit Sub
End If
'=============================================
Dim As decfloat approx, ans, logx
approx = x
factor = 4096
fpsqr approx, approx
fpsqr approx, approx
fpsqr approx, approx
fpsqr approx, approx
fpsqr approx, approx
fpsqr approx, approx
fpsqr approx, approx
fpsqr approx, approx
fpsqr approx, approx
fpsqr approx, approx
fpsqr approx, approx
fpsqr approx, approx
fplogTaylor logx, approx
fpmul_si ans, logx, factor
result = ans
End Sub
' sqrt(num)
Sub fpsqr (result As decfloat, num As decfloat)
Dim As decfloat r, r2, tmp, n, half
Dim As Integer ex, k
Dim As String ts, v
Dim As Double x
'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
Si2fp tmp, 0
If fpcmp&(n, tmp) = 0 Then
Si2fp r, 0
result = r
Exit Sub
End If
Si2fp tmp, 1
If fpcmp&(n, tmp) = 0 Then
Si2fp r, 1
result = r
Exit Sub
End If
Si2fp tmp, 0
If fpcmp&(n, tmp) < 0 Then
Si2fp r, 0
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 &H7FFF) - BIAS - 1
Else
ex = 0
End If
ts = _Trim$(Str$(n.M0))
If Len(ts) < 9 Then
ts = ts + String$(9 - Len(ts), "0")
End If
v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
ts = _Trim$(Str$(n.M1))
If Len(ts) < 9 Then
ts = String$(9 - 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
Si2fp r, 1
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, ".")
str2fp r, v
r.exponent = ex \ 2 + BIAS + 1
If Len(v) > 1 And k = 0 Then r.exponent = r.exponent + 1
'half.M0=500000000
'half.exponent = BIAS
'half.sign = 0
str2fp half, ".5"
'=====================================================================
'Newton-Raphson method
For k = 1 To 2
fpdiv tmp, n, r
fpadd r2, r, tmp
fpmul r, r2, half
Next
result = r
End Sub
Sub fpnroot (result As decfloat, x As decfloat, p_in As Long)
Dim As Long p, psign
Dim As decfloat ry, tmp, tmp2
Dim As Double t, t2
Dim As Long i, ex, l
x.sign = 0
psign = Sgn(p_in)
p = Abs(p_in)
l = 2 'Log((NUM_DIGITS + 8) * 0.0625) * 1.5 'calculate the number of iterations needed
ex = (x.exponent And &H7FFF) - BIAS - 1 'extract the exponent
t = x.M0 + x.M1 / 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 / 100000000 '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)
str2fp ry, Str$(t2) 'convert the double t2 to decfloat in ry
t = Log(10) * ex / p
t2 = Exp(t - Fix(t))
str2fp tmp, Str$(t2) 'convert the double t2 to decfloat in tmp
fpmul ry, ry, tmp 'ry = ry * Log(10) * ex / p
str2fp tmp, "2.7182818284590452353602874713527"
fpipow tmp, tmp, Fix(t)
fpmul ry, ry, tmp
fpipow tmp, ry, p - 1 'tmp = ry ^ (p-1)
fpdiv tmp, x, tmp 'tmp = x * tmp
fpmul_si tmp2, ry, p - 1 'tmp2 = ry * (p-1)
fpadd tmp2, tmp2, tmp 'tmp2 = tmp2 + tmp
fpdiv_si ry, tmp2, p 'ry = tmp2 / p
For i = 1 To l + 1
fpipow tmp, ry, p - 1 'tmp = ry^(p-1)
fpdiv tmp, x, tmp 'tmp = x/tmp
fpmul_si tmp2, ry, p - 1 'tmp2 = ry*(p-1)
fpadd tmp2, tmp2, tmp 'tmp2 = tmp2+tmp
fpdiv_si ry, tmp2, p 'ry = tmp2/p
Next
If psign < 0 Then
Si2fp tmp, 1
fpdiv ry, tmp, ry
End If
result = ry
End Sub
Sub fpipow (result As decfloat, x As decfloat, e As _Integer64)
'take x to an Long power
Dim As decfloat y
Dim As decfloat z
Dim As _Integer64 n, c
c = 0
y = x
n = Abs(e)
z.sign = 0
z.exponent = (BIAS + 1)
z.M0 = 100000000
While n > 0
While (n And 1) = 0
n = n \ 2
fpmul y, y, y
c = c + 1
Wend
n = n - 1
fpmul z, y, z
c = c + 1
Wend
If e < 0 Then
Si2fp y, 1
fpdiv z, y, z
End If
result = z
End Sub
Sub fpneg (result As decfloat, x As decfloat)
result = x
result.sign = result.sign Xor &H8000
End Sub
Sub fpabs (result As decfloat, x As decfloat)
result = x
result.sign = 0
End Sub
Function fp2str$ (num As decfloat, digits As Integer)
Dim As decfloat n, one
Dim As Long ex, ln
Dim As String s, sd, sign
Dim As _Unsigned Integer xpn
Dim As Integer signn
n = num
xpn = n.exponent
signn = n.sign
'round up
If ndigit(n, digits + 1) > 4 Then
If digits < 1 Then digits = 1
If digits > 58 Then digits = 58
n.exponent = digits + BIAS
n.sign = 0
one.M0 = 100000000
one.exponent = 1 + BIAS
fpadd n, n, one
If n.exponent > (digits + BIAS) Then
xpn = xpn + (n.exponent - (digits + BIAS))
End If
End If
n.exponent = xpn
n.sign = signn
If n.exponent > 0 Then
ex = (n.exponent And &H7FFF) - BIAS - 1
Else
If n.exponent = 0 Then
fp2str$ = " 0"
Exit Function
Else
fp2str$ = "Exponent overflow"
Exit Function
End If
End If
If n.sign Then sign = "-" Else sign = " "
sd = ""
s = _Trim$(Str$(n.M0))
sd = sd + s
s = _Trim$(Str$(n.M1))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
s = _Trim$(Str$(n.M2))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
s = _Trim$(Str$(n.M3))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
s = _Trim$(Str$(n.M4))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
s = _Trim$(Str$(n.M5))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
s = _Trim$(Str$(n.M6))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
sd = Left$(sd, digits)
ln = Len(sd)
If ex >= 0 Then
If ex < digits Then
sd = Left$(sd, ex + 1) + "." + Mid$(sd, ex + 2)
ElseIf ex > digits Then
sd = Left$(sd, 1) + "." + Mid$(sd, 2) + "e" + _Trim$(Str$(ex))
End If
ElseIf ex < 0 Then
If ex > (-5) Then
sd = "." + String$(Abs(ex) - 1, "0") + sd
Else
sd = Left$(sd, 1) + "." + Mid$(sd, 2) + "e" + _Trim$(Str$(ex))
End If
End If
fp2str$ = sign + sd
End Function
Function ndigit& (num As decfloat, digit As Long)
Dim As Long ex, ex2, j, k
Dim As _Unsigned Integer xp
xp = num.exponent
num.exponent = (digit + BIAS)
ex = (num.exponent And &H7FFF) - BIAS
Dim As _Unsigned Long numa(6)
numa(0) = num.M0
numa(1) = num.M1
numa(2) = num.M2
numa(3) = num.M3
numa(4) = num.M4
numa(5) = num.M5
numa(6) = num.M6
ex = ex - 1
ex2 = ex \ 9
k = ex2
j = ex Mod 9
If j = 8 Then
ndigit& = numa(k) Mod 10
Else
ndigit& = (numa(k) \ shift_constants(7 - j)) Mod 10
End If
num.exponent = xp
End Function
Function fp2str_exp$ (n As decfloat, places_in As Long)
Dim As Long ex, places
Dim As String v, f, ts
places = places_in
If places > 54 Then places = 54
places = 62 'places-1
If n.exponent > 0 Then
ex = (n.exponent And &H7FFF) - BIAS - 1
Else
ex = 0
End If
If n.sign Then
v = "-"
Else
v = " "
End If
ts = Str$(n.M0)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = ts + String$(9 - Len(ts), "0")
End If
v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
ts = Str$(n.M1)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
ts = Str$(n.M2)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
ts = Str$(n.M3)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
ts = Str$(n.M4)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
ts = Str$(n.M5)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
ts = Str$(n.M6)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
v = Left$(v, places + 3)
f = _Trim$(Str$(Abs(ex)))
'f = string$(9 - Len(f), "0") + f
If ex < 0 Then
v = v + "e-"
Else
v = v + "e+"
End If
v = v + f
fp2str_exp$ = v
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) > (63 + 1) Then
f1 = Mid$(f1, 1, (63 + 1))
End If
While Len(f1) < (63 + 1)
f1 = f1 + "0"
Wend
j = 1
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M0 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M1 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M2 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M3 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M4 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M5 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M6 = ulng
If ulng <> 0 Then fp = 1
If fp Then n.exponent = (ex + BIAS + 1) Else n.exponent = 0
result = n
End Sub
Sub Si2fp (result As decfloat, m As _Integer64)
Dim As decfloat fac1
Dim As _Integer64 n
n = Abs(m)
If n > 9999999999999999~&& Then
str2fp result, Str$(m)
Exit Sub
End If
fac1.M0 = 0
fac1.M1 = 0
fac1.M2 = 0
fac1.M3 = 0
fac1.M4 = 0
fac1.M5 = 0
fac1.M6 = 0
If m = 0 Then
fac1.exponent = 0
fac1.sign = 0
result = fac1
Exit Sub
End If
fac1.exponent = BIAS
If n < 1000000000 Then
If n < 10 Then
fac1.M0 = n * 100000000
fac1.exponent = fac1.exponent + 1
ElseIf n < 100 Then
fac1.M0 = n * 10000000
fac1.exponent = fac1.exponent + 2
ElseIf n < 1000 Then
fac1.M0 = n * 1000000
fac1.exponent = fac1.exponent + 3
ElseIf n < 10000 Then
fac1.M0 = n * 100000
fac1.exponent = fac1.exponent + 4
ElseIf n < 100000 Then
fac1.M0 = n * 10000
fac1.exponent = fac1.exponent + 5
ElseIf n < 1000000 Then
fac1.M0 = n * 1000
fac1.exponent = fac1.exponent + 6
ElseIf n < 10000000 Then
fac1.M0 = n * 100
fac1.exponent = fac1.exponent + 7
ElseIf n < 100000000 Then
fac1.M0 = n * 10
fac1.exponent = fac1.exponent + 8
ElseIf n < 1000000000 Then
fac1.M0 = n
fac1.exponent = fac1.exponent + 9
End If
End If
If n > 999999999 Then
fac1.exponent = fac1.exponent + 9
If n < 10000000000 Then
fac1.M0 = n \ 10
fac1.M1 = (n Mod 10) * 100000000
fac1.exponent = fac1.exponent + 1
ElseIf n < 100000000000 Then
fac1.M0 = n \ 100
fac1.M1 = (n Mod 100) * 10000000
fac1.exponent = fac1.exponent + 2
ElseIf n < 1000000000000 Then
fac1.M0 = n \ 1000
fac1.M1 = (n Mod 1000) * 1000000
fac1.exponent = fac1.exponent + 3
ElseIf n < 10000000000000 Then
fac1.M0 = n \ 10000
fac1.M1 = (n Mod 10000) * 100000
fac1.exponent = fac1.exponent + 4
ElseIf n < 100000000000000 Then
fac1.M0 = n \ 100000
fac1.M1 = (n Mod 100000) * 10000
fac1.exponent = fac1.exponent + 5
ElseIf n < 1000000000000000 Then
fac1.M0 = n \ 1000000
fac1.M1 = (n Mod 1000000) * 1000
fac1.exponent = fac1.exponent + 6
ElseIf n < 10000000000000000 Then
fac1.M0 = n \ 10000000
fac1.M1 = (n Mod 10000000) * 100
fac1.exponent = fac1.exponent + 7
ElseIf n < 100000000000000000 Then
fac1.M0 = n \ 100000000
fac1.M1 = (n Mod 100000000) * 10
fac1.exponent = fac1.exponent + 8
ElseIf n < 1000000000000000000 Then
fac1.M0 = n \ 1000000000
fac1.M1 = n Mod 1000000000
fac1.exponent = fac1.exponent + 9
End If
End If
If m < 0 Then
fac1.sign = &H8000
Else
fac1.sign = 0
End If
result = fac1
End Sub
Sub RSHIFT_n (mantissa As decfloat, n As Long)
If n = 9 Then
mantissa.M6 = mantissa.M5
mantissa.M5 = mantissa.M4
mantissa.M4 = mantissa.M3
mantissa.M3 = mantissa.M2
mantissa.M2 = mantissa.M1
mantissa.M1 = mantissa.M0
mantissa.M0 = 0
Exit Sub
Else
Dim As _Unsigned Long v1, v2, c1, c2
c1 = shift_constants(n - 1)
c2 = shift_constants(8 - n)
v1 = mantissa.M6 \ c1
v2 = mantissa.M5 Mod c1
v2 = v2 * c2 + v1
mantissa.M6 = v2
v1 = mantissa.M5 \ c1
v2 = mantissa.M4 Mod c1
v2 = v2 * c2 + v1
mantissa.M5 = v2
v1 = mantissa.M4 \ c1
v2 = mantissa.M3 Mod c1
v2 = v2 * c2 + v1
mantissa.M4 = v2
v1 = mantissa.M3 \ c1
v2 = mantissa.M2 Mod c1
v2 = v2 * c2 + v1
mantissa.M3 = v2
v1 = mantissa.M2 \ c1
v2 = mantissa.M1 Mod c1
v2 = v2 * c2 + v1
mantissa.M2 = v2
v1 = mantissa.M1 \ c1
v2 = mantissa.M0 Mod c1
v2 = v2 * c2 + v1
mantissa.M1 = v2
mantissa.M0 = mantissa.M0 \ c1
End If
End Sub
Sub LSHIFT_n (mantissa As decfloat, n As Long)
If n = 9 Then
mantissa.M0 = mantissa.M1
mantissa.M1 = mantissa.M2
mantissa.M2 = mantissa.M3
mantissa.M3 = mantissa.M4
mantissa.M4 = mantissa.M5
mantissa.M5 = mantissa.M6
mantissa.M6 = 0
Exit Sub
Else
Dim As _Unsigned Long v1, v2, c1, c2
c1 = shift_constants(n - 1)
c2 = shift_constants(8 - n)
v1 = mantissa.M0 Mod c2
v2 = mantissa.M1 \ c2
mantissa.M0 = v1 * c1 + v2
mantissa.M1 = mantissa.M1 Mod c2
v1 = mantissa.M1 Mod c2
v2 = mantissa.M2 \ c2
mantissa.M1 = v1 * c1 + v2
mantissa.M2 = mantissa.M2 Mod c2
v1 = mantissa.M2 Mod c2
v2 = mantissa.M3 \ c2
mantissa.M2 = v1 * c1 + v2
mantissa.M3 = mantissa.M3 Mod c2
v1 = mantissa.M3 Mod c2
v2 = mantissa.M4 \ c2
mantissa.M3 = v1 * c1 + v2
mantissa.M4 = mantissa.M4 Mod c2
v1 = mantissa.M4 Mod c2
v2 = mantissa.M5 \ c2
mantissa.M4 = v1 * c1 + v2
mantissa.M5 = mantissa.M5 Mod c2
v1 = mantissa.M5 Mod c2
v2 = mantissa.M6 \ c2
mantissa.M5 = v1 * c1 + v2
mantissa.M6 = mantissa.M6 Mod c2
mantissa.M6 = c1 * (mantissa.M6 Mod c2)
End If
End Sub
Sub LSHIFT_a (mantissa() As _Unsigned Long, n As Long, k As Long)
Dim As _Unsigned Long v1, v2, c1, c2
Dim As Long i
c1 = shift_constants(k - 1)
c2 = shift_constants(8 - k)
For i = 0 To n
v1 = mantissa(i) Mod c2
v2 = mantissa(i + 1) \ c2
mantissa(i) = v1 * c1 + v2
mantissa(i + 1) = mantissa(i + 1) Mod c2
Next
mantissa(n) = c1 * (mantissa(n) Mod c2)
End Sub
Sub RSHIFT_a (mantissa() As _Unsigned Long, n As Long, k As Long)
Dim As Long i
If k = 9 Then
For i = n To 1 Step -1
mantissa(i) = mantissa(i - 1)
Next
mantissa(0) = 0
Exit Sub
Else
Dim As _Unsigned Long v1, v2, c1, c2
c1 = shift_constants(k - 1)
c2 = shift_constants(8 - k)
For i = n To 1 Step -1
v1 = mantissa(i) \ c1
v2 = mantissa(i - 1) Mod c1
v2 = v2 * c2 + v1
mantissa(i) = v2
Next
mantissa(0) = mantissa(0) \ c1
End If
End Sub
Sub LSHIFT_da (mantissa() As Double, n As Long, k As Long)
Dim As _Unsigned Long v1, v2, c1, c2
Dim As Long i
c1 = shift_constants(k - 1)
c2 = shift_constants(3 - k)
For i = 2 To n - 1
v1 = mantissa(i) Mod c2
v2 = mantissa(i + 1) \ c2
mantissa(i) = v1 * c1 + v2
mantissa(i + 1) = Fix(mantissa(i + 1)) Mod c2
Next
mantissa(n) = c1 * (mantissa(n) Mod c2)
End Sub
Function fpcmp& (x As decfloat, y As decfloat)
Dim As Long c
If x.sign < y.sign Then
fpcmp& = -1: Exit Function
End If
If x.sign > y.sign Then
fpcmp& = 1: Exit Function
End If
If x.exponent < y.exponent Then
If x.sign = 0 Then
fpcmp& = -1: Exit Function
Else
fpcmp& = 1: Exit Function
End If
End If
If x.exponent > y.exponent Then
If x.sign = 0 Then
fpcmp& = 1: Exit Function
Else
fpcmp& = -1: Exit Function
End If
End If
c = x.M0 - y.M0
If c <> 0 Then GoTo fpcmpcompare
c = x.M1 - y.M1
If c <> 0 Then GoTo fpcmpcompare
c = x.M2 - y.M2
If c <> 0 Then GoTo fpcmpcompare
c = x.M3 - y.M3
If c <> 0 Then GoTo fpcmpcompare
c = x.M4 - y.M4
If c <> 0 Then GoTo fpcmpcompare
c = x.M5 - y.M5
If c <> 0 Then GoTo fpcmpcompare
c = x.M6 - y.M6
fpcmpcompare:
If c = 0 Then
fpcmp& = 0: Exit Function
End If
If c < 0 Then
If x.sign = 0 Then
fpcmp& = -1: Exit Function
Else
fpcmp& = 1: Exit Function
End If
End If
If c > 0 Then
If x.sign = 0 Then
fpcmp& = 1: Exit Function
Else
fpcmp& = -1: Exit Function
End If
End If
End Function
Sub NORM_FAC1 (fac1 As decfloat)
Dim As Long er, f
' 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.
er = 0: f = 0
If fac1.M0 > 0 Then f = 1
If fac1.M1 > 0 Then f = 1
If fac1.M2 > 0 Then f = 1
If fac1.M3 > 0 Then f = 1
If fac1.M4 > 0 Then f = 1
If fac1.M5 > 0 Then f = 1
If fac1.M6 > 0 Then f = 1
If f = 0 Then
fac1.exponent = 0
Exit Sub
'if the highmost Digit in fac1_man is nonzero,
'shift the mantissa right 1 Digit and
'increment the exponent
ElseIf fac1.M0 > 999999999 Then
RSHIFT_n fac1, 1
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 fac1.M0 = 0
LSHIFT_n fac1, 9
fac1.exponent = fac1.exponent - 9
If fac1.exponent <= 0 Then
Print "NORM_FAC1=EXPU_ERR"
Exit Sub
End If
Wend
If fac1.M0 < 10 Then
LSHIFT_n fac1, 8
fac1.exponent = fac1.exponent - 8
ElseIf fac1.M0 < 100 Then
LSHIFT_n fac1, 7
fac1.exponent = fac1.exponent - 7
ElseIf fac1.M0 < 1000 Then
LSHIFT_n fac1, 6
fac1.exponent = fac1.exponent - 6
ElseIf fac1.M0 < 10000 Then
LSHIFT_n fac1, 5
fac1.exponent = fac1.exponent - 5
ElseIf fac1.M0 < 100000 Then
LSHIFT_n fac1, 4
fac1.exponent = fac1.exponent - 4
ElseIf fac1.M0 < 1000000 Then
LSHIFT_n fac1, 3
fac1.exponent = fac1.exponent - 3
ElseIf fac1.M0 < 10000000 Then
LSHIFT_n fac1, 2
fac1.exponent = fac1.exponent - 2
ElseIf fac1.M0 < 100000000 Then
LSHIFT_n fac1, 1
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)
Dim As Long v, c
c = 0
v = fac2.M6 + fac1.M6 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M6 = v
v = fac2.M5 + fac1.M5 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M5 = v
v = fac2.M4 + fac1.M4 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M4 = v
v = fac2.M3 + fac1.M3 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M3 = v
v = fac2.M2 + fac1.M2 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M2 = v
v = fac2.M1 + fac1.M1 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M1 = v
v = fac1.M0 + fac2.M0 + c
fac1.M0 = v
NORM_FAC1 fac1
End Sub
Sub fpsub_aux (fac1 As decfloat, fac2 As decfloat)
Dim As Long v, c
c = 0
v = fac1.M6 - fac2.M6 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M6 = v
v = fac1.M5 - fac2.M5 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M5 = v
v = fac1.M4 - fac2.M4 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M4 = v
v = fac1.M3 - fac2.M3 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M3 = v
v = fac1.M2 - fac2.M2 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M2 = v
v = fac1.M1 - fac2.M1 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M1 = v
v = fac1.M0 - fac2.M0 - c
fac1.M0 = v
NORM_FAC1 fac1
End Sub
Sub fpadd (result As decfloat, x As decfloat, y As decfloat)
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)
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 &H7FFF) - BIAS - 1) - ((fac2.exponent And &H7FFF) - BIAS - 1)
If t < 63 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 < 63 Then 'shift
i = t \ 9
While i > 0
RSHIFT_n fac2, 9
t = t - 9
i = i - 1
Wend
If t = 8 Then
RSHIFT_n fac2, 8
ElseIf t = 7 Then
RSHIFT_n fac2, 7
ElseIf t = 6 Then
RSHIFT_n fac2, 6
ElseIf t = 5 Then
RSHIFT_n fac2, 5
ElseIf t = 4 Then
RSHIFT_n fac2, 4
ElseIf t = 3 Then
RSHIFT_n fac2, 3
ElseIf t = 2 Then
RSHIFT_n fac2, 2
ElseIf t = 1 Then
RSHIFT_n fac2, 1
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
fpadd_aux fac1, fac2
Else
fpsub_aux fac1, fac2
End If
End If
result = fac1
'now check exponent of result, watch for overflow.
If result.exponent > 32768 Then
'er = EXPO_ERR
Print "add: Exponent over/under flow"
End If
End Sub
Sub fpsub (result As decfloat, x As decfloat, y As decfloat)
Dim As decfloat fac1, fac2
fac1 = x
fac2 = y
fac2.sign = fac2.sign Xor &H8000
fpadd fac1, fac1, fac2
result = fac1
'now check exponent of result, watch for overflow.
If result.exponent > 32768 Then
'er = EXPO_ERR
Print "sub: Exponent over/under flow"
End If
End Sub
Sub fpmul (result As decfloat, x As decfloat, y As decfloat)
'Dim As decfloat fac1,fac2
Dim As Integer i, j, ex, den, num
Dim As _Integer64 digit, carry, prod
Dim As _Unsigned Long fac3(0 To DWORDS * 2 + 1)
Dim As Long fac1(0 To DWORDS - 1), fac2(0 To DWORDS - 1)
Dim As _Unsigned Integer fac1exponent, fac2exponent
Dim As Integer fac1sign, fac2sign
' fac1=x
' fac2=y
'check exponents. if either is zero,
'the result is zero
If x.exponent = 0 Or y.exponent = 0 Then 'result is zero...clear fac1.
result.exponent = 0
result.sign = 0
result.M0 = 0
result.M1 = 0
result.M2 = 0
result.M3 = 0
result.M4 = 0
result.M5 = 0
result.M6 = 0
Exit Sub
Else
fac1(0) = x.M0
fac1(1) = x.M1
fac1(2) = x.M2
fac1(3) = x.M3
fac1(4) = x.M4
fac1(5) = x.M5
fac1(6) = x.M6
fac1exponent = x.exponent
fac1sign = x.sign
fac2(0) = y.M0
fac2(1) = y.M1
fac2(2) = y.M2
fac2(3) = y.M3
fac2(4) = y.M4
fac2(5) = y.M5
fac2(6) = y.M6
fac2exponent = y.exponent
fac2sign = y.sign
'clear fac3 mantissa
For i = 0 To DWORDS * 2 + 1
fac3(i) = 0
Next
den = DWORDS - 1
While fac2(den) = 0
den = den - 1
Wend
num = DWORDS - 1
While fac1(num) = 0
num = num - 1
Wend
If num < den Then
'Swap fac1, fac2
fac2(0) = x.M0
fac2(1) = x.M1
fac2(2) = x.M2
fac2(3) = x.M3
fac2(4) = x.M4
fac2(5) = x.M5
fac2(6) = x.M6
fac2exponent = x.exponent
fac2sign = x.sign
fac1(0) = y.M0
fac1(1) = y.M1
fac1(2) = y.M2
fac1(3) = y.M3
fac1(4) = y.M4
fac1(5) = y.M5
fac1(6) = y.M6
fac1exponent = y.exponent
fac1sign = y.sign
Swap den, num
End If
For j = den To 0 Step -1
carry = 0
digit = fac2(j)
For i = num To 0 Step -1
prod = fac3(i + j + 1) + digit * fac1(i) + carry
carry = prod \ 1000000000
fac3(i + j + 1) = (prod Mod 1000000000)
Next
fac3(j) = carry
Next
j = 0
If fac3(0) < 10~& Then
j = 8
ElseIf fac3(0) < 100~& Then
j = 7
ElseIf fac3(0) < 1000~& Then
j = 6
ElseIf fac3(0) < 10000~& Then
j = 5
ElseIf fac3(0) < 100000~& Then
j = 4
ElseIf fac3(0) < 1000000~& Then
j = 3
ElseIf fac3(0) < 10000000~& Then
j = 2
ElseIf fac3(0) < 100000000~& Then
j = 1
End If
If j > 0 Then
LSHIFT_a fac3(), 7, j
End If
result.M0 = fac3(0)
result.M1 = fac3(1)
result.M2 = fac3(2)
result.M3 = fac3(3)
result.M4 = fac3(4)
result.M5 = fac3(5)
result.M6 = fac3(6)
End If
'now determine exponent of result.
'as you do...watch for overflow.
ex = x.exponent - BIAS + y.exponent
result.exponent = ex - j
If fac3(DWORDS) >= 500000000~& Then
Dim As decfloat fac
fac.exponent = ex - j
fac.sign = 0
fac.M0 = 0
fac.M1 = 0
fac.M2 = 0
fac.M3 = 0
fac.M4 = 0
fac.M5 = 0
fac.M6 = 1
fpadd_aux result, fac
End If
'determine the sign of the product
result.sign = x.sign Xor y.sign
NORM_FAC1 result
'now check exponent of result, watch for overflow.
If result.exponent > 32768 Then
'er = EXPO_ERR
Print "mul: Exponent over/under flow"
End If
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_out As decfloat, x As decfloat, y As decfloat)
Dim As Long fac1(6), fac2(6), r
Dim As Long i, er, is_power_of_ten, rndup
Dim As _Unsigned Integer fac1exponent, fac2exponent
Dim As Integer fac1sign, fac2sign
fac1(0) = x.M0
fac1(1) = x.M1
fac1(2) = x.M2
fac1(3) = x.M3
fac1(4) = x.M4
fac1(5) = x.M5
fac1(6) = x.M6
fac1exponent = x.exponent
fac1sign = x.sign
fac2(0) = y.M0
fac2(1) = y.M1
fac2(2) = y.M2
fac2(3) = y.M3
fac2(4) = y.M4
fac2(5) = y.M5
fac2(6) = y.M6
fac2exponent = y.exponent
fac2sign = y.sign
If fac2exponent = 0 Then ' if fac2 = 0, return
' a divide-by-zero error and
' bail out.
result_out.M0 = 999999999
result_out.M1 = 999999999
result_out.M2 = 999999999
result_out.M3 = 999999999
result_out.M4 = 999999999
result_out.M5 = 999999999
result_out.M6 = 999999999
result_out.exponent = 9999 + BIAS + 1
'er = DIVZ_ERR
Print "division by zero"
Exit Sub
ElseIf fac1exponent = 0 Then 'fact1=0, just return
er = 0
result_out.M0 = 0
result_out.M1 = 0
result_out.M2 = 0
result_out.M3 = 0
result_out.M4 = 0
result_out.M5 = 0
result_out.M6 = 0
result_out.exponent = 0
result_out.sign = 0
Exit Sub
Else
'check to see if fac2 is a power of ten
is_power_of_ten = 0
If fac2(0) = 100000000 Then
is_power_of_ten = 1
For i = 1 To 6
If fac2(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
result_out.sign = fac1sign Xor fac2sign
result_out.exponent = fac1exponent - fac2exponent + BIAS + 1
result_out.M0 = fac1(0)
result_out.M1 = fac1(1)
result_out.M2 = fac1(2)
result_out.M3 = fac1(3)
result_out.M4 = fac1(4)
result_out.M5 = fac1(5)
result_out.M6 = fac1(6)
Exit Sub
End If
Const dm = 18
Dim As Double result(1 To dm), n(1 To dm), d(1 To dm)
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 6
' n(2 * j + 2) = fac1(j) \ 10000
' n(2 * j + 3) = fac1(j) Mod 10000
' d(2 * j + 2) = fac2(j) \ 10000
' d(2 * j + 3) = fac2(j) Mod 10000
' Next
n(2) = x.M0 \ 100000: r = x.M0 Mod 100000
n(3) = r \ 10: r = r Mod 10
n(4) = r * 1000 + x.M1 \ 1000000: r = x.M1 Mod 1000000
n(5) = r \ 100: r = r Mod 100
n(6) = r * 100 + x.M2 \ 10000000: r = x.M2 Mod 10000000
n(7) = r \ 1000: r = r Mod 1000
n(8) = r * 10 + x.M3 \ 100000000: r = x.M3 Mod 100000000
n(9) = r \ 10000: r = r Mod 10000
n(10) = r
n(11) = x.M4 \ 100000: r = x.M4 Mod 100000
n(12) = r \ 10: r = r Mod 10
n(13) = r * 1000 + x.M5 \ 1000000: r = x.M5 Mod 1000000
n(14) = r \ 100: r = r Mod 100
n(15) = r * 100 + x.M6 \ 10000000: r = x.M6 Mod 10000000
n(16) = r \ 1000: r = r Mod 1000
n(17) = r * 10
d(2) = y.M0 \ 100000: r = y.M0 Mod 100000
d(3) = r \ 10: r = r Mod 10
d(4) = r * 1000 + y.M1 \ 1000000: r = y.M1 Mod 1000000
d(5) = r \ 100: r = r Mod 100
d(6) = r * 100 + y.M2 \ 10000000: r = y.M2 Mod 10000000
d(7) = r \ 1000: r = r Mod 1000
d(8) = r * 10 + y.M3 \ 100000000: r = y.M3 Mod 100000000
d(9) = r \ 10000: r = r Mod 10000
d(10) = r
d(11) = y.M4 \ 100000: r = y.M4 Mod 100000
d(12) = r \ 10: r = r Mod 10
d(13) = r * 1000 + y.M5 \ 1000000: r = y.M5 Mod 1000000
d(14) = r \ 100: r = r Mod 100
d(15) = r * 100 + y.M6 \ 10000000: r = y.M6 Mod 10000000
d(16) = r \ 1000: r = r Mod 1000
d(17) = r * 10
n(1) = (fac1exponent And &H7FFF) - BIAS - 1
d(1) = (fac2exponent And &H7FFF) - 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))
subtract w(), q, d(), (stp + 2), last
normalize w(), (stp + 2), q
Next
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)
j = 0
If result(2) < 10~& Then
j = 3
ElseIf result(2) < 100~& Then
j = 2
ElseIf result(2) < 1000~& Then
j = 1
End If
If j > 0 Then
LSHIFT_da result(), dm, j
End If
j = 2
result_out.M0 = (result(j) * 10000 + result(j + 1)) * 10 + result(j + 2) \ 1000: r = result(j + 2) Mod 1000: j = 5
result_out.M1 = (r * 10000 + result(j)) * 100 + result(j + 1) \ 100: r = result(j + 1) Mod 100: j = 7
result_out.M2 = (r * 10000 + result(j)) * 1000 + result(j + 1) \ 10: r = result(j + 1) Mod 10: j = 9
result_out.M3 = (r * 10000 + result(j)) * 10000 + result(j + 1): j = 11
result_out.M4 = (result(j) * 10000 + result(j + 1)) * 10 + result(j + 2) \ 1000: r = result(j + 2) Mod 1000: j = 14
result_out.M5 = (r * 10000 + result(j)) * 100 + result(j + 1) \ 100: r = result(j + 1) Mod 100: j = 16
result_out.M6 = (r * 10000 + result(j)) * 1000 + result(j + 1) \ 10
rndup = 0
If (result(j + 1) Mod 10) > 4 Then
Dim As decfloat ru
ru.M0 = 100000000
ru.exponent = 1 + BIAS
result_out.exponent = 63 + BIAS
fpadd result_out, result_out, ru
If result_out.exponent > (63 + BIAS) Then
rndup = 1
End If
End If
NORM_FAC1 result_out
result_out.exponent = (result(1) + BIAS) + rndup
End If
result_out.sign = fac1sign Xor fac2sign
End Sub
Sub fpdiv_si (result As decfloat, num As decfloat, den As Long)
Dim As Long fac1(6)
Dim As _Unsigned Integer fac1exponent
Dim As Integer fac1sign
Dim As _Unsigned _Integer64 carry, remder
Dim As _Integer64 i, divisor
Dim As _Integer64 quotient
remder = 0
divisor = Abs(den)
fac1(0) = num.M0
fac1(1) = num.M1
fac1(2) = num.M2
fac1(3) = num.M3
fac1(4) = num.M4
fac1(5) = num.M5
fac1(6) = num.M6
fac1exponent = num.exponent
fac1sign = num.sign
result.M0 = num.M0
result.M1 = num.M1
result.M2 = num.M2
result.M3 = num.M3
result.M4 = num.M4
result.M5 = num.M5
result.M6 = num.M6
result.exponent = num.exponent
result.sign = num.sign
If divisor = 0 Then
Print "error: divisor = 0"
Exit Sub
End If
If divisor > 2147483647 Then
Print "error: divisor too large"
Exit Sub
End If
For i = 0 To 6
quotient = fac1(i) + remder * 1000000000
remder = quotient Mod divisor
fac1(i) = quotient \ divisor
Next
quotient = remder * 1000000000
quotient = quotient \ divisor
carry = fac1(0)
result.M0 = fac1(0)
result.M1 = fac1(1)
result.M2 = fac1(2)
result.M3 = fac1(3)
result.M4 = fac1(4)
result.M5 = fac1(5)
result.M6 = fac1(6)
result.exponent = fac1exponent
result.sign = fac1sign
If carry = 0 Then
LSHIFT_n result, 9
result.exponent = result.exponent - 9
result.M6 = result.M6 + quotient
ElseIf carry < 10 Then
LSHIFT_n result, 8
result.exponent = result.exponent - 8
result.M6 = result.M6 + quotient \ 10
ElseIf carry < 100 Then
LSHIFT_n result, 7
result.exponent = result.exponent - 7
result.M6 = result.M6 + quotient \ 100
ElseIf carry < 1000 Then
LSHIFT_n result, 6
result.exponent = result.exponent - 6
result.M6 = result.M6 + quotient \ 1000
ElseIf carry < 10000 Then
LSHIFT_n result, 5
result.exponent = result.exponent - 5
result.M6 = result.M6 + quotient \ 10000
ElseIf carry < 100000 Then
LSHIFT_n result, 4
result.exponent = result.exponent - 4
result.M6 = result.M6 + quotient \ 100000
ElseIf carry < 1000000 Then
LSHIFT_n result, 3
result.exponent = result.exponent - 3
result.M6 = result.M6 + quotient \ 1000000
ElseIf carry < 10000000 Then
LSHIFT_n result, 2
result.exponent = result.exponent - 2
result.M6 = result.M6 + quotient \ 10000000
ElseIf carry < 100000000 Then
LSHIFT_n result, 1
result.exponent = result.exponent - 1
result.M6 = result.M6 + quotient \ 100000000
End If
'NORM_FAC1(fac1)
result.sign = fac1sign
If den < 0 Then
result.sign = fac1sign Xor &H8000
End If
'now check exponent of result, watch for overflow.
If result.exponent > 32768 Then
'er = EXPO_ERR
Print "div_si: Exponent over/under flow"
End If
End Sub
Sub fpmul_si (result As decfloat, x As decfloat, y As _Integer64)
Dim As decfloat fac1, fac2
Dim As Long i
Dim As _Integer64 carry, digit, prod, value
Dim As _Unsigned Long fac3(7)
fac1 = x
digit = Abs(y)
If digit > 999999999 Then
Si2fp fac2, y
fpmul fac1, fac1, fac2
result = fac1: Exit Sub
End If
If fac1.exponent = 0 Or y = 0 Then 'result is zero...clear fac1.
fac1.sign = 0
fac1.exponent = 0
fac1.M0 = 0
fac1.M1 = 0
fac1.M2 = 0
fac1.M3 = 0
fac1.M4 = 0
fac1.M5 = 0
fac1.M6 = 0
NORM_FAC1 fac1
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
fac3(0) = fac1.M0
fac3(1) = fac1.M1
fac3(2) = fac1.M2
fac3(3) = fac1.M3
fac3(4) = fac1.M4
fac3(5) = fac1.M5
fac3(6) = fac1.M6
fac3(7) = 0
carry = 0
For i = 6 To 0 Step -1
prod = digit * fac3(i) + carry
value = prod Mod 1000000000
fac3(i) = value
carry = prod \ 1000000000
Next
If carry < 10 Then
RSHIFT_a fac3(), 7, 1
fac1.exponent = fac1.exponent + 1
fac3(0) = fac3(0) + carry * 100000000
ElseIf carry < 100 Then
RSHIFT_a fac3(), 7, 2
fac1.exponent = fac1.exponent + 2
fac3(0) = fac3(0) + carry * 10000000
ElseIf carry < 1000 Then
RSHIFT_a fac3(), 7, 3
fac1.exponent = fac1.exponent + 3
fac3(0) = fac3(0) + carry * 1000000
ElseIf carry < 10000 Then
RSHIFT_a fac3(), 7, 4
fac1.exponent = fac1.exponent + 4
fac3(0) = fac3(0) + carry * 100000
ElseIf carry < 100000 Then
RSHIFT_a fac3(), 7, 5
fac1.exponent = fac1.exponent + 5
fac3(0) = fac3(0) + carry * 10000
ElseIf carry < 1000000 Then
RSHIFT_a fac3(), 7, 6
fac1.exponent = fac1.exponent + 6
fac3(0) = fac3(0) + carry * 1000
ElseIf carry < 10000000 Then
RSHIFT_a fac3(), 7, 7
fac1.exponent = fac1.exponent + 7
fac3(0) = fac3(0) + carry * 100
ElseIf carry < 100000000 Then
RSHIFT_a fac3(), 7, 8
fac1.exponent = fac1.exponent + 8
fac3(0) = fac3(0) + carry * 10
ElseIf carry < 1000000000 Then
RSHIFT_a fac3(), 7, 9
fac1.exponent = fac1.exponent + 9
fac3(0) = fac3(0) + carry
End If
End If
fac1.M0 = fac3(0)
fac1.M1 = fac3(1)
fac1.M2 = fac3(2)
fac1.M3 = fac3(3)
fac1.M4 = fac3(4)
fac1.M5 = fac3(5)
fac1.M6 = fac3(6)
If fac3(7) >= 500000000~& Then
Dim As decfloat fac
fac.exponent = fac1.exponent
fac.sign = 0
fac.M0 = 0
fac.M1 = 0
fac.M2 = 0
fac.M3 = 0
fac.M4 = 0
fac.M5 = 0
fac.M6 = 1
fpadd_aux fac1, fac
End If
NORM_FAC1 fac1
If y < 0 Then
fac1.sign = fac1.sign Xor &H8000
End If
result = fac1
'now check exponent of result, watch for overflow.
If result.exponent > 32768 Then
'er = EXPO_ERR
Print "mul_si: Exponent over/under flow"
End If
End Sub
'integer part of num
Sub fpfix (result As decfloat, num As decfloat)
Dim As decfloat ips
Dim As Long ex, ex2, j, k
ex = (num.exponent And &H7FFF) - BIAS
If ex < 1 Then
result = ips
Exit Sub
End If
If ex >= 63 Then
result = num
Exit Sub
End If
Dim As _Unsigned Long ip(6), numa(6)
numa(0) = num.M0
numa(1) = num.M1
numa(2) = num.M2
numa(3) = num.M3
numa(4) = num.M4
numa(5) = num.M5
numa(6) = num.M6
ex2 = ex \ 9
k = ex2
j = ex Mod 9
While ex2 > 0
ex2 = ex2 - 1
ip(ex2) = numa(ex2)
Wend
If j = 1 Then
ip(k) = 100000000 * (numa(k) \ 100000000)
ElseIf j = 2 Then
ip(k) = 10000000 * (numa(k) \ 10000000)
ElseIf j = 3 Then
ip(k) = 1000000 * (numa(k) \ 1000000)
ElseIf j = 4 Then
ip(k) = 100000 * (numa(k) \ 100000)
ElseIf j = 5 Then
ip(k) = 10000 * (numa(k) \ 10000)
ElseIf j = 6 Then
ip(k) = 1000 * (numa(k) \ 1000)
ElseIf j = 7 Then
ip(k) = 100 * (numa(k) \ 100)
ElseIf j = 8 Then
ip(k) = 10 * (numa(k) \ 10)
ElseIf j = 9 Then
ip(k) = numa(k)
End If
result.exponent = ex + BIAS
result.sign = num.sign
result.M0 = ip(0)
result.M1 = ip(1)
result.M2 = ip(2)
result.M3 = ip(3)
result.M4 = ip(4)
result.M5 = ip(5)
result.M6 = ip(6)
End Sub
Sub fpfrac (result As decfloat, num As decfloat)
Dim As decfloat ip, fp
fpfix ip, num
fpsub fp, num, ip
result = fp
End Sub
Function fpfix_is_odd& (num As decfloat)
Dim As Long ex, ex2, j, k
ex = (num.exponent And &H7FFF) - BIAS
If ex < 0 Then
fpfix_is_odd = 0
Exit Function
End If
If ex >= 63 Then
fpfix_is_odd = 0
Exit Function
End If
Dim As _Unsigned Long numa(6)
numa(0) = num.M0
numa(1) = num.M1
numa(2) = num.M2
numa(3) = num.M3
numa(4) = num.M4
numa(5) = num.M5
numa(6) = num.M6
ex = ex - 1
ex2 = ex \ 9
k = ex2
j = ex Mod 9
If j = 0 Then
fpfix_is_odd = (numa(k) \ 100000000) And 1: Exit Function
ElseIf j = 1 Then
fpfix_is_odd = (numa(k) \ 10000000) And 1: Exit Function
ElseIf j = 2 Then
fpfix_is_odd = (numa(k) \ 1000000) And 1: Exit Function
ElseIf j = 3 Then
fpfix_is_odd = (numa(k) \ 100000) And 1: Exit Function
ElseIf j = 4 Then
fpfix_is_odd = (numa(k) \ 10000) And 1: Exit Function
ElseIf j = 5 Then
fpfix_is_odd = (numa(k) \ 1000) And 1: Exit Function
ElseIf j = 6 Then
fpfix_is_odd = (numa(k) \ 100) And 1: Exit Function
ElseIf j = 7 Then
fpfix_is_odd = (numa(k) \ 10) And 1: Exit Function
ElseIf j = 8 Then
fpfix_is_odd = numa(k) And 1: Exit Function
End If
fpfix_is_odd = 0
End Function
sample run
Code: (Select All) 1 1.000000000000000000000000000000000000000000000000000000000 rel. error 0
2 1.414213562373095048801688724209698078569671875376948073177 rel. error 5.000000000e-63
3 1.732050807568877293527446341505872366942805253810380628056 rel. error 3.333333333e-63
4 2.000000000000000000000000000000000000000000000000000000000 rel. error 0
5 2.236067977499789696409173668731276235440618359611525724271 rel. error 6.000000000e-63
6 2.449489742783178098197284074705891391965947480656670128433 rel. error 1.666666667e-63
7 2.645751311064590590501615753639260425710259183082450180368 rel. error 0
8 2.828427124746190097603377448419396157139343750753896146353 rel. error 3.750000000e-63
9 3.000000000000000000000000000000000000000000000000000000000 rel. error 0
10 3.162277660168379331998893544432718533719555139325216826858 rel. error 0
11 3.316624790355399849114932736670686683927088545589353597059 rel. error 0
12 3.464101615137754587054892683011744733885610507620761256112 rel. error 0
13 3.605551275463989293119221267470495946251296573845246212710 rel. error 0
14 3.741657386773941385583748732316549301756019807778726946304 rel. error 0
15 3.872983346207416885179265399782399610832921705291590826588 rel. error 0
16 4.000000000000000000000000000000000000000000000000000000000 rel. error 0
17 4.123105625617660549821409855974077025147199225373620434399 rel. error 5.882352941e-63
18 4.242640687119285146405066172629094235709015626130844219530 rel. error 5.555555556e-63
19 4.358898943540673552236981983859615659137003925232444936890 rel. error 0
20 4.472135954999579392818347337462552470881236719223051448542 rel. error 0
Posts: 2,186
Threads: 222
Joined: Apr 2022
Reputation:
104
Hi Jack,
Bit shifting is something I haven't made familiar to me. What is familiar to me is optimization large routines. Since you are getting into lots of lines of code, would this modification be of any benefit to you?
Swap out IF/THEN routines like this...
Code: (Select All) If carry = 0 Then
LSHIFT_n result, 9
result.exponent = result.exponent - 9
result.M6 = result.M6 + quotient
ElseIf carry < 10 Then
LSHIFT_n result, 8
result.exponent = result.exponent - 8
result.M6 = result.M6 + quotient \ 10
ElseIf carry < 100 Then
LSHIFT_n result, 7
result.exponent = result.exponent - 7
result.M6 = result.M6 + quotient \ 100
ElseIf carry < 1000 Then
LSHIFT_n result, 6
result.exponent = result.exponent - 6
result.M6 = result.M6 + quotient \ 1000
ElseIf carry < 10000 Then
LSHIFT_n result, 5
result.exponent = result.exponent - 5
result.M6 = result.M6 + quotient \ 10000
ElseIf carry < 100000 Then
LSHIFT_n result, 4
result.exponent = result.exponent - 4
result.M6 = result.M6 + quotient \ 100000
ElseIf carry < 1000000 Then
LSHIFT_n result, 3
result.exponent = result.exponent - 3
result.M6 = result.M6 + quotient \ 1000000
ElseIf carry < 10000000 Then
LSHIFT_n result, 2
result.exponent = result.exponent - 2
result.M6 = result.M6 + quotient \ 10000000
ElseIf carry < 100000000 Then
LSHIFT_n result, 1
result.exponent = result.exponent - 1
result.M6 = result.M6 + quotient \ 100000000
End If
With an algorithm like this...
Code: (Select All) j& = LEN(LTRIM$(STR$(carry)))
LSHIFT_n result, 10 - j&
result.exponent = result.exponent - 10 - j&
result.M6 = result.M6 + quotient \ carry
Example of number output...
Code: (Select All) DIM carry AS DOUBLE
FOR i = 0 TO 8
carry = 10 ^ i
j = LEN(LTRIM$(STR$(carry)))
'LSHIFT_n result, 10 - j
PRINT "carry ="; carry;: LOCATE , 27: PRINT "j ="; j, "10 - j ="; 10 - j
result.exponent = result.exponent - 10 - j
result.M6 = result.M6 + quotient \ LEN(LTRIM$(STR$(carry)))
NEXT
The flip side is the longer method is easier to read through, but for speed sake, the algorithm would be faster than running through all the those conditional statements.
Sorry I couldn't test my code out in the routine; what you are doing is a bit flipping beyond me. I think I code it correctly, but if you want to use the idea, please test it yourself. I wouldn't want to set back your progress.
Thanks for all your help with my iteration routines. You really have some mad math skills!
Pete
Posts: 422
Threads: 27
Joined: Apr 2022
Reputation:
26
10-09-2022, 07:21 PM
(This post was last modified: 10-09-2022, 07:26 PM by Jack.)
thanks Pete for your suggestions 
right now I am busy with other things, will look into it another day
P.S.
using select case might be better
Posts: 422
Threads: 27
Joined: Apr 2022
Reputation:
26
10-12-2022, 01:00 AM
(This post was last modified: 10-12-2022, 05:35 AM by Jack.)
I am trying to optimize the routines to get the maximum accuracy but I am not satisfied with the result so far
if I roundup the multiply sub then the result tends to be a bit larger than it should be, and if I don't roundup then it's less than it should be
please give this a test run and see what you think
Code: (Select All) _Title "defloat-32 v2 by Jack"
$NoPrefix
$Console:Only
Dest Console
Dim Shared As Long shift_constants(7)
shift_constants(0) = 10
shift_constants(1) = 100
shift_constants(2) = 1000
shift_constants(3) = 10000
shift_constants(4) = 100000
shift_constants(5) = 1000000
shift_constants(6) = 10000000
shift_constants(7) = 100000000
Const BIAS = 16384 '2 ^ 14
Const DWORDS = 7
Type decfloat
As _Unsigned Integer exponent
As Integer sign
As _Unsigned Long M0
As _Unsigned Long M1
As _Unsigned Long M2
As _Unsigned Long M3
As _Unsigned Long M4
As _Unsigned Long M5
As _Unsigned Long M6
End Type
Dim Shared As decfloat pi_dec, pi2_dec, pi_dec_half, xpn, xpd
Dim Shared As Long multrndup
str2fp pi_dec, "3.14159265358979323846264338327950288419716939937510582097494459"
str2fp pi2_dec, "6.28318530717958647692528676655900576839433879875021164194988918"
str2fp pi_dec_half, "1.57079632679489661923132169163975144209858469968755291048747230"
' Error definitions
Const DIVZ_ERR = 1 'Divide by zero
Const EXPO_ERR = 2 'Exponent overflow error
Const EXPU_ERR = 3 'Exponent underflow error
Dim As decfloat x, y, z, p
Dim As Long j
Dim As Double t
Print "without rounding up multiply"
multrndup = 0
str2fp x, "2"
Print "x = 2"
fplog y, x
t = Timer(.0001)
For j = 1 To 10000
fpexp z, y
Next
t = Timer(.0001) - t
Print "log(x) = "; fp2str_exp(y, 63)
Print "exp(log(x)) = "; fp2str_exp(z, 63)
fpsub p, z, x
fpdiv p, p, x
Print "rel. error of exp(log(x)) = "; fp2str(p, 19)
Print "elapsed time for 10000 exp calls = "; t
Print
str2fp x, "999999999123456789123456789123456789123456789123456789123456789"
Print "x = 999999999123456789123456789123456789123456789123456789123456789"
fplog y, x
t = Timer(.0001)
For j = 1 To 10000
fpexp z, y
Next
t = Timer(.0001) - t
Print "log(x) = "; fp2str_exp(y, 63)
Print "exp(log(x)) = "; fp2str_exp(z, 63)
fpsub p, z, x
fpdiv p, p, x
Print "rel. error of exp(log(x)) = "; fp2str(p, 19)
Print "elapsed time for 10000 exp calls = "; t
Print
Print "with rounding up multiply"
multrndup = 1
str2fp x, "2"
Print "x = 2"
fplog y, x
t = Timer(.0001)
For j = 1 To 10000
fpexp z, y
Next
t = Timer(.0001) - t
Print "log(x) = "; fp2str_exp(y, 63)
Print "exp(log(x)) = "; fp2str_exp(z, 63)
fpsub p, z, x
fpdiv p, p, x
Print "rel. error of exp(log(x)) = "; fp2str(p, 19)
Print "elapsed time for 10000 exp calls = "; t
Print
str2fp x, "999999999123456789123456789123456789123456789123456789123456789"
Print "x = 999999999123456789123456789123456789123456789123456789123456789"
fplog y, x
t = Timer(.0001)
For j = 1 To 10000
fpexp z, y
Next
t = Timer(.0001) - t
Print "log(x) = "; fp2str_exp(y, 63)
Print "exp(log(x)) = "; fp2str_exp(z, 63)
fpsub p, z, x
fpdiv p, p, x
Print "rel. error of exp(log(x)) = "; fp2str(p, 19)
Print "elapsed time for 10000 exp calls = "; t
Sub fpexp_cf (result As decfloat, z As decfloat)
Dim As decfloat x, s, x1, x2, tmp, tmp2
Dim As Long ex, i, ip, n
Static As decfloat one, euler
n = 10
If euler.exponent = 0 Then
str2fp euler, "2.71828182845904523536028747135266249775724709369995957496696763"
one.sign = 0
one.exponent = (BIAS + 1)
one.M0 = 100000000
one.M1 = 0
one.M2 = 0
one.M3 = 0
one.M4 = 0
one.M5 = 0
one.M6 = 0
End If
' 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)
If z.exponent = 0 Then
result = one
Exit Sub
Else
x = z
fpfix tmp, x
ex = (tmp.exponent And &H7FFF) - BIAS - 1
If ex = 0 Then
ip = tmp.M0 \ 100000000
ElseIf ex = 1 Then
ip = tmp.M0 \ 10000000
ElseIf ex = 2 Then
ip = tmp.M0 \ 1000000
ElseIf ex = 3 Then
ip = tmp.M0 \ 100000
ElseIf ex = 4 Then
ip = tmp.M0 \ 10000
ElseIf ex = 5 Then
ip = tmp.M0 \ 1000
ElseIf ex = 6 Then
ip = tmp.M0 \ 100
ElseIf ex = 7 Then
ip = tmp.M0 \ 10
ElseIf ex = 8 Then
ip = tmp.M0
End If
fpfrac x, x
fpdiv_si x1, x, 16 '32768 ' 2^13
fpmul x2, x1, x1
Si2fp s, 4 * n + 6
For i = n To 0 Step -1
Si2fp tmp, 4 * i + 2
fpdiv s, x2, s
fpadd s, s, tmp
Next
fpdiv s, x1, s
Si2fp tmp, 1
tmp2 = tmp
fpadd tmp, tmp, s
fpsub tmp2, tmp2, s
fpdiv s, tmp, tmp2
For i = 1 To 4 '15 ' s^13
fpmul s, s, s
Next
fpipow tmp, euler, ip
fpmul result, s, tmp
End If
End Sub
Sub fpexp (result As decfloat, z As decfloat)
Dim As decfloat x, xpd, xpn
Dim As Long ex, i, ip
Static As decfloat one, euler, exponentialn(17), exponentiald(17)
If exponentialn(0).exponent = 0 Then
Si2fp exponentialn(0), 1
str2fp exponentialn(1), ".503787763085141830426316836835640413180555373602013787765506507"
str2fp exponentialn(2), ".123113754811115580240199391979941449439554252687431414889261917"
str2fp exponentialn(3), "0.194014786172100868397379733957773068817466004597117759916094387e-1"
str2fp exponentialn(4), "0.221053514053621305631679326515226926261161756634139701475334846e-2"
str2fp exponentialn(5), "0.193454865206744920585608214573628003845560417585149728091593248e-3"
str2fp exponentialn(6), "0.134817206102830919302898061303157717441762386334371916554674483e-4"
str2fp exponentialn(7), "7.65161710091203441510338717292779761121416735865848016278269163e-7"
str2fp exponentialn(8), "3.58548435070740017294732467393249181951311755864779732794287718e-8"
str2fp exponentialn(9), "1.39715881029815224422817956869312878658017689986058208796945882e-9"
str2fp exponentialn(10), "4.53464900275857744649901664225234895410940311133637948725625465e-11"
str2fp exponentialn(11), "1.22102800612904874865986466356114289547925418884289566733570583e-12"
str2fp exponentialn(12), "2.69942119086103755038141670655438202935697783465198117242288620e-14"
str2fp exponentialn(13), "4.80709158088050415577996944372246571546305491654413869758830148e-16"
str2fp exponentialn(14), "6.67446218117009199081966862645625972417483845805581243134455550e-18"
str2fp exponentialn(15), "6.82626745056408867477242403742610973442937848360726389682334564e-20"
str2fp exponentialn(16), "4.60485274414644615989466638049841829236889098591288845615263752e-22"
str2fp exponentialn(17), "1.54768708764447533766179612797166908030142040749901411508788783e-24"
Si2fp exponentiald(0), 1
str2fp exponentiald(1), "-.496212236914858169573683163164359586819444626397986212234493493"
str2fp exponentiald(2), ".119325991725973749813882555144301036258998879085417627123755410"
str2fp exponentiald(3), "-0.184850613180012448539696668330106026341966320933794116815659028e-1"
str2fp exponentiald(4), "0.206797341469361126562570972007228490388291451667636384019144944e-2"
str2fp exponentiald(5), "-0.177476640029445154798895481195212175295325133388353948891798456e-3"
str2fp exponentiald(6), "0.121119697791530382246392392593837579456791651154395855360677559e-4"
str2fp exponentiald(7), "-6.72136000621469478770704356909331255584862479068643438675858294e-7"
str2fp exponentiald(8), "3.07422101412642477250911345710710651015389275273741476097162835e-8"
str2fp exponentiald(9), "-1.16701132447665083402247823313377310891548863251014843319140333e-9"
str2fp exponentiald(10), "3.68187090359605484256890447350797602918878757385065579230141815e-11"
str2fp exponentiald(11), "-9.61338101757297931016053534655678826462590722321717850532044698e-13"
str2fp exponentiald(12), "2.05509049962341068403280550130851518395842965866907693627957626e-14"
str2fp exponentiald(13), "-3.52746564887524616342404030739724544161545203402406781919425177e-16"
str2fp exponentiald(14), "4.70347116684618761386137496578607769970491970387898680511311425e-18"
str2fp exponentiald(15), "-4.59998458732742996655211581291170107043066864993354030929388480e-20"
str2fp exponentiald(16), "2.95256671673728810659480905519132991941158386504354019400480885e-22"
str2fp exponentiald(17), "-9.38719670297727932196732935370899807766127816686343783729700392e-25"
str2fp euler, "2.71828182845904523536028747135266249775724709369995957496696763"
one.sign = 0
one.exponent = (BIAS + 1)
one.M0 = 100000000
one.M1 = 0
one.M2 = 0
one.M3 = 0
one.M4 = 0
one.M5 = 0
one.M6 = 0
End If
If z.exponent = 0 Then
result = one
Exit Sub
Else
x = z
fpfix xpn, x
ex = (xpn.exponent And &H7FFF) - BIAS - 1
If ex = 0 Then
ip = xpn.M0 \ 100000000
ElseIf ex = 1 Then
ip = xpn.M0 \ 10000000
ElseIf ex = 2 Then
ip = xpn.M0 \ 1000000
ElseIf ex = 3 Then
ip = xpn.M0 \ 100000
ElseIf ex = 4 Then
ip = xpn.M0 \ 10000
ElseIf ex = 5 Then
ip = xpn.M0 \ 1000
ElseIf ex = 6 Then
ip = xpn.M0 \ 100
ElseIf ex = 7 Then
ip = xpn.M0 \ 10
ElseIf ex = 8 Then
ip = xpn.M0
End If
fpfrac x, x
fpmul xpn, exponentialn(17), x
For i = 16 To 0 Step -1
fpadd xpn, xpn, exponentialn(i)
fpmul xpn, xpn, x
Next
fpmul xpd, exponentiald(17), x
For i = 16 To 0 Step -1
fpadd xpd, xpd, exponentiald(i)
fpmul xpd, xpd, x
Next
fpdiv xpn, xpn, xpd
fpipow xpd, euler, ip
fpmul result, xpn, xpd
End If
End Sub
Sub fplog (result As decfloat, z As decfloat)
Dim As decfloat x, xpd, xpn
Dim As Long ex, i
Static As decfloat fpln10, one, logarithmn(22), logarithmd(24)
If logarithmn(0).exponent = 0 Then
Si2fp logarithmn(0), 1
str2fp logarithmn(1), "10.2310520037359891285036755418469258827514540594420472830660753"
str2fp logarithmn(2), "48.6627011559331606160440612789640717584111988048005196804661728"
str2fp logarithmn(3), "142.883839025492099568442001118949943253750601849543916975383500"
str2fp logarithmn(4), "290.042754824482811137859430225980858456412815030133662556507508"
str2fp logarithmn(5), "431.988540379436291654167001685688951066167323915161481314938683"
str2fp logarithmn(6), "488.988532243997467759781233533239006273973991246082219125463314"
str2fp logarithmn(7), "429.964620903161621905726323536985462923884859411935413887634588"
str2fp logarithmn(8), "297.667076072840840028860885452874923891756123601208220339342769"
str2fp logarithmn(9), "163.496246139169410681932569760493099792107706750131034678156585"
str2fp logarithmn(10), "71.4689010653151414969354268941968334288999401577488778627880921"
str2fp logarithmn(11), "24.8454129712280032841853780591175947041028424111335032188809328"
str2fp logarithmn(12), "6.83955915425638170286810308584881831872981167696540451581896838"
str2fp logarithmn(13), "1.47913861498769669443157158341616767519966697200814961672841044"
str2fp logarithmn(14), ".248298871909766565290925597836167509193223133304509724924096267"
str2fp logarithmn(15), "0.318133082830696859180749013242710465363599980360532595003499144e-1"
str2fp logarithmn(16), "0.304020922309269050902540112467577787426750224236109865380535113e-2"
str2fp logarithmn(17), "0.209978185892809569953723710672234372599610486457723316925724789e-3"
str2fp logarithmn(18), "0.100311469452988931955197299915204806738105549379720367796921701e-4"
str2fp logarithmn(19), "3.11027381043420800967758659489427925931042335875083315411300509e-7"
str2fp logarithmn(20), "5.67148960215418390372633356121445540918492594737567360786766486e-9"
str2fp logarithmn(21), "5.12133417069955237659371135908115173603073989914114172007865451e-11"
str2fp logarithmn(22), "1.55966557818637624120176613242050057417746497195348156692036175e-13"
Si2fp logarithmd(0), 1
str2fp logarithmd(1), "10.7310520037359891285036755418469258827514540594420472830660738"
str2fp logarithmd(2), "53.6948938244678218469625657165542013664535925011882099886674915"
str2fp logarithmd(3), "166.404268603147347449088725463278068642726913413657339541359027"
str2fp logarithmd(4), "357.829354185501208195542189937563557126312937751426774299756481"
str2fp logarithmd(5), "566.879307659840870904948427820813871913477928966397711410863037"
str2fp logarithmd(6), "685.659474589076294970895474467750778199315371532367251137728620"
str2fp logarithmd(7), "647.552215345138140652347527859581928707051002114597464686830520"
str2fp logarithmd(8), "484.337597665677116149955858143764741744286855179474277002221246"
str2fp logarithmd(9), "289.339023581167087024317711351001662002419858544547812026145386"
str2fp logarithmd(10), "138.626741697453253940901869606998742571923585752594114076259249"
str2fp logarithmd(11), "53.2967984633656431994464301417545208285735718517685211650739340"
str2fp logarithmd(12), "16.3977970191605190152676712715525604280407449895571363588426976"
str2fp logarithmd(13), "4.01344869231909159178580772978304524142055416630126527640586347"
str2fp logarithmd(14), ".774140683278008004255775706220863488435679762354478168050645942"
str2fp logarithmd(15), ".116116794765397457080324911252454433970856146870870944104060064"
str2fp logarithmd(16), "0.132995134740536235985194287064770583314074985150347963060323352e-1"
str2fp logarithmd(17), "0.113495871890204061922407499675091077013725252521536694012719530e-2"
str2fp logarithmd(18), "0.697980393999574131276345270647881218831498205731172278517042043e-4"
str2fp logarithmd(19), "0.295301241701723056394124042476367677312456922839751324605065650e-5"
str2fp logarithmd(20), "8.03411378004236615745112193270089533877927829488633592988420379e-8"
str2fp logarithmd(21), "1.26487301511055907218913028822339261374104341757811651124170618e-9"
str2fp logarithmd(22), "9.55430735765744135751947902497685515801773089075935394013525279e-12"
str2fp logarithmd(23), "2.22906425441965268150705233939758319540146232049504069112205822e-14"
str2fp logarithmd(24), "-3.94659794323870669763249350331131008555188524939917666031562314e-18"
str2fp fpln10, "2.30258509299404568401799145468436420760110148862877297603332790"
one.sign = 0
one.exponent = (BIAS + 1)
one.M0 = 100000000
one.M1 = 0
one.M2 = 0
one.M3 = 0
one.M4 = 0
one.M5 = 0
one.M6 = 0
End If
x = z
If fpcmp&(x, one) = 0 Then
result.sign = 0
result.exponent = 0
result.M0 = 0
result.M1 = 0
result.M2 = 0
result.M3 = 0
result.M4 = 0
result.M5 = 0
result.M6 = 0
Exit Sub
End If
If x.exponent <> 0 Then
ex = (x.exponent And &H7FFF) - BIAS - 1
x.exponent = BIAS + 1
fpsqr x, x
fpsqr x, x
fpsub x, x, one
fpmul xpn, logarithmn(22), x
For i = 21 To 0 Step -1
fpadd xpn, xpn, logarithmn(i)
fpmul xpn, xpn, x
Next
fpmul xpn, xpn, x
fpmul xpd, logarithmd(24), x
For i = 23 To 0 Step -1
fpadd xpd, xpd, logarithmd(i)
fpmul xpd, xpd, x
Next
fpdiv xpn, xpn, xpd
fpmul_si xpn, xpn, 4
fpmul_si xpd, fpln10, ex
fpadd result, xpn, xpd
End If
End Sub
Sub fpatn (result As decfloat, x As decfloat)
Dim As Long sign(3): sign(3) = 1
Dim As Long z, 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
one.sign = 0
one.exponent = (BIAS + 1)
one.M0 = 100000000
one.M1 = 0
one.M2 = 0
one.M3 = 0
one.M4 = 0
one.M5 = 0
one.M6 = 0
If fpcmp(decnum2, one) = 0 Then
fpdiv_si decnum, pi_dec, 4
decnum.sign = x.sign
result = decnum
Exit Sub
End If
decnum2.sign = x.sign
Dim As Long lm2: lm2 = 16
Si2fp factor, _ShL(2, (lm2 - 1))
For z = 1 To lm2
fpmul decnum, decnum2, decnum2
fpadd decnum, decnum, one
fpsqr decnum, decnum
fpadd decnum, decnum, one
fpdiv decnum, decnum2, decnum
decnum2 = decnum
Next z
mt = decnum
x_2 = decnum
p = decnum
fpmul XX, x_2, x_2
Do
c = c + 2
mt2 = mt
Si2fp strC, c
fpmul p, p, XX
fpdiv Term, p, strC
If sign(c And 3) Then
fpsub Accum, mt, Term
Else
fpadd Accum, mt, Term
End If
Swap mt, Accum
Loop Until fpcmp(mt, mt2) = 0
fpmul result, factor, mt
End Sub
Sub fpasin (result As decfloat, x As decfloat)
Dim As Double num
Dim As decfloat one, T, B, term1, minusone
' ARCSIN = ATN(x / SQR(-x * x + 1))
'============= ARCSIN GUARD =========
num = Val(fp2str(x, 16))
If num > 1 Then
result = one
Exit Sub
End If
If num < -1 Then
result = one
Exit Sub
End If
'========================
one.sign = 0
one.exponent = (BIAS + 1)
one.M0 = 100000000
one.M1 = 0
one.M2 = 0
one.M3 = 0
one.M4 = 0
one.M5 = 0
one.M6 = 0
minusone = one: minusone.sign = &H8000
T = x
fpmul B, x, x 'x*x
'for 1 and -1
If fpcmp(B, one) = 0 Then
Dim As decfloat two, atn1
two.sign = 0
two.exponent = (BIAS + 1)
two.M0 = 200000000
two.M0 = 0
two.M0 = 0
two.M0 = 0
two.M0 = 0
two.M0 = 0
two.M0 = 0
fpdiv_si atn1, pi_dec, 4
If fpcmp(x, minusone) = 0 Then
fpmul two, two, atn1
fpmul result, two, minusone
Exit Sub
Else
result = pi_dec_half
Exit Sub
End If
End If
fpsub B, one, B '1-x*x
fpsqr B, B 'sqr(1-x*x)
fpdiv term1, T, B
fpatn result, term1
End Sub
Sub fpcos (result As decfloat, z As decfloat)
Dim As decfloat x_2
fpsub x_2, pi_dec_half, z
fpsin result, x_2
End Sub
Sub fptan (result As decfloat, z As decfloat)
Dim As decfloat x_2, s, c
x_2 = z
fpsin s, x_2
x_2 = z
fpcos c, x_2
fpdiv result, s, c
End Sub
Sub fpsin (result As decfloat, x As decfloat)
Dim As decfloat XX, Term, Accum, p, temp2, fac, x_2
Dim As decfloat pi2, circ, Ab
x_2 = x
pi2 = pi_dec
circ = pi2_dec
fpabs Ab, x
If fpcmp(Ab, circ) > 0 Then
'======== CENTRALIZE ==============
'floor/ceil to centralize
Dim As decfloat tmp, tmp2
pi2 = pi2_dec 'got 2*pi
fpdiv tmp, x_2, pi2
tmp2 = tmp
fpfix tmp, tmp 'int part
fpsub tmp, tmp2, tmp 'frac part
fpmul tmp, tmp, pi2
x_2 = tmp
End If
Dim As Long lm, lm2, i
Dim As decfloat factor
lm = 63
lm2 = 1 + Int(-0.45344993886092585968 + 0.022333002852398072433 * lm + 5.0461814408333079844E-7 * lm * lm - 4.2338453039804235772E-11 * lm * lm * lm)
If lm2 < 0 Then lm2 = 0
Si2fp factor, 5
fpipow factor, factor, lm2
fpdiv x_2, x_2, factor 'x_=x_/5^lm2
'==================================
Dim As Long sign(3): sign(3) = 1
Dim As Long c: c = 1
Accum = x_2
Si2fp fac, 1
p = x_2
fpmul XX, x_2, x_2
Do
c = c + 2
temp2 = Accum
fpmul_si fac, fac, c * (c - 1)
fpmul p, p, XX
fpdiv Term, p, fac
If sign(c And 3) Then
fpsub Accum, temp2, Term
Else
fpadd Accum, temp2, Term
End If
Loop Until fpcmp(Accum, temp2) = 0
'multiply the result by 5^lm2
For i = 1 To lm2
fpmul p, Accum, Accum
fpmul temp2, Accum, p
'*** sin(5*x) = 5 * sin(x) - 20 * sin(x)^3 + 16 * sin(x)^5
fpmul_si Accum, Accum, 5
fpmul_si Term, temp2, 20
fpmul_si XX, temp2, 16
fpmul XX, XX, p
fpsub Accum, Accum, Term
fpadd Accum, Accum, XX
Next i
result = Accum
End Sub
Sub fppow (result As decfloat, lhs As decfloat, rhs As decfloat)
Dim As decfloat lhs2
fplog lhs2, lhs
fpmul lhs2, lhs2, rhs
fpexp result, lhs2
End Sub
' sqrt(num)
Sub fpsqr (result As decfloat, num As decfloat)
Dim As decfloat r, r2, tmp, n, half
Dim As Integer ex, k
Dim As String ts, v
Dim As Double x
'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
Si2fp tmp, 0
If fpcmp&(n, tmp) = 0 Then
Si2fp r, 0
result = r
Exit Sub
End If
Si2fp tmp, 1
If fpcmp&(n, tmp) = 0 Then
Si2fp r, 1
result = r
Exit Sub
End If
Si2fp tmp, 0
If fpcmp&(n, tmp) < 0 Then
Si2fp r, 0
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 &H7FFF) - BIAS - 1
Else
ex = 0
End If
ts = _Trim$(Str$(n.M0))
If Len(ts) < 9 Then
ts = ts + String$(9 - Len(ts), "0")
End If
v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
ts = _Trim$(Str$(n.M1))
If Len(ts) < 9 Then
ts = String$(9 - 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
Si2fp r, 1
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, ".")
str2fp r, v
r.exponent = ex \ 2 + BIAS + 1
If Len(v) > 1 And k = 0 Then r.exponent = r.exponent + 1
'half.M0=500000000
'half.exponent = BIAS
'half.sign = 0
str2fp half, ".5"
'=====================================================================
'Newton-Raphson method
For k = 1 To 2
fpdiv tmp, n, r
fpadd r2, r, tmp
fpmul r, r2, half
Next
result = r
End Sub
Sub fpnroot (result As decfloat, x As decfloat, p_in As Long)
Dim As Long p, psign
Dim As decfloat ry, tmp, tmp2
Dim As Double t, t2
Dim As Long i, ex, l
x.sign = 0
psign = Sgn(p_in)
p = Abs(p_in)
l = 2 'Log((NUM_DIGITS + 8) * 0.0625) * 1.5 'calculate the number of iterations needed
ex = (x.exponent And &H7FFF) - BIAS - 1 'extract the exponent
t = x.M0 + x.M1 / 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 / 100000000 '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)
str2fp ry, Str$(t2) 'convert the double t2 to decfloat in ry
t = Log(10) * ex / p
t2 = Exp(t - Fix(t))
str2fp tmp, Str$(t2) 'convert the double t2 to decfloat in tmp
fpmul ry, ry, tmp 'ry = ry * Log(10) * ex / p
str2fp tmp, "2.7182818284590452353602874713527"
fpipow tmp, tmp, Fix(t)
fpmul ry, ry, tmp
fpipow tmp, ry, p - 1 'tmp = ry ^ (p-1)
fpdiv tmp, x, tmp 'tmp = x * tmp
fpmul_si tmp2, ry, p - 1 'tmp2 = ry * (p-1)
fpadd tmp2, tmp2, tmp 'tmp2 = tmp2 + tmp
fpdiv_si ry, tmp2, p 'ry = tmp2 / p
For i = 1 To l '+ 1
fpipow tmp, ry, p - 1 'tmp = ry^(p-1)
fpdiv tmp, x, tmp 'tmp = x/tmp
fpmul_si tmp2, ry, p - 1 'tmp2 = ry*(p-1)
fpadd tmp2, tmp2, tmp 'tmp2 = tmp2+tmp
fpdiv_si ry, tmp2, p 'ry = tmp2/p
Next
If psign < 0 Then
Si2fp tmp, 1
fpdiv ry, tmp, ry
End If
result = ry
End Sub
Sub fpipow (result As decfloat, x As decfloat, e As _Integer64)
'take x to an Long power
Dim As decfloat y
Dim As decfloat z
Dim As _Integer64 n, c
c = 0
y = x
n = Abs(e)
z.sign = 0
z.exponent = (BIAS + 1)
z.M0 = 100000000
While n > 0
While (n And 1) = 0
n = n \ 2
fpmul y, y, y
c = c + 1
Wend
n = n - 1
fpmul z, y, z
c = c + 1
Wend
If e < 0 Then
Si2fp y, 1
fpdiv z, y, z
End If
result = z
End Sub
Sub fpneg (result As decfloat, x As decfloat)
result = x
result.sign = result.sign Xor &H8000
End Sub
Sub fpabs (result As decfloat, x As decfloat)
result = x
result.sign = 0
End Sub
Function fp2str$ (num As decfloat, digits As Integer)
Dim As decfloat n, one
Dim As Long ex, ln
Dim As String s, sd, sign
Dim As _Unsigned Integer xpn
Dim As Integer signn
n = num
xpn = n.exponent
signn = n.sign
If digits < 1 Then digits = 1
If digits > 58 Then digits = 58
'round up
If ndigit(n, digits + 1) > 4 Then
n.exponent = digits + BIAS
n.sign = 0
one.M0 = 100000000
one.exponent = 1 + BIAS
fpadd n, n, one
If n.exponent > (digits + BIAS) Then
xpn = xpn + (n.exponent - (digits + BIAS))
End If
End If
n.exponent = xpn
n.sign = signn
If n.exponent > 0 Then
ex = (n.exponent And &H7FFF) - BIAS - 1
Else
If n.exponent = 0 Then
fp2str$ = " 0"
Exit Function
Else
fp2str$ = "Exponent overflow"
Exit Function
End If
End If
If n.sign Then sign = "-" Else sign = " "
sd = ""
s = _Trim$(Str$(n.M0))
sd = sd + s
s = _Trim$(Str$(n.M1))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
s = _Trim$(Str$(n.M2))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
s = _Trim$(Str$(n.M3))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
s = _Trim$(Str$(n.M4))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
s = _Trim$(Str$(n.M5))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
s = _Trim$(Str$(n.M6))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
sd = Left$(sd, digits)
ln = Len(sd)
If ex >= 0 Then
If ex < digits Then
sd = Left$(sd, ex + 1) + "." + Mid$(sd, ex + 2)
ElseIf ex > digits Then
sd = Left$(sd, 1) + "." + Mid$(sd, 2) + "e" + _Trim$(Str$(ex))
End If
ElseIf ex < 0 Then
If ex > (-5) Then
sd = "." + String$(Abs(ex) - 1, "0") + sd
Else
sd = Left$(sd, 1) + "." + Mid$(sd, 2) + "e" + _Trim$(Str$(ex))
End If
End If
fp2str$ = sign + sd
End Function
Function ndigit& (num As decfloat, digit As Long)
Dim As Long ex, ex2, j, k
Dim As _Unsigned Integer xp
xp = num.exponent
num.exponent = (digit + BIAS)
ex = (num.exponent And &H7FFF) - BIAS
Dim As _Unsigned Long numa(6)
numa(0) = num.M0
numa(1) = num.M1
numa(2) = num.M2
numa(3) = num.M3
numa(4) = num.M4
numa(5) = num.M5
numa(6) = num.M6
ex = ex - 1
ex2 = ex \ 9
k = ex2
j = ex Mod 9
If j = 8 Then
ndigit& = numa(k) Mod 10
Else
ndigit& = (numa(k) \ shift_constants(7 - j)) Mod 10
End If
num.exponent = xp
End Function
Function fp2str_exp$ (n As decfloat, places_in As Long)
Dim As Long ex, places
Dim As String v, f, ts
places = places_in
If places > 54 Then places = 54
places = 62 'places-1
If n.exponent > 0 Then
ex = (n.exponent And &H7FFF) - BIAS - 1
Else
ex = 0
End If
If n.sign Then
v = "-"
Else
v = " "
End If
ts = Str$(n.M0)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = ts + String$(9 - Len(ts), "0")
End If
v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
ts = Str$(n.M1)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
ts = Str$(n.M2)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
ts = Str$(n.M3)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
ts = Str$(n.M4)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
ts = Str$(n.M5)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
ts = Str$(n.M6)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
v = Left$(v, places + 3)
f = _Trim$(Str$(Abs(ex)))
'f = string$(9 - Len(f), "0") + f
If ex < 0 Then
v = v + "e-"
Else
v = v + "e+"
End If
v = v + f
fp2str_exp$ = v
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) > (63 + 1) Then
f1 = Mid$(f1, 1, (63 + 1))
End If
While Len(f1) < (63 + 1)
f1 = f1 + "0"
Wend
j = 1
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M0 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M1 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M2 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M3 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M4 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M5 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M6 = ulng
If ulng <> 0 Then fp = 1
If fp Then n.exponent = (ex + BIAS + 1) Else n.exponent = 0
result = n
End Sub
Sub Si2fp (result As decfloat, m As _Integer64)
Dim As decfloat fac1
Dim As _Integer64 n
n = Abs(m)
If n > 9999999999999999~&& Then
str2fp result, Str$(m)
Exit Sub
End If
fac1.M0 = 0
fac1.M1 = 0
fac1.M2 = 0
fac1.M3 = 0
fac1.M4 = 0
fac1.M5 = 0
fac1.M6 = 0
If m = 0 Then
fac1.exponent = 0
fac1.sign = 0
result = fac1
Exit Sub
End If
fac1.exponent = BIAS
If n < 1000000000 Then
If n < 10 Then
fac1.M0 = n * 100000000
fac1.exponent = fac1.exponent + 1
ElseIf n < 100 Then
fac1.M0 = n * 10000000
fac1.exponent = fac1.exponent + 2
ElseIf n < 1000 Then
fac1.M0 = n * 1000000
fac1.exponent = fac1.exponent + 3
ElseIf n < 10000 Then
fac1.M0 = n * 100000
fac1.exponent = fac1.exponent + 4
ElseIf n < 100000 Then
fac1.M0 = n * 10000
fac1.exponent = fac1.exponent + 5
ElseIf n < 1000000 Then
fac1.M0 = n * 1000
fac1.exponent = fac1.exponent + 6
ElseIf n < 10000000 Then
fac1.M0 = n * 100
fac1.exponent = fac1.exponent + 7
ElseIf n < 100000000 Then
fac1.M0 = n * 10
fac1.exponent = fac1.exponent + 8
ElseIf n < 1000000000 Then
fac1.M0 = n
fac1.exponent = fac1.exponent + 9
End If
End If
If n > 999999999 Then
fac1.exponent = fac1.exponent + 9
If n < 10000000000 Then
fac1.M0 = n \ 10
fac1.M1 = (n Mod 10) * 100000000
fac1.exponent = fac1.exponent + 1
ElseIf n < 100000000000 Then
fac1.M0 = n \ 100
fac1.M1 = (n Mod 100) * 10000000
fac1.exponent = fac1.exponent + 2
ElseIf n < 1000000000000 Then
fac1.M0 = n \ 1000
fac1.M1 = (n Mod 1000) * 1000000
fac1.exponent = fac1.exponent + 3
ElseIf n < 10000000000000 Then
fac1.M0 = n \ 10000
fac1.M1 = (n Mod 10000) * 100000
fac1.exponent = fac1.exponent + 4
ElseIf n < 100000000000000 Then
fac1.M0 = n \ 100000
fac1.M1 = (n Mod 100000) * 10000
fac1.exponent = fac1.exponent + 5
ElseIf n < 1000000000000000 Then
fac1.M0 = n \ 1000000
fac1.M1 = (n Mod 1000000) * 1000
fac1.exponent = fac1.exponent + 6
ElseIf n < 10000000000000000 Then
fac1.M0 = n \ 10000000
fac1.M1 = (n Mod 10000000) * 100
fac1.exponent = fac1.exponent + 7
ElseIf n < 100000000000000000 Then
fac1.M0 = n \ 100000000
fac1.M1 = (n Mod 100000000) * 10
fac1.exponent = fac1.exponent + 8
ElseIf n < 1000000000000000000 Then
fac1.M0 = n \ 1000000000
fac1.M1 = n Mod 1000000000
fac1.exponent = fac1.exponent + 9
End If
End If
If m < 0 Then
fac1.sign = &H8000
Else
fac1.sign = 0
End If
result = fac1
End Sub
Sub RSHIFT_n (mantissa As decfloat, n As Long)
If n = 9 Then
mantissa.M6 = mantissa.M5
mantissa.M5 = mantissa.M4
mantissa.M4 = mantissa.M3
mantissa.M3 = mantissa.M2
mantissa.M2 = mantissa.M1
mantissa.M1 = mantissa.M0
mantissa.M0 = 0
Exit Sub
Else
Dim As _Unsigned Long v1, v2, c1, c2
c1 = shift_constants(n - 1)
c2 = shift_constants(8 - n)
v1 = mantissa.M6 \ c1
v2 = mantissa.M5 Mod c1
v2 = v2 * c2 + v1
mantissa.M6 = v2
v1 = mantissa.M5 \ c1
v2 = mantissa.M4 Mod c1
v2 = v2 * c2 + v1
mantissa.M5 = v2
v1 = mantissa.M4 \ c1
v2 = mantissa.M3 Mod c1
v2 = v2 * c2 + v1
mantissa.M4 = v2
v1 = mantissa.M3 \ c1
v2 = mantissa.M2 Mod c1
v2 = v2 * c2 + v1
mantissa.M3 = v2
v1 = mantissa.M2 \ c1
v2 = mantissa.M1 Mod c1
v2 = v2 * c2 + v1
mantissa.M2 = v2
v1 = mantissa.M1 \ c1
v2 = mantissa.M0 Mod c1
v2 = v2 * c2 + v1
mantissa.M1 = v2
mantissa.M0 = mantissa.M0 \ c1
End If
End Sub
Sub LSHIFT_n (mantissa As decfloat, n As Long)
If n = 9 Then
mantissa.M0 = mantissa.M1
mantissa.M1 = mantissa.M2
mantissa.M2 = mantissa.M3
mantissa.M3 = mantissa.M4
mantissa.M4 = mantissa.M5
mantissa.M5 = mantissa.M6
mantissa.M6 = 0
Exit Sub
Else
Dim As _Unsigned Long v1, v2, c1, c2
c1 = shift_constants(n - 1)
c2 = shift_constants(8 - n)
v1 = mantissa.M0 Mod c2
v2 = mantissa.M1 \ c2
mantissa.M0 = v1 * c1 + v2
mantissa.M1 = mantissa.M1 Mod c2
v1 = mantissa.M1 Mod c2
v2 = mantissa.M2 \ c2
mantissa.M1 = v1 * c1 + v2
mantissa.M2 = mantissa.M2 Mod c2
v1 = mantissa.M2 Mod c2
v2 = mantissa.M3 \ c2
mantissa.M2 = v1 * c1 + v2
mantissa.M3 = mantissa.M3 Mod c2
v1 = mantissa.M3 Mod c2
v2 = mantissa.M4 \ c2
mantissa.M3 = v1 * c1 + v2
mantissa.M4 = mantissa.M4 Mod c2
v1 = mantissa.M4 Mod c2
v2 = mantissa.M5 \ c2
mantissa.M4 = v1 * c1 + v2
mantissa.M5 = mantissa.M5 Mod c2
v1 = mantissa.M5 Mod c2
v2 = mantissa.M6 \ c2
mantissa.M5 = v1 * c1 + v2
mantissa.M6 = mantissa.M6 Mod c2
mantissa.M6 = c1 * (mantissa.M6 Mod c2)
End If
End Sub
Sub LSHIFT_a (mantissa() As _Unsigned Long, n As Long, k As Long)
Dim As _Unsigned Long v1, v2, c1, c2
Dim As Long i
c1 = shift_constants(k - 1)
c2 = shift_constants(8 - k)
For i = 0 To n
v1 = mantissa(i) Mod c2
v2 = mantissa(i + 1) \ c2
mantissa(i) = v1 * c1 + v2
mantissa(i + 1) = mantissa(i + 1) Mod c2
Next
mantissa(n) = c1 * (mantissa(n) Mod c2)
End Sub
Sub RSHIFT_a (mantissa() As _Unsigned Long, n As Long, k As Long)
Dim As Long i
If k = 9 Then
For i = n To 1 Step -1
mantissa(i) = mantissa(i - 1)
Next
mantissa(0) = 0
Exit Sub
Else
Dim As _Unsigned Long v1, v2, c1, c2
c1 = shift_constants(k - 1)
c2 = shift_constants(8 - k)
For i = n To 1 Step -1
v1 = mantissa(i) \ c1
v2 = mantissa(i - 1) Mod c1
v2 = v2 * c2 + v1
mantissa(i) = v2
Next
mantissa(0) = mantissa(0) \ c1
End If
End Sub
Sub LSHIFT_da (mantissa() As Double, n As Long, k As Long)
Dim As _Unsigned Long v1, v2, c1, c2
Dim As Long i
c1 = shift_constants(k - 1)
c2 = shift_constants(3 - k)
For i = 2 To n - 1
v1 = mantissa(i) Mod c2
v2 = mantissa(i + 1) \ c2
mantissa(i) = v1 * c1 + v2
mantissa(i + 1) = Fix(mantissa(i + 1)) Mod c2
Next
mantissa(n) = c1 * (mantissa(n) Mod c2)
End Sub
Function fpcmp& (x As decfloat, y As decfloat)
Dim As Long c
If x.sign < y.sign Then
fpcmp& = -1: Exit Function
End If
If x.sign > y.sign Then
fpcmp& = 1: Exit Function
End If
If x.exponent < y.exponent Then
If x.sign = 0 Then
fpcmp& = -1: Exit Function
Else
fpcmp& = 1: Exit Function
End If
End If
If x.exponent > y.exponent Then
If x.sign = 0 Then
fpcmp& = 1: Exit Function
Else
fpcmp& = -1: Exit Function
End If
End If
c = x.M0 - y.M0
If c <> 0 Then GoTo fpcmpcompare
c = x.M1 - y.M1
If c <> 0 Then GoTo fpcmpcompare
c = x.M2 - y.M2
If c <> 0 Then GoTo fpcmpcompare
c = x.M3 - y.M3
If c <> 0 Then GoTo fpcmpcompare
c = x.M4 - y.M4
If c <> 0 Then GoTo fpcmpcompare
c = x.M5 - y.M5
If c <> 0 Then GoTo fpcmpcompare
c = x.M6 - y.M6
fpcmpcompare:
If c = 0 Then
fpcmp& = 0: Exit Function
End If
If c < 0 Then
If x.sign = 0 Then
fpcmp& = -1: Exit Function
Else
fpcmp& = 1: Exit Function
End If
End If
If c > 0 Then
If x.sign = 0 Then
fpcmp& = 1: Exit Function
Else
fpcmp& = -1: Exit Function
End If
End If
End Function
Sub NORM_FAC1 (fac1 As decfloat)
Dim As Long er, f
' 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.
er = 0: f = 0
If fac1.M0 > 0 Then f = 1
If fac1.M1 > 0 Then f = 1
If fac1.M2 > 0 Then f = 1
If fac1.M3 > 0 Then f = 1
If fac1.M4 > 0 Then f = 1
If fac1.M5 > 0 Then f = 1
If fac1.M6 > 0 Then f = 1
If f = 0 Then
fac1.exponent = 0
Exit Sub
'if the highmost Digit in fac1_man is nonzero,
'shift the mantissa right 1 Digit and
'increment the exponent
ElseIf fac1.M0 > 999999999 Then
RSHIFT_n fac1, 1
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 fac1.M0 = 0
LSHIFT_n fac1, 9
fac1.exponent = fac1.exponent - 9
If fac1.exponent <= 0 Then
Print "NORM_FAC1=EXPU_ERR"
Exit Sub
End If
Wend
If fac1.M0 < 10 Then
LSHIFT_n fac1, 8
fac1.exponent = fac1.exponent - 8
ElseIf fac1.M0 < 100 Then
LSHIFT_n fac1, 7
fac1.exponent = fac1.exponent - 7
ElseIf fac1.M0 < 1000 Then
LSHIFT_n fac1, 6
fac1.exponent = fac1.exponent - 6
ElseIf fac1.M0 < 10000 Then
LSHIFT_n fac1, 5
fac1.exponent = fac1.exponent - 5
ElseIf fac1.M0 < 100000 Then
LSHIFT_n fac1, 4
fac1.exponent = fac1.exponent - 4
ElseIf fac1.M0 < 1000000 Then
LSHIFT_n fac1, 3
fac1.exponent = fac1.exponent - 3
ElseIf fac1.M0 < 10000000 Then
LSHIFT_n fac1, 2
fac1.exponent = fac1.exponent - 2
ElseIf fac1.M0 < 100000000 Then
LSHIFT_n fac1, 1
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, carryover As Long)
Dim As Long v, c
c = carryover
v = fac2.M6 + fac1.M6 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M6 = v
v = fac2.M5 + fac1.M5 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M5 = v
v = fac2.M4 + fac1.M4 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M4 = v
v = fac2.M3 + fac1.M3 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M3 = v
v = fac2.M2 + fac1.M2 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M2 = v
v = fac2.M1 + fac1.M1 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M1 = v
v = fac1.M0 + fac2.M0 + c
fac1.M0 = v
NORM_FAC1 fac1
End Sub
Sub fpsub_aux (fac1 As decfloat, fac2 As decfloat, carryover As Long)
Dim As Long v, c
c = carryover
v = fac1.M6 - fac2.M6 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M6 = v
v = fac1.M5 - fac2.M5 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M5 = v
v = fac1.M4 - fac2.M4 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M4 = v
v = fac1.M3 - fac2.M3 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M3 = v
v = fac1.M2 - fac2.M2 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M2 = v
v = fac1.M1 - fac2.M1 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M1 = v
v = fac1.M0 - fac2.M0 - c
fac1.M0 = v
NORM_FAC1 fac1
End Sub
Sub fpadd (result As decfloat, x As decfloat, y As decfloat)
Dim As decfloat fac1, fac2
Dim As Long i, t, c, xsign, ysign, m
Dim As Unsigned Long fac3(DWORDS)
xsign = x.sign: x.sign = 0
ysign = y.sign: y.sign = 0
c = fpcmp(x, y)
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 &H7FFF) - BIAS - 1) - ((fac2.exponent And &H7FFF) - BIAS - 1)
If t < 63 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
fac3(0) = fac2.M0
fac3(1) = fac2.M1
fac3(2) = fac2.M2
fac3(3) = fac2.M3
fac3(4) = fac2.M4
fac3(5) = fac2.M5
fac3(6) = fac2.M6
If t > 0 And t < 63 Then 'shift
i = t \ 9
While i > 0
RSHIFT_a fac3(), DWORDS, 9
'RSHIFT_n fac2, 9
t = t - 9
i = i - 1
Wend
If t = 8 Then
RSHIFT_a fac3(), DWORDS, 8
'RSHIFT_n fac2, 8
ElseIf t = 7 Then
RSHIFT_a fac3(), DWORDS, 7
'RSHIFT_n fac2, 7
ElseIf t = 6 Then
RSHIFT_a fac3(), DWORDS, 6
'RSHIFT_n fac2, 6
ElseIf t = 5 Then
RSHIFT_a fac3(), DWORDS, 5
'RSHIFT_n fac2, 5
ElseIf t = 4 Then
RSHIFT_a fac3(), DWORDS, 4
'RSHIFT_n fac2, 4
ElseIf t = 3 Then
RSHIFT_a fac3(), DWORDS, 3
'RSHIFT_n fac2, 3
ElseIf t = 2 Then
RSHIFT_a fac3(), DWORDS, 2
'RSHIFT_n fac2, 2
ElseIf t = 1 Then
RSHIFT_a fac3(), DWORDS, 1
'RSHIFT_n fac2, 1
End If
End If
m = 0
If fac3(DWORDS) >= 500000000 Then m = 1
fac2.M0 = fac3(0)
fac2.M1 = fac3(1)
fac2.M2 = fac3(2)
fac2.M3 = fac3(3)
fac2.M4 = fac3(4)
fac2.M5 = fac3(5)
fac2.M6 = fac3(6)
'See if the signs of the two numbers
'are the same. If so, add; if not, subtract.
If fac1.sign = fac2.sign Then 'add
fpadd_aux fac1, fac2, m
Else
fpsub_aux fac1, fac2, m
End If
End If
result = fac1
'now check exponent of result, watch for overflow.
If result.exponent > 32768 Then
'er = EXPO_ERR
Print "add: Exponent over/under flow"
End If
End Sub
Sub fpsub (result As decfloat, x As decfloat, y As decfloat)
Dim As decfloat fac1, fac2
fac1 = x
fac2 = y
fac2.sign = fac2.sign Xor &H8000
fpadd fac1, fac1, fac2
result = fac1
'now check exponent of result, watch for overflow.
If result.exponent > 32768 Then
'er = EXPO_ERR
Print "sub: Exponent over/under flow"
End If
End Sub
Sub fpmul (result As decfloat, x As decfloat, y As decfloat)
'Dim As decfloat fac1,fac2
Dim As Integer i, j, ex, den, num
Dim As _Integer64 digit, carry, prod
Dim As _Unsigned Long fac3(0 To DWORDS * 2 + 1)
Dim As Long fac1(0 To DWORDS - 1), fac2(0 To DWORDS - 1)
Dim As _Unsigned Integer fac1exponent, fac2exponent
Dim As Integer fac1sign, fac2sign
' fac1=x
' fac2=y
'check exponents. if either is zero,
'the result is zero
If x.exponent = 0 Or y.exponent = 0 Then 'result is zero...clear fac1.
result.exponent = 0
result.sign = 0
result.M0 = 0
result.M1 = 0
result.M2 = 0
result.M3 = 0
result.M4 = 0
result.M5 = 0
result.M6 = 0
Exit Sub
Else
fac1(0) = x.M0
fac1(1) = x.M1
fac1(2) = x.M2
fac1(3) = x.M3
fac1(4) = x.M4
fac1(5) = x.M5
fac1(6) = x.M6
fac1exponent = x.exponent
fac1sign = x.sign
fac2(0) = y.M0
fac2(1) = y.M1
fac2(2) = y.M2
fac2(3) = y.M3
fac2(4) = y.M4
fac2(5) = y.M5
fac2(6) = y.M6
fac2exponent = y.exponent
fac2sign = y.sign
'clear fac3 mantissa
For i = 0 To DWORDS * 2 + 1
fac3(i) = 0
Next
den = DWORDS - 1
While fac2(den) = 0
den = den - 1
Wend
num = DWORDS - 1
While fac1(num) = 0
num = num - 1
Wend
If num < den Then
'Swap fac1, fac2
fac2(0) = x.M0
fac2(1) = x.M1
fac2(2) = x.M2
fac2(3) = x.M3
fac2(4) = x.M4
fac2(5) = x.M5
fac2(6) = x.M6
fac2exponent = x.exponent
fac2sign = x.sign
fac1(0) = y.M0
fac1(1) = y.M1
fac1(2) = y.M2
fac1(3) = y.M3
fac1(4) = y.M4
fac1(5) = y.M5
fac1(6) = y.M6
fac1exponent = y.exponent
fac1sign = y.sign
Swap den, num
End If
For j = den To 0 Step -1
carry = 0
digit = fac2(j)
For i = num To 0 Step -1
prod = fac3(i + j + 1) + digit * fac1(i) + carry
carry = prod \ 1000000000
fac3(i + j + 1) = (prod Mod 1000000000)
Next
fac3(j) = carry
Next
j = 0
If fac3(0) < 10~& Then
j = 8
ElseIf fac3(0) < 100~& Then
j = 7
ElseIf fac3(0) < 1000~& Then
j = 6
ElseIf fac3(0) < 10000~& Then
j = 5
ElseIf fac3(0) < 100000~& Then
j = 4
ElseIf fac3(0) < 1000000~& Then
j = 3
ElseIf fac3(0) < 10000000~& Then
j = 2
ElseIf fac3(0) < 100000000~& Then
j = 1
End If
If j > 0 Then
LSHIFT_a fac3(), 7, j
End If
result.M0 = fac3(0)
result.M1 = fac3(1)
result.M2 = fac3(2)
result.M3 = fac3(3)
result.M4 = fac3(4)
result.M5 = fac3(5)
result.M6 = fac3(6)
End If
'now determine exponent of result.
'as you do...watch for overflow.
ex = x.exponent - BIAS + y.exponent
result.exponent = ex - j
If multrndup Then
If fac3(DWORDS) >= 500000000~& Then
Dim As decfloat fac
fac.exponent = ex - j
fac.sign = 0
fac.M0 = 0
fac.M1 = 0
fac.M2 = 0
fac.M3 = 0
fac.M4 = 0
fac.M5 = 0
fac.M6 = 1
fpadd_aux result, fac, 0
End If
End If
'determine the sign of the product
result.sign = x.sign Xor y.sign
NORM_FAC1 result
'now check exponent of result, watch for overflow.
If result.exponent > 32768 Then
'er = EXPO_ERR
Print "mul: Exponent over/under flow"
End If
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_out As decfloat, x As decfloat, y As decfloat)
Dim As Long fac1(6), fac2(6), r
Dim As Long i, er, is_power_of_ten, rndup
Dim As _Unsigned Integer fac1exponent, fac2exponent
Dim As Integer fac1sign, fac2sign
fac1(0) = x.M0
fac1(1) = x.M1
fac1(2) = x.M2
fac1(3) = x.M3
fac1(4) = x.M4
fac1(5) = x.M5
fac1(6) = x.M6
fac1exponent = x.exponent
fac1sign = x.sign
fac2(0) = y.M0
fac2(1) = y.M1
fac2(2) = y.M2
fac2(3) = y.M3
fac2(4) = y.M4
fac2(5) = y.M5
fac2(6) = y.M6
fac2exponent = y.exponent
fac2sign = y.sign
If fac2exponent = 0 Then ' if fac2 = 0, return
' a divide-by-zero error and
' bail out.
result_out.M0 = 999999999
result_out.M1 = 999999999
result_out.M2 = 999999999
result_out.M3 = 999999999
result_out.M4 = 999999999
result_out.M5 = 999999999
result_out.M6 = 999999999
result_out.exponent = 9999 + BIAS + 1
'er = DIVZ_ERR
Print "division by zero"
Exit Sub
ElseIf fac1exponent = 0 Then 'fact1=0, just return
er = 0
result_out.M0 = 0
result_out.M1 = 0
result_out.M2 = 0
result_out.M3 = 0
result_out.M4 = 0
result_out.M5 = 0
result_out.M6 = 0
result_out.exponent = 0
result_out.sign = 0
Exit Sub
Else
'check to see if fac2 is a power of ten
is_power_of_ten = 0
If fac2(0) = 100000000 Then
is_power_of_ten = 1
For i = 1 To 6
If fac2(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
result_out.sign = fac1sign Xor fac2sign
result_out.exponent = fac1exponent - fac2exponent + BIAS + 1
result_out.M0 = fac1(0)
result_out.M1 = fac1(1)
result_out.M2 = fac1(2)
result_out.M3 = fac1(3)
result_out.M4 = fac1(4)
result_out.M5 = fac1(5)
result_out.M6 = fac1(6)
Exit Sub
End If
Const dm = 18
Dim As Double result(1 To dm), n(1 To dm), d(1 To dm)
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 6
' n(2 * j + 2) = fac1(j) \ 10000
' n(2 * j + 3) = fac1(j) Mod 10000
' d(2 * j + 2) = fac2(j) \ 10000
' d(2 * j + 3) = fac2(j) Mod 10000
' Next
n(2) = x.M0 \ 100000: r = x.M0 Mod 100000
n(3) = r \ 10: r = r Mod 10
n(4) = r * 1000 + x.M1 \ 1000000: r = x.M1 Mod 1000000
n(5) = r \ 100: r = r Mod 100
n(6) = r * 100 + x.M2 \ 10000000: r = x.M2 Mod 10000000
n(7) = r \ 1000: r = r Mod 1000
n(8) = r * 10 + x.M3 \ 100000000: r = x.M3 Mod 100000000
n(9) = r \ 10000: r = r Mod 10000
n(10) = r
n(11) = x.M4 \ 100000: r = x.M4 Mod 100000
n(12) = r \ 10: r = r Mod 10
n(13) = r * 1000 + x.M5 \ 1000000: r = x.M5 Mod 1000000
n(14) = r \ 100: r = r Mod 100
n(15) = r * 100 + x.M6 \ 10000000: r = x.M6 Mod 10000000
n(16) = r \ 1000: r = r Mod 1000
n(17) = r * 10
d(2) = y.M0 \ 100000: r = y.M0 Mod 100000
d(3) = r \ 10: r = r Mod 10
d(4) = r * 1000 + y.M1 \ 1000000: r = y.M1 Mod 1000000
d(5) = r \ 100: r = r Mod 100
d(6) = r * 100 + y.M2 \ 10000000: r = y.M2 Mod 10000000
d(7) = r \ 1000: r = r Mod 1000
d(8) = r * 10 + y.M3 \ 100000000: r = y.M3 Mod 100000000
d(9) = r \ 10000: r = r Mod 10000
d(10) = r
d(11) = y.M4 \ 100000: r = y.M4 Mod 100000
d(12) = r \ 10: r = r Mod 10
d(13) = r * 1000 + y.M5 \ 1000000: r = y.M5 Mod 1000000
d(14) = r \ 100: r = r Mod 100
d(15) = r * 100 + y.M6 \ 10000000: r = y.M6 Mod 10000000
d(16) = r \ 1000: r = r Mod 1000
d(17) = r * 10
n(1) = (fac1exponent And &H7FFF) - BIAS - 1
d(1) = (fac2exponent And &H7FFF) - 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))
subtract w(), q, d(), (stp + 2), last
normalize w(), (stp + 2), q
Next
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)
j = 0
If result(2) < 10~& Then
j = 3
ElseIf result(2) < 100~& Then
j = 2
ElseIf result(2) < 1000~& Then
j = 1
End If
If j > 0 Then
LSHIFT_da result(), dm, j
End If
j = 2
result_out.M0 = (result(j) * 10000 + result(j + 1)) * 10 + result(j + 2) \ 1000: r = result(j + 2) Mod 1000: j = 5
result_out.M1 = (r * 10000 + result(j)) * 100 + result(j + 1) \ 100: r = result(j + 1) Mod 100: j = 7
result_out.M2 = (r * 10000 + result(j)) * 1000 + result(j + 1) \ 10: r = result(j + 1) Mod 10: j = 9
result_out.M3 = (r * 10000 + result(j)) * 10000 + result(j + 1): j = 11
result_out.M4 = (result(j) * 10000 + result(j + 1)) * 10 + result(j + 2) \ 1000: r = result(j + 2) Mod 1000: j = 14
result_out.M5 = (r * 10000 + result(j)) * 100 + result(j + 1) \ 100: r = result(j + 1) Mod 100: j = 16
result_out.M6 = (r * 10000 + result(j)) * 1000 + result(j + 1) \ 10
rndup = 0
If (result(j + 1) Mod 10) > 4 Then
Dim As decfloat ru
ru.M0 = 100000000
ru.exponent = 1 + BIAS
result_out.exponent = 63 + BIAS
fpadd result_out, result_out, ru
If result_out.exponent > (63 + BIAS) Then
rndup = 1
End If
End If
NORM_FAC1 result_out
result_out.exponent = (result(1) + BIAS) + rndup
End If
result_out.sign = fac1sign Xor fac2sign
End Sub
Sub fpdiv_si (result As decfloat, num As decfloat, den As Long)
Dim As Long fac1(6)
Dim As _Unsigned Integer fac1exponent
Dim As Integer fac1sign
Dim As _Unsigned _Integer64 carry, remder
Dim As _Integer64 i, divisor
Dim As _Integer64 quotient
remder = 0
divisor = Abs(den)
fac1(0) = num.M0
fac1(1) = num.M1
fac1(2) = num.M2
fac1(3) = num.M3
fac1(4) = num.M4
fac1(5) = num.M5
fac1(6) = num.M6
fac1exponent = num.exponent
fac1sign = num.sign
result.M0 = num.M0
result.M1 = num.M1
result.M2 = num.M2
result.M3 = num.M3
result.M4 = num.M4
result.M5 = num.M5
result.M6 = num.M6
result.exponent = num.exponent
result.sign = num.sign
If divisor = 0 Then
Print "error: divisor = 0"
Exit Sub
End If
If divisor > 2147483647 Then
Print "error: divisor too large"
Exit Sub
End If
For i = 0 To 6
quotient = fac1(i) + remder * 1000000000
remder = quotient Mod divisor
fac1(i) = quotient \ divisor
Next
quotient = remder * 1000000000
quotient = quotient \ divisor
carry = fac1(0)
result.M0 = fac1(0)
result.M1 = fac1(1)
result.M2 = fac1(2)
result.M3 = fac1(3)
result.M4 = fac1(4)
result.M5 = fac1(5)
result.M6 = fac1(6)
result.exponent = fac1exponent
result.sign = fac1sign
If carry = 0 Then
LSHIFT_n result, 9
result.exponent = result.exponent - 9
result.M6 = result.M6 + quotient
ElseIf carry < 10 Then
LSHIFT_n result, 8
result.exponent = result.exponent - 8
result.M6 = result.M6 + quotient \ 10
ElseIf carry < 100 Then
LSHIFT_n result, 7
result.exponent = result.exponent - 7
result.M6 = result.M6 + quotient \ 100
ElseIf carry < 1000 Then
LSHIFT_n result, 6
result.exponent = result.exponent - 6
result.M6 = result.M6 + quotient \ 1000
ElseIf carry < 10000 Then
LSHIFT_n result, 5
result.exponent = result.exponent - 5
result.M6 = result.M6 + quotient \ 10000
ElseIf carry < 100000 Then
LSHIFT_n result, 4
result.exponent = result.exponent - 4
result.M6 = result.M6 + quotient \ 100000
ElseIf carry < 1000000 Then
LSHIFT_n result, 3
result.exponent = result.exponent - 3
result.M6 = result.M6 + quotient \ 1000000
ElseIf carry < 10000000 Then
LSHIFT_n result, 2
result.exponent = result.exponent - 2
result.M6 = result.M6 + quotient \ 10000000
ElseIf carry < 100000000 Then
LSHIFT_n result, 1
result.exponent = result.exponent - 1
result.M6 = result.M6 + quotient \ 100000000
End If
'NORM_FAC1(fac1)
result.sign = fac1sign
If den < 0 Then
result.sign = fac1sign Xor &H8000
End If
'now check exponent of result, watch for overflow.
If result.exponent > 32768 Then
'er = EXPO_ERR
Print "div_si: Exponent over/under flow"
End If
End Sub
Sub fpmul_si (result As decfloat, x As decfloat, y As _Integer64)
Dim As decfloat fac1, fac2
Dim As Long i
Dim As _Integer64 carry, digit, prod, value
Dim As _Unsigned Long fac3(7)
fac1 = x
digit = Abs(y)
If digit > 999999999 Then
Si2fp fac2, y
fpmul fac1, fac1, fac2
result = fac1: Exit Sub
End If
If fac1.exponent = 0 Or y = 0 Then 'result is zero...clear fac1.
fac1.sign = 0
fac1.exponent = 0
fac1.M0 = 0
fac1.M1 = 0
fac1.M2 = 0
fac1.M3 = 0
fac1.M4 = 0
fac1.M5 = 0
fac1.M6 = 0
NORM_FAC1 fac1
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
fac3(0) = fac1.M0
fac3(1) = fac1.M1
fac3(2) = fac1.M2
fac3(3) = fac1.M3
fac3(4) = fac1.M4
fac3(5) = fac1.M5
fac3(6) = fac1.M6
fac3(7) = 0
carry = 0
For i = 6 To 0 Step -1
prod = digit * fac3(i) + carry
value = prod Mod 1000000000
fac3(i) = value
carry = prod \ 1000000000
Next
If carry < 10 Then
RSHIFT_a fac3(), 7, 1
fac1.exponent = fac1.exponent + 1
fac3(0) = fac3(0) + carry * 100000000
ElseIf carry < 100 Then
RSHIFT_a fac3(), 7, 2
fac1.exponent = fac1.exponent + 2
fac3(0) = fac3(0) + carry * 10000000
ElseIf carry < 1000 Then
RSHIFT_a fac3(), 7, 3
fac1.exponent = fac1.exponent + 3
fac3(0) = fac3(0) + carry * 1000000
ElseIf carry < 10000 Then
RSHIFT_a fac3(), 7, 4
fac1.exponent = fac1.exponent + 4
fac3(0) = fac3(0) + carry * 100000
ElseIf carry < 100000 Then
RSHIFT_a fac3(), 7, 5
fac1.exponent = fac1.exponent + 5
fac3(0) = fac3(0) + carry * 10000
ElseIf carry < 1000000 Then
RSHIFT_a fac3(), 7, 6
fac1.exponent = fac1.exponent + 6
fac3(0) = fac3(0) + carry * 1000
ElseIf carry < 10000000 Then
RSHIFT_a fac3(), 7, 7
fac1.exponent = fac1.exponent + 7
fac3(0) = fac3(0) + carry * 100
ElseIf carry < 100000000 Then
RSHIFT_a fac3(), 7, 8
fac1.exponent = fac1.exponent + 8
fac3(0) = fac3(0) + carry * 10
ElseIf carry < 1000000000 Then
RSHIFT_a fac3(), 7, 9
fac1.exponent = fac1.exponent + 9
fac3(0) = fac3(0) + carry
End If
End If
fac1.M0 = fac3(0)
fac1.M1 = fac3(1)
fac1.M2 = fac3(2)
fac1.M3 = fac3(3)
fac1.M4 = fac3(4)
fac1.M5 = fac3(5)
fac1.M6 = fac3(6)
If fac3(7) >= 500000000~& Then
Dim As decfloat fac
fac.exponent = fac1.exponent
fac.sign = 0
fac.M0 = 0
fac.M1 = 0
fac.M2 = 0
fac.M3 = 0
fac.M4 = 0
fac.M5 = 0
fac.M6 = 1
fpadd_aux fac1, fac, 0
End If
NORM_FAC1 fac1
If y < 0 Then
fac1.sign = fac1.sign Xor &H8000
End If
result = fac1
'now check exponent of result, watch for overflow.
If result.exponent > 32768 Then
'er = EXPO_ERR
Print "mul_si: Exponent over/under flow"
End If
End Sub
'integer part of num
Sub fpfix (result As decfloat, num As decfloat)
Dim As decfloat ips
Dim As Long ex, ex2, j, k
ex = (num.exponent And &H7FFF) - BIAS
If ex < 1 Then
result = ips
Exit Sub
End If
If ex >= 63 Then
result = num
Exit Sub
End If
Dim As _Unsigned Long ip(6), numa(6)
numa(0) = num.M0
numa(1) = num.M1
numa(2) = num.M2
numa(3) = num.M3
numa(4) = num.M4
numa(5) = num.M5
numa(6) = num.M6
ex2 = ex \ 9
k = ex2
j = ex Mod 9
While ex2 > 0
ex2 = ex2 - 1
ip(ex2) = numa(ex2)
Wend
If j = 1 Then
ip(k) = 100000000 * (numa(k) \ 100000000)
ElseIf j = 2 Then
ip(k) = 10000000 * (numa(k) \ 10000000)
ElseIf j = 3 Then
ip(k) = 1000000 * (numa(k) \ 1000000)
ElseIf j = 4 Then
ip(k) = 100000 * (numa(k) \ 100000)
ElseIf j = 5 Then
ip(k) = 10000 * (numa(k) \ 10000)
ElseIf j = 6 Then
ip(k) = 1000 * (numa(k) \ 1000)
ElseIf j = 7 Then
ip(k) = 100 * (numa(k) \ 100)
ElseIf j = 8 Then
ip(k) = 10 * (numa(k) \ 10)
ElseIf j = 9 Then
ip(k) = numa(k)
End If
result.exponent = ex + BIAS
result.sign = num.sign
result.M0 = ip(0)
result.M1 = ip(1)
result.M2 = ip(2)
result.M3 = ip(3)
result.M4 = ip(4)
result.M5 = ip(5)
result.M6 = ip(6)
End Sub
Sub fpfrac (result As decfloat, num As decfloat)
Dim As decfloat ip, fp
fpfix ip, num
fpsub fp, num, ip
result = fp
End Sub
Function fpfix_is_odd& (num As decfloat)
Dim As Long ex, ex2, j, k
ex = (num.exponent And &H7FFF) - BIAS
If ex < 0 Then
fpfix_is_odd = 0
Exit Function
End If
If ex >= 63 Then
fpfix_is_odd = 0
Exit Function
End If
Dim As _Unsigned Long numa(6)
numa(0) = num.M0
numa(1) = num.M1
numa(2) = num.M2
numa(3) = num.M3
numa(4) = num.M4
numa(5) = num.M5
numa(6) = num.M6
ex = ex - 1
ex2 = ex \ 9
k = ex2
j = ex Mod 9
If j = 0 Then
fpfix_is_odd = (numa(k) \ 100000000) And 1: Exit Function
ElseIf j = 1 Then
fpfix_is_odd = (numa(k) \ 10000000) And 1: Exit Function
ElseIf j = 2 Then
fpfix_is_odd = (numa(k) \ 1000000) And 1: Exit Function
ElseIf j = 3 Then
fpfix_is_odd = (numa(k) \ 100000) And 1: Exit Function
ElseIf j = 4 Then
fpfix_is_odd = (numa(k) \ 10000) And 1: Exit Function
ElseIf j = 5 Then
fpfix_is_odd = (numa(k) \ 1000) And 1: Exit Function
ElseIf j = 6 Then
fpfix_is_odd = (numa(k) \ 100) And 1: Exit Function
ElseIf j = 7 Then
fpfix_is_odd = (numa(k) \ 10) And 1: Exit Function
ElseIf j = 8 Then
fpfix_is_odd = numa(k) And 1: Exit Function
End If
fpfix_is_odd = 0
End Function
Posts: 2,186
Threads: 222
Joined: Apr 2022
Reputation:
104
On my ad/subtract/multiply/divide calculator, I just kept all results as fractions, and displayed the results as a decimal. It works great for simple operations. Very easy to go in one direction, back in the other, and end up with same number you started with. However, I discovered that method is daunting when it comes to dealing with logs, decimal powers and decimal roots. Rounding is the alternative, but I haven't found any definitive source to build a reliable routine around. It's weird, but I suppose some results that look like rounding to 2 would be correct, but how about in a case where 1.99999999999999999902 is the correct answer? So we know we can't just make rounding rules based on assumptions, and more importantly, for a high precision calculator to be accepted, it would have to conform to some apparent "universal" rounding standards, whatever the heck they are. Unfortunately, my last name isn't Casio. Have you found any authoritative source to work from?
Pete
Posts: 422
Threads: 27
Joined: Apr 2022
Reputation:
26
10-12-2022, 01:43 AM
(This post was last modified: 10-12-2022, 01:54 AM by Jack.)
Pete, at the moment the conversion routine fp2str_exp is set to print 63 digits, it's so that I can spot the accuracy of the calculations
but the fp2str routine is set to display a maximum of 58 digits and it will correctly roundup the number.
BTW, with this code, setting the optimization to -O3 makes a significant speed improvement
Posts: 2,186
Threads: 222
Joined: Apr 2022
Reputation:
104
Ah, so you do have a reliable rounding method. kudos!
Pete
Shoot first and shoot people who ask questions, later.
Posts: 422
Threads: 27
Joined: Apr 2022
Reputation:
26
10-14-2022, 12:04 AM
(This post was last modified: 10-14-2022, 08:49 PM by Jack.)
decfloat with a super-simple eval, at the "?" prompt enter a numerical expresion, example
1-1/3!+1/5!-1/7!+1/9!-1/11!+1/13!-1/15!+1/17!-1/19!+1/21!-1/23!+1/25!-1/27!+1/29!-1/31!+1/33!-1/35!+1/37!-1/39!+1/41!-1/43!+1/45!-1/47!+1/49!
result .8414709848078965066525023216302989996225630607983710656728
P.S. added root(x,n) function, example: root(2,3) --> 1.2599210498948731647672106072782...
Code: (Select All) _Title "defloat-32 v2 by Jack"
$NoPrefix
$Console:Only
Dest Console
Dim Shared As Long shift_constants(7)
shift_constants(0) = 10
shift_constants(1) = 100
shift_constants(2) = 1000
shift_constants(3) = 10000
shift_constants(4) = 100000
shift_constants(5) = 1000000
shift_constants(6) = 10000000
shift_constants(7) = 100000000
Const BIAS = 16384 '2 ^ 14
Const DWORDS = 7
Type decfloat
As _Unsigned Integer exponent
As Integer sign
As _Unsigned Long M0
As _Unsigned Long M1
As _Unsigned Long M2
As _Unsigned Long M3
As _Unsigned Long M4
As _Unsigned Long M5
As _Unsigned Long M6
End Type
Dim Shared As decfloat pi_dec, pi2_dec, pi_dec_half, xpn, xpd
Dim Shared As Long multrndup
multrndup = 1
str2fp pi_dec, "3.14159265358979323846264338327950288419716939937510582097494459"
str2fp pi2_dec, "6.28318530717958647692528676655900576839433879875021164194988918"
str2fp pi_dec_half, "1.57079632679489661923132169163975144209858469968755291048747230"
' Error definitions
Const DIVZ_ERR = 1 'Divide by zero
Const EXPO_ERR = 2 'Exponent overflow error
Const EXPU_ERR = 3 'Exponent underflow error
'op_stack used by eval
ReDim Shared op_stack(10) As String * 6
Dim Shared As Long op_stack_pointer, op_stack_ubound
op_stack_pointer = 0: op_stack_ubound = 10
'value_stack used by eval
ReDim Shared value_stack(10) As decfloat
Dim Shared As Long value_stack_pointer, value_stack_ubound
value_stack_pointer = 0: value_stack_ubound = 10
'variables used by eval
Dim Shared As Integer expr_index, expr_length
Dim Shared As String op, number, expression, char, rpn
Dim Shared As decfloat factorial, lhs_value, rhs_value
'variables used by demo
Dim As String math_expr
Dim As decfloat result_value
' eval demo
math_expr = " "
While math_expr <> ""
Line Input "? ", math_expr
eval math_expr, result_value
Print fp2str(result_value, 58)
Print "rpn = "; rpn
Wend
' end of demo
'================================
Sub op_stack_push (n As String)
op_stack(op_stack_pointer) = n
If op_stack_pointer = op_stack_ubound Then
op_stack_ubound = op_stack_ubound + 10
ReDim _Preserve As String * 6 op_stack(op_stack_ubound)
End If
op_stack_pointer = op_stack_pointer + 1
End Sub
Function op_stack_pop$ ()
If op_stack_pointer = 0 Then
Print "op_stack is empty"
op_stack_pop = ""
End If
op_stack_pointer = op_stack_pointer - 1
op_stack_pop = op_stack(op_stack_pointer)
End Function
Sub value_stack_push (n As decfloat)
value_stack(value_stack_pointer) = n
If value_stack_pointer = value_stack_ubound Then
value_stack_ubound = value_stack_ubound + 10
ReDim _Preserve As decfloat value_stack(value_stack_ubound)
End If
value_stack_pointer = value_stack_pointer + 1
End Sub
Sub value_stack_pop (result As decfloat)
If value_stack_pointer = 0 Then
Print "value_stack is empty"
Exit Sub
End If
value_stack_pointer = value_stack_pointer - 1
result = value_stack(value_stack_pointer)
End Sub
'eval
Sub eval (exprn As String, result As decfloat)
expression = UCase$(exprn)
If Len(expression) = 0 Then expression = "0"
rpn = ""
expr_index = 1: expr_length = Len(expression)
value_stack_pointer = 0
op_stack_pointer = 0
scan
expr
If char <> " " Then
Print
Print "Syntax Error"
Print
End If
Print
value_stack_pop result 'value_stack ( 0 )
End Sub
Sub scan
If expr_index > expr_length Then
char = " "
Exit Sub
End If
char = Mid$(expression, expr_index, 1)
expr_index = expr_index + 1
If char = " " Then scan
End Sub
Sub unary
Dim As decfloat tmp
If char = "-" Or char = "+" Then
op_stack_push char
scan
term
op = op_stack_pop
If _Trim$(op) <> "-" Then Exit Sub
value_stack_pop lhs_value
tmp = lhs_value
tmp.sign = tmp.sign Xor &H8000
value_stack_push tmp
rpn = rpn + "NEG "
Exit Sub
End If
factor
End Sub
Sub gamma
unary
While char = "!"
value_stack_pop lhs_value
fn_factorial
value_stack_push factorial
scan
rpn = rpn + "! "
Wend
End Sub
Sub expon
Dim As decfloat tmp
gamma
While char = "^"
scan
gamma
value_stack_pop rhs_value
value_stack_pop lhs_value
fpfrac tmp, rhs_value
If tmp.exponent = 0 Then
fpfix tmp, rhs_value
fpipow tmp, lhs_value, fp2long(tmp)
Else
fplog tmp, lhs_value
fpmul tmp, tmp, rhs_value
fpexp tmp, tmp
End If
value_stack_push tmp
rpn = rpn + "^ "
Wend
End Sub
Sub term
Dim As decfloat tmp
expon
While (char = "*" Or char = "/")
op_stack_push char
scan
expon
op = op_stack_pop
If _Trim$(op) = "*" Then
value_stack_pop rhs_value
value_stack_pop lhs_value
fpmul tmp, lhs_value, rhs_value
value_stack_push tmp
rpn = rpn + "* "
End If
If _Trim$(op) = "/" Then
value_stack_pop rhs_value
value_stack_pop lhs_value
fpdiv tmp, lhs_value, rhs_value
value_stack_push tmp
rpn = rpn + "/ "
End If
Wend
End Sub
Sub expr
Dim As decfloat tmp
term
While (char = "-" Or char = "+")
op_stack_push char
scan
term
op = op_stack_pop
If _Trim$(op) = "-" Then
value_stack_pop rhs_value
value_stack_pop lhs_value
fpsub tmp, lhs_value, rhs_value
value_stack_push tmp
rpn = rpn + "- "
End If
If _Trim$(op) = "+" Then
value_stack_pop rhs_value
value_stack_pop lhs_value
fpadd tmp, lhs_value, rhs_value
value_stack_push tmp
rpn = rpn + "+ "
End If
Wend
End Sub
Sub factor
Dim As decfloat tmp
If InStr(".0123456789", char) Then
number = ""
While char <> "." And InStr("0123456789", char)
number = number + char
scan
Wend
If char = "." Then
number = number + char
scan
End If
While InStr("0123456789", char)
number = number + char
scan
Wend
If char = "E" Then
number = number + char
scan
If char = "-" Or char = "+" Then
number = number + char
scan
End If
If InStr("0123456789", char) Then
While InStr("0123456789", char)
number = number + char
scan
Wend
Else
number = number + "0"
End If
End If
str2fp tmp, number
value_stack_push tmp
rpn = rpn + number + " "
Exit Sub
End If
If char = "(" Then
scan
expr
If char = "," Then
op_stack_push char
scan
expr
op = op_stack_pop
End If
If char = ")" Then
scan
Else
Print
Print "Missing ')'"
End If
End If
'functions would be added here
'here's a very crude example just to illustrate
If char = "S" Then
If Mid$(expression, expr_index - 1, 4) = "SQR(" Then
expr_index = expr_index + 2 'advance pointer to just before "("
scan
factor
If _Trim$(op) = "," Then
Print "too many arguments to function Sqr"
value_stack_pointer = value_stack_pointer - 1
Else
fpsqr tmp, value_stack(value_stack_pointer - 1)
value_stack(value_stack_pointer - 1) = tmp
rpn = rpn + "SQR "
End If
End If
ElseIf char = "R" Then '' check for more functions
If Mid$(expression, expr_index - 1, 5) = "ROOT(" Then
expr_index = expr_index + 3
scan
factor
If _Trim$(op) = "," Then
value_stack_pointer = value_stack_pointer - 1
fpnroot value_stack(value_stack_pointer - 1), value_stack(value_stack_pointer - 1), fp2long(value_stack(value_stack_pointer))
rpn = rpn + "ROOT "
op = ""
Else
Print "missing second argument in function Root"
End If
End If
End If
End Sub
Function fp2long& (n As decfloat)
Dim As Long ex, ip, s
ex = (n.exponent And &H7FFF) - BIAS - 1
If ex = 0 Then
ip = n.M0 \ 100000000
ElseIf ex = 1 Then
ip = n.M0 \ 10000000
ElseIf ex = 2 Then
ip = n.M0 \ 1000000
ElseIf ex = 3 Then
ip = n.M0 \ 100000
ElseIf ex = 4 Then
ip = n.M0 \ 10000
ElseIf ex = 5 Then
ip = n.M0 \ 1000
ElseIf ex = 6 Then
ip = n.M0 \ 100
ElseIf ex = 7 Then
ip = n.M0 \ 10
ElseIf ex = 8 Then
ip = n.M0
End If
s = n.sign
If s <> 0 Then ip = -ip
fp2long = ip
End Function
Sub fn_factorial
Dim As Long i
Si2fp factorial, 1
For i = 1 To fp2long(lhs_value)
fpmul_si factorial, factorial, i
Next i
End Sub
Sub fpexp_cf (result As decfloat, z As decfloat)
Dim As decfloat x, s, x1, x2, tmp, tmp2
Dim As Long ex, i, ip, n
Static As decfloat one, euler
n = 10
If euler.exponent = 0 Then
str2fp euler, "2.71828182845904523536028747135266249775724709369995957496696763"
one.sign = 0
one.exponent = (BIAS + 1)
one.M0 = 100000000
one.M1 = 0
one.M2 = 0
one.M3 = 0
one.M4 = 0
one.M5 = 0
one.M6 = 0
End If
' 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)
If z.exponent = 0 Then
result = one
Exit Sub
Else
x = z
fpfix tmp, x
ex = (tmp.exponent And &H7FFF) - BIAS - 1
If ex = 0 Then
ip = tmp.M0 \ 100000000
ElseIf ex = 1 Then
ip = tmp.M0 \ 10000000
ElseIf ex = 2 Then
ip = tmp.M0 \ 1000000
ElseIf ex = 3 Then
ip = tmp.M0 \ 100000
ElseIf ex = 4 Then
ip = tmp.M0 \ 10000
ElseIf ex = 5 Then
ip = tmp.M0 \ 1000
ElseIf ex = 6 Then
ip = tmp.M0 \ 100
ElseIf ex = 7 Then
ip = tmp.M0 \ 10
ElseIf ex = 8 Then
ip = tmp.M0
End If
fpfrac x, x
fpdiv_si x1, x, 16 '32768 ' 2^13
fpmul x2, x1, x1
Si2fp s, 4 * n + 6
For i = n To 0 Step -1
Si2fp tmp, 4 * i + 2
fpdiv s, x2, s
fpadd s, s, tmp
Next
fpdiv s, x1, s
Si2fp tmp, 1
tmp2 = tmp
fpadd tmp, tmp, s
fpsub tmp2, tmp2, s
fpdiv s, tmp, tmp2
For i = 1 To 4 '15 ' s^13
fpmul s, s, s
Next
fpipow tmp, euler, ip
fpmul result, s, tmp
End If
End Sub
Sub fpexp (result As decfloat, z As decfloat)
Dim As decfloat x, xpd, xpn
Dim As Long ex, i, ip
Static As decfloat one, euler, exponentialn(17), exponentiald(17)
If exponentialn(0).exponent = 0 Then
Si2fp exponentialn(0), 1
str2fp exponentialn(1), ".503787763085141830426316836835640413180555373602013787765506507"
str2fp exponentialn(2), ".123113754811115580240199391979941449439554252687431414889261917"
str2fp exponentialn(3), "0.194014786172100868397379733957773068817466004597117759916094387e-1"
str2fp exponentialn(4), "0.221053514053621305631679326515226926261161756634139701475334846e-2"
str2fp exponentialn(5), "0.193454865206744920585608214573628003845560417585149728091593248e-3"
str2fp exponentialn(6), "0.134817206102830919302898061303157717441762386334371916554674483e-4"
str2fp exponentialn(7), "7.65161710091203441510338717292779761121416735865848016278269163e-7"
str2fp exponentialn(8), "3.58548435070740017294732467393249181951311755864779732794287718e-8"
str2fp exponentialn(9), "1.39715881029815224422817956869312878658017689986058208796945882e-9"
str2fp exponentialn(10), "4.53464900275857744649901664225234895410940311133637948725625465e-11"
str2fp exponentialn(11), "1.22102800612904874865986466356114289547925418884289566733570583e-12"
str2fp exponentialn(12), "2.69942119086103755038141670655438202935697783465198117242288620e-14"
str2fp exponentialn(13), "4.80709158088050415577996944372246571546305491654413869758830148e-16"
str2fp exponentialn(14), "6.67446218117009199081966862645625972417483845805581243134455550e-18"
str2fp exponentialn(15), "6.82626745056408867477242403742610973442937848360726389682334564e-20"
str2fp exponentialn(16), "4.60485274414644615989466638049841829236889098591288845615263752e-22"
str2fp exponentialn(17), "1.54768708764447533766179612797166908030142040749901411508788783e-24"
Si2fp exponentiald(0), 1
str2fp exponentiald(1), "-.496212236914858169573683163164359586819444626397986212234493493"
str2fp exponentiald(2), ".119325991725973749813882555144301036258998879085417627123755410"
str2fp exponentiald(3), "-0.184850613180012448539696668330106026341966320933794116815659028e-1"
str2fp exponentiald(4), "0.206797341469361126562570972007228490388291451667636384019144944e-2"
str2fp exponentiald(5), "-0.177476640029445154798895481195212175295325133388353948891798456e-3"
str2fp exponentiald(6), "0.121119697791530382246392392593837579456791651154395855360677559e-4"
str2fp exponentiald(7), "-6.72136000621469478770704356909331255584862479068643438675858294e-7"
str2fp exponentiald(8), "3.07422101412642477250911345710710651015389275273741476097162835e-8"
str2fp exponentiald(9), "-1.16701132447665083402247823313377310891548863251014843319140333e-9"
str2fp exponentiald(10), "3.68187090359605484256890447350797602918878757385065579230141815e-11"
str2fp exponentiald(11), "-9.61338101757297931016053534655678826462590722321717850532044698e-13"
str2fp exponentiald(12), "2.05509049962341068403280550130851518395842965866907693627957626e-14"
str2fp exponentiald(13), "-3.52746564887524616342404030739724544161545203402406781919425177e-16"
str2fp exponentiald(14), "4.70347116684618761386137496578607769970491970387898680511311425e-18"
str2fp exponentiald(15), "-4.59998458732742996655211581291170107043066864993354030929388480e-20"
str2fp exponentiald(16), "2.95256671673728810659480905519132991941158386504354019400480885e-22"
str2fp exponentiald(17), "-9.38719670297727932196732935370899807766127816686343783729700392e-25"
str2fp euler, "2.71828182845904523536028747135266249775724709369995957496696763"
one.sign = 0
one.exponent = (BIAS + 1)
one.M0 = 100000000
one.M1 = 0
one.M2 = 0
one.M3 = 0
one.M4 = 0
one.M5 = 0
one.M6 = 0
End If
If z.exponent = 0 Then
result = one
Exit Sub
Else
x = z
fpfix xpn, x
ex = (xpn.exponent And &H7FFF) - BIAS - 1
If ex = 0 Then
ip = xpn.M0 \ 100000000
ElseIf ex = 1 Then
ip = xpn.M0 \ 10000000
ElseIf ex = 2 Then
ip = xpn.M0 \ 1000000
ElseIf ex = 3 Then
ip = xpn.M0 \ 100000
ElseIf ex = 4 Then
ip = xpn.M0 \ 10000
ElseIf ex = 5 Then
ip = xpn.M0 \ 1000
ElseIf ex = 6 Then
ip = xpn.M0 \ 100
ElseIf ex = 7 Then
ip = xpn.M0 \ 10
ElseIf ex = 8 Then
ip = xpn.M0
End If
fpfrac x, x
fpmul xpn, exponentialn(17), x
For i = 16 To 0 Step -1
fpadd xpn, xpn, exponentialn(i)
fpmul xpn, xpn, x
Next
fpmul xpd, exponentiald(17), x
For i = 16 To 0 Step -1
fpadd xpd, xpd, exponentiald(i)
fpmul xpd, xpd, x
Next
fpdiv xpn, xpn, xpd
fpipow xpd, euler, ip
fpmul result, xpn, xpd
End If
End Sub
Sub fplog (result As decfloat, z As decfloat)
Dim As decfloat x, xpd, xpn
Dim As Long ex, i
Static As decfloat fpln10, one, logarithmn(22), logarithmd(24)
If logarithmn(0).exponent = 0 Then
Si2fp logarithmn(0), 1
str2fp logarithmn(1), "10.2310520037359891285036755418469258827514540594420472830660753"
str2fp logarithmn(2), "48.6627011559331606160440612789640717584111988048005196804661728"
str2fp logarithmn(3), "142.883839025492099568442001118949943253750601849543916975383500"
str2fp logarithmn(4), "290.042754824482811137859430225980858456412815030133662556507508"
str2fp logarithmn(5), "431.988540379436291654167001685688951066167323915161481314938683"
str2fp logarithmn(6), "488.988532243997467759781233533239006273973991246082219125463314"
str2fp logarithmn(7), "429.964620903161621905726323536985462923884859411935413887634588"
str2fp logarithmn(8), "297.667076072840840028860885452874923891756123601208220339342769"
str2fp logarithmn(9), "163.496246139169410681932569760493099792107706750131034678156585"
str2fp logarithmn(10), "71.4689010653151414969354268941968334288999401577488778627880921"
str2fp logarithmn(11), "24.8454129712280032841853780591175947041028424111335032188809328"
str2fp logarithmn(12), "6.83955915425638170286810308584881831872981167696540451581896838"
str2fp logarithmn(13), "1.47913861498769669443157158341616767519966697200814961672841044"
str2fp logarithmn(14), ".248298871909766565290925597836167509193223133304509724924096267"
str2fp logarithmn(15), "0.318133082830696859180749013242710465363599980360532595003499144e-1"
str2fp logarithmn(16), "0.304020922309269050902540112467577787426750224236109865380535113e-2"
str2fp logarithmn(17), "0.209978185892809569953723710672234372599610486457723316925724789e-3"
str2fp logarithmn(18), "0.100311469452988931955197299915204806738105549379720367796921701e-4"
str2fp logarithmn(19), "3.11027381043420800967758659489427925931042335875083315411300509e-7"
str2fp logarithmn(20), "5.67148960215418390372633356121445540918492594737567360786766486e-9"
str2fp logarithmn(21), "5.12133417069955237659371135908115173603073989914114172007865451e-11"
str2fp logarithmn(22), "1.55966557818637624120176613242050057417746497195348156692036175e-13"
Si2fp logarithmd(0), 1
str2fp logarithmd(1), "10.7310520037359891285036755418469258827514540594420472830660738"
str2fp logarithmd(2), "53.6948938244678218469625657165542013664535925011882099886674915"
str2fp logarithmd(3), "166.404268603147347449088725463278068642726913413657339541359027"
str2fp logarithmd(4), "357.829354185501208195542189937563557126312937751426774299756481"
str2fp logarithmd(5), "566.879307659840870904948427820813871913477928966397711410863037"
str2fp logarithmd(6), "685.659474589076294970895474467750778199315371532367251137728620"
str2fp logarithmd(7), "647.552215345138140652347527859581928707051002114597464686830520"
str2fp logarithmd(8), "484.337597665677116149955858143764741744286855179474277002221246"
str2fp logarithmd(9), "289.339023581167087024317711351001662002419858544547812026145386"
str2fp logarithmd(10), "138.626741697453253940901869606998742571923585752594114076259249"
str2fp logarithmd(11), "53.2967984633656431994464301417545208285735718517685211650739340"
str2fp logarithmd(12), "16.3977970191605190152676712715525604280407449895571363588426976"
str2fp logarithmd(13), "4.01344869231909159178580772978304524142055416630126527640586347"
str2fp logarithmd(14), ".774140683278008004255775706220863488435679762354478168050645942"
str2fp logarithmd(15), ".116116794765397457080324911252454433970856146870870944104060064"
str2fp logarithmd(16), "0.132995134740536235985194287064770583314074985150347963060323352e-1"
str2fp logarithmd(17), "0.113495871890204061922407499675091077013725252521536694012719530e-2"
str2fp logarithmd(18), "0.697980393999574131276345270647881218831498205731172278517042043e-4"
str2fp logarithmd(19), "0.295301241701723056394124042476367677312456922839751324605065650e-5"
str2fp logarithmd(20), "8.03411378004236615745112193270089533877927829488633592988420379e-8"
str2fp logarithmd(21), "1.26487301511055907218913028822339261374104341757811651124170618e-9"
str2fp logarithmd(22), "9.55430735765744135751947902497685515801773089075935394013525279e-12"
str2fp logarithmd(23), "2.22906425441965268150705233939758319540146232049504069112205822e-14"
str2fp logarithmd(24), "-3.94659794323870669763249350331131008555188524939917666031562314e-18"
str2fp fpln10, "2.30258509299404568401799145468436420760110148862877297603332790"
one.sign = 0
one.exponent = (BIAS + 1)
one.M0 = 100000000
one.M1 = 0
one.M2 = 0
one.M3 = 0
one.M4 = 0
one.M5 = 0
one.M6 = 0
End If
x = z
If fpcmp&(x, one) = 0 Then
result.sign = 0
result.exponent = 0
result.M0 = 0
result.M1 = 0
result.M2 = 0
result.M3 = 0
result.M4 = 0
result.M5 = 0
result.M6 = 0
Exit Sub
End If
If x.exponent <> 0 Then
ex = (x.exponent And &H7FFF) - BIAS - 1
x.exponent = BIAS + 1
fpsqr x, x
fpsqr x, x
fpsub x, x, one
fpmul xpn, logarithmn(22), x
For i = 21 To 0 Step -1
fpadd xpn, xpn, logarithmn(i)
fpmul xpn, xpn, x
Next
fpmul xpn, xpn, x
fpmul xpd, logarithmd(24), x
For i = 23 To 0 Step -1
fpadd xpd, xpd, logarithmd(i)
fpmul xpd, xpd, x
Next
fpdiv xpn, xpn, xpd
fpmul_si xpn, xpn, 4
fpmul_si xpd, fpln10, ex
fpadd result, xpn, xpd
End If
End Sub
Sub fpatn (result As decfloat, x As decfloat)
Dim As Long sign(3): sign(3) = 1
Dim As Long z, c: c = 1
Dim As decfloat XX, Termn, Accum, strC, x_2, mt, mt2, p
Dim As decfloat decnum, one, decnum2, factorn
decnum2 = x
decnum2.sign = 0
one.sign = 0
one.exponent = (BIAS + 1)
one.M0 = 100000000
one.M1 = 0
one.M2 = 0
one.M3 = 0
one.M4 = 0
one.M5 = 0
one.M6 = 0
If fpcmp(decnum2, one) = 0 Then
fpdiv_si decnum, pi_dec, 4
decnum.sign = x.sign
result = decnum
Exit Sub
End If
decnum2.sign = x.sign
Dim As Long lm2: lm2 = 16
Si2fp factorn, _ShL(2, (lm2 - 1))
For z = 1 To lm2
fpmul decnum, decnum2, decnum2
fpadd decnum, decnum, one
fpsqr decnum, decnum
fpadd decnum, decnum, one
fpdiv decnum, decnum2, decnum
decnum2 = decnum
Next z
mt = decnum
x_2 = decnum
p = decnum
fpmul XX, x_2, x_2
Do
c = c + 2
mt2 = mt
Si2fp strC, c
fpmul p, p, XX
fpdiv Termn, p, strC
If sign(c And 3) Then
fpsub Accum, mt, Termn
Else
fpadd Accum, mt, Termn
End If
Swap mt, Accum
Loop Until fpcmp(mt, mt2) = 0
fpmul result, factorn, mt
End Sub
Sub fpasin (result As decfloat, x As decfloat)
Dim As Double num
Dim As decfloat one, T, B, term1, minusone
' ARCSIN = ATN(x / SQR(-x * x + 1))
'============= ARCSIN GUARD =========
num = Val(fp2str(x, 16))
If num > 1 Then
result = one
Exit Sub
End If
If num < -1 Then
result = one
Exit Sub
End If
'========================
one.sign = 0
one.exponent = (BIAS + 1)
one.M0 = 100000000
one.M1 = 0
one.M2 = 0
one.M3 = 0
one.M4 = 0
one.M5 = 0
one.M6 = 0
minusone = one: minusone.sign = &H8000
T = x
fpmul B, x, x 'x*x
'for 1 and -1
If fpcmp(B, one) = 0 Then
Dim As decfloat two, atn1
two.sign = 0
two.exponent = (BIAS + 1)
two.M0 = 200000000
two.M0 = 0
two.M0 = 0
two.M0 = 0
two.M0 = 0
two.M0 = 0
two.M0 = 0
fpdiv_si atn1, pi_dec, 4
If fpcmp(x, minusone) = 0 Then
fpmul two, two, atn1
fpmul result, two, minusone
Exit Sub
Else
result = pi_dec_half
Exit Sub
End If
End If
fpsub B, one, B '1-x*x
fpsqr B, B 'sqr(1-x*x)
fpdiv term1, T, B
fpatn result, term1
End Sub
Sub fpcos (result As decfloat, z As decfloat)
Dim As decfloat x_2
fpsub x_2, pi_dec_half, z
fpsin result, x_2
End Sub
Sub fptan (result As decfloat, z As decfloat)
Dim As decfloat x_2, s, c
x_2 = z
fpsin s, x_2
x_2 = z
fpcos c, x_2
fpdiv result, s, c
End Sub
Sub fpsin (result As decfloat, x As decfloat)
Dim As decfloat XX, Termn, Accum, p, temp2, fac, x_2
Dim As decfloat pi2, circ, Ab
x_2 = x
pi2 = pi_dec
circ = pi2_dec
fpabs Ab, x
If fpcmp(Ab, circ) > 0 Then
'======== CENTRALIZE ==============
'floor/ceil to centralize
Dim As decfloat tmp, tmp2
pi2 = pi2_dec 'got 2*pi
fpdiv tmp, x_2, pi2
tmp2 = tmp
fpfix tmp, tmp 'int part
fpsub tmp, tmp2, tmp 'frac part
fpmul tmp, tmp, pi2
x_2 = tmp
End If
Dim As Long lm, lm2, i
Dim As decfloat factorn
lm = 63
lm2 = 1 + Int(-0.45344993886092585968 + 0.022333002852398072433 * lm + 5.0461814408333079844E-7 * lm * lm - 4.2338453039804235772E-11 * lm * lm * lm)
If lm2 < 0 Then lm2 = 0
Si2fp factorn, 5
fpipow factorn, factorn, lm2
fpdiv x_2, x_2, factorn 'x_=x_/5^lm2
'==================================
Dim As Long sign(3): sign(3) = 1
Dim As Long c: c = 1
Accum = x_2
Si2fp fac, 1
p = x_2
fpmul XX, x_2, x_2
Do
c = c + 2
temp2 = Accum
fpmul_si fac, fac, c * (c - 1)
fpmul p, p, XX
fpdiv Termn, p, fac
If sign(c And 3) Then
fpsub Accum, temp2, Termn
Else
fpadd Accum, temp2, Termn
End If
Loop Until fpcmp(Accum, temp2) = 0
'multiply the result by 5^lm2
For i = 1 To lm2
fpmul p, Accum, Accum
fpmul temp2, Accum, p
'*** sin(5*x) = 5 * sin(x) - 20 * sin(x)^3 + 16 * sin(x)^5
fpmul_si Accum, Accum, 5
fpmul_si Termn, temp2, 20
fpmul_si XX, temp2, 16
fpmul XX, XX, p
fpsub Accum, Accum, Termn
fpadd Accum, Accum, XX
Next i
result = Accum
End Sub
Sub fppow (result As decfloat, lhs As decfloat, rhs As decfloat)
Dim As decfloat lhs2
fplog lhs2, lhs
fpmul lhs2, lhs2, rhs
fpexp result, lhs2
End Sub
' sqrt(num)
Sub fpsqr (result As decfloat, num As decfloat)
Dim As decfloat r, r2, tmp, n, half
Dim As Integer ex, k
Dim As String ts, v
Dim As Double x
'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
Si2fp tmp, 0
If fpcmp&(n, tmp) = 0 Then
Si2fp r, 0
result = r
Exit Sub
End If
Si2fp tmp, 1
If fpcmp&(n, tmp) = 0 Then
Si2fp r, 1
result = r
Exit Sub
End If
Si2fp tmp, 0
If fpcmp&(n, tmp) < 0 Then
Si2fp r, 0
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 &H7FFF) - BIAS - 1
Else
ex = 0
End If
ts = _Trim$(Str$(n.M0))
If Len(ts) < 9 Then
ts = ts + String$(9 - Len(ts), "0")
End If
v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
ts = _Trim$(Str$(n.M1))
If Len(ts) < 9 Then
ts = String$(9 - 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
Si2fp r, 1
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, ".")
str2fp r, v
r.exponent = ex \ 2 + BIAS + 1
If Len(v) > 1 And k = 0 Then r.exponent = r.exponent + 1
'half.M0=500000000
'half.exponent = BIAS
'half.sign = 0
str2fp half, ".5"
'=====================================================================
'Newton-Raphson method
For k = 1 To 2
fpdiv tmp, n, r
fpadd r2, r, tmp
fpmul r, r2, half
Next
result = r
End Sub
Sub fpnroot (result As decfloat, x_in As decfloat, p_in As Long)
Dim As Long p, psign
Dim As decfloat ry, tmp, tmp2, x
Dim As Double t, t2
Dim As Long i, ex, l, s
x = x_in
If x.sign <> 0 And (p_in And 1) = 0 Then
Print "can't extract root of negative number to an even power"
Exit Sub
End If
s = x.sign
x.sign = 0
psign = Sgn(p_in)
p = Abs(p_in)
l = 2 'Log((NUM_DIGITS + 8) * 0.0625) * 1.5 'calculate the number of iterations needed
ex = (x.exponent And &H7FFF) - BIAS - 1 'extract the exponent
t = x.M0 + x.M1 / 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 / 100000000 '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)
str2fp ry, Str$(t2) 'convert the double t2 to decfloat in ry
t = Log(10) * ex / p
t2 = Exp(t - Fix(t))
str2fp tmp, Str$(t2) 'convert the double t2 to decfloat in tmp
fpmul ry, ry, tmp 'ry = ry * Log(10) * ex / p
str2fp tmp, "2.7182818284590452353602874713527"
fpipow tmp, tmp, Fix(t)
fpmul ry, ry, tmp
fpipow tmp, ry, p - 1 'tmp = ry ^ (p-1)
fpdiv tmp, x, tmp 'tmp = x * tmp
fpmul_si tmp2, ry, p - 1 'tmp2 = ry * (p-1)
fpadd tmp2, tmp2, tmp 'tmp2 = tmp2 + tmp
fpdiv_si ry, tmp2, p 'ry = tmp2 / p
For i = 1 To l '+ 1
fpipow tmp, ry, p - 1 'tmp = ry^(p-1)
fpdiv tmp, x, tmp 'tmp = x/tmp
fpmul_si tmp2, ry, p - 1 'tmp2 = ry*(p-1)
fpadd tmp2, tmp2, tmp 'tmp2 = tmp2+tmp
fpdiv_si ry, tmp2, p 'ry = tmp2/p
Next
If psign < 0 Then
Si2fp tmp, 1
fpdiv ry, tmp, ry
End If
result = ry
If s <> 0 Then
result.sign = &H8000
End If
End Sub
Sub fpipow (result As decfloat, x As decfloat, e As _Integer64)
'take x to an Long power
Dim As decfloat y
Dim As decfloat z
Dim As _Integer64 n, c
c = 0
y = x
n = Abs(e)
z.sign = 0
z.exponent = (BIAS + 1)
z.M0 = 100000000
While n > 0
While (n And 1) = 0
n = n \ 2
fpmul y, y, y
c = c + 1
Wend
n = n - 1
fpmul z, y, z
c = c + 1
Wend
If e < 0 Then
Si2fp y, 1
fpdiv z, y, z
End If
result = z
End Sub
Sub fpneg (result As decfloat, x As decfloat)
result = x
result.sign = result.sign Xor &H8000
End Sub
Sub fpabs (result As decfloat, x As decfloat)
result = x
result.sign = 0
End Sub
Function fp2str$ (num As decfloat, digits As Integer)
Dim As decfloat n, one
Dim As Long ex, ln
Dim As String s, sd, sign
Dim As _Unsigned Integer xpn
Dim As Integer signn
n = num
xpn = n.exponent
signn = n.sign
If digits < 1 Then digits = 1
If digits > 58 Then digits = 58
'round up
If ndigit(n, digits + 1) > 4 Then
n.exponent = digits + BIAS
n.sign = 0
one.M0 = 100000000
one.exponent = 1 + BIAS
fpadd n, n, one
If n.exponent > (digits + BIAS) Then
xpn = xpn + (n.exponent - (digits + BIAS))
End If
End If
n.exponent = xpn
n.sign = signn
If n.exponent > 0 Then
ex = (n.exponent And &H7FFF) - BIAS - 1
Else
If n.exponent = 0 Then
fp2str$ = " 0"
Exit Function
Else
fp2str$ = "Exponent overflow"
Exit Function
End If
End If
If n.sign Then sign = "-" Else sign = " "
sd = ""
s = _Trim$(Str$(n.M0))
sd = sd + s
s = _Trim$(Str$(n.M1))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
s = _Trim$(Str$(n.M2))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
s = _Trim$(Str$(n.M3))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
s = _Trim$(Str$(n.M4))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
s = _Trim$(Str$(n.M5))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
s = _Trim$(Str$(n.M6))
ln = Len(s)
If ln < 9 Then
s = String$(9 - ln, "0") + s
End If
sd = sd + s
sd = Left$(sd, digits)
ln = Len(sd)
If ex >= 0 Then
If ex < digits Then
sd = Left$(sd, ex + 1) + "." + Mid$(sd, ex + 2)
ElseIf ex > digits Then
sd = Left$(sd, 1) + "." + Mid$(sd, 2) + "e" + _Trim$(Str$(ex))
End If
ElseIf ex < 0 Then
If ex > (-5) Then
sd = "." + String$(Abs(ex) - 1, "0") + sd
Else
sd = Left$(sd, 1) + "." + Mid$(sd, 2) + "e" + _Trim$(Str$(ex))
End If
End If
fp2str$ = sign + sd
End Function
Function ndigit& (num As decfloat, digit As Long)
Dim As Long ex, ex2, j, k
Dim As _Unsigned Integer xp
xp = num.exponent
num.exponent = (digit + BIAS)
ex = (num.exponent And &H7FFF) - BIAS
Dim As _Unsigned Long numa(6)
numa(0) = num.M0
numa(1) = num.M1
numa(2) = num.M2
numa(3) = num.M3
numa(4) = num.M4
numa(5) = num.M5
numa(6) = num.M6
ex = ex - 1
ex2 = ex \ 9
k = ex2
j = ex Mod 9
If j = 8 Then
ndigit& = numa(k) Mod 10
Else
ndigit& = (numa(k) \ shift_constants(7 - j)) Mod 10
End If
num.exponent = xp
End Function
Function fp2str_exp$ (n As decfloat, places_in As Long)
Dim As Long ex, places
Dim As String v, f, ts
places = places_in
If places > 54 Then places = 54
places = 62 'places-1
If n.exponent > 0 Then
ex = (n.exponent And &H7FFF) - BIAS - 1
Else
ex = 0
End If
If n.sign Then
v = "-"
Else
v = " "
End If
ts = Str$(n.M0)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = ts + String$(9 - Len(ts), "0")
End If
v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
ts = Str$(n.M1)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
ts = Str$(n.M2)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
ts = Str$(n.M3)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
ts = Str$(n.M4)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
ts = Str$(n.M5)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
ts = Str$(n.M6)
ts = _Trim$(ts)
If Len(ts) < 9 Then
ts = String$(9 - Len(ts), "0") + ts
End If
v = v + ts
v = Left$(v, places + 3)
f = _Trim$(Str$(Abs(ex)))
'f = string$(9 - Len(f), "0") + f
If ex < 0 Then
v = v + "e-"
Else
v = v + "e+"
End If
v = v + f
fp2str_exp$ = v
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) > (63 + 1) Then
f1 = Mid$(f1, 1, (63 + 1))
End If
While Len(f1) < (63 + 1)
f1 = f1 + "0"
Wend
j = 1
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M0 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M1 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M2 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M3 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M4 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M5 = ulng
If ulng <> 0 Then fp = 1
j = j + 9
ts = Mid$(f1, j, 9)
ulng = Val(ts)
n.M6 = ulng
If ulng <> 0 Then fp = 1
If fp Then n.exponent = (ex + BIAS + 1) Else n.exponent = 0
result = n
End Sub
Sub Si2fp (result As decfloat, m As _Integer64)
Dim As decfloat fac1
Dim As _Integer64 n
n = Abs(m)
If n > 9999999999999999~&& Then
str2fp result, Str$(m)
Exit Sub
End If
fac1.M0 = 0
fac1.M1 = 0
fac1.M2 = 0
fac1.M3 = 0
fac1.M4 = 0
fac1.M5 = 0
fac1.M6 = 0
If m = 0 Then
fac1.exponent = 0
fac1.sign = 0
result = fac1
Exit Sub
End If
fac1.exponent = BIAS
If n < 1000000000 Then
If n < 10 Then
fac1.M0 = n * 100000000
fac1.exponent = fac1.exponent + 1
ElseIf n < 100 Then
fac1.M0 = n * 10000000
fac1.exponent = fac1.exponent + 2
ElseIf n < 1000 Then
fac1.M0 = n * 1000000
fac1.exponent = fac1.exponent + 3
ElseIf n < 10000 Then
fac1.M0 = n * 100000
fac1.exponent = fac1.exponent + 4
ElseIf n < 100000 Then
fac1.M0 = n * 10000
fac1.exponent = fac1.exponent + 5
ElseIf n < 1000000 Then
fac1.M0 = n * 1000
fac1.exponent = fac1.exponent + 6
ElseIf n < 10000000 Then
fac1.M0 = n * 100
fac1.exponent = fac1.exponent + 7
ElseIf n < 100000000 Then
fac1.M0 = n * 10
fac1.exponent = fac1.exponent + 8
ElseIf n < 1000000000 Then
fac1.M0 = n
fac1.exponent = fac1.exponent + 9
End If
End If
If n > 999999999 Then
fac1.exponent = fac1.exponent + 9
If n < 10000000000 Then
fac1.M0 = n \ 10
fac1.M1 = (n Mod 10) * 100000000
fac1.exponent = fac1.exponent + 1
ElseIf n < 100000000000 Then
fac1.M0 = n \ 100
fac1.M1 = (n Mod 100) * 10000000
fac1.exponent = fac1.exponent + 2
ElseIf n < 1000000000000 Then
fac1.M0 = n \ 1000
fac1.M1 = (n Mod 1000) * 1000000
fac1.exponent = fac1.exponent + 3
ElseIf n < 10000000000000 Then
fac1.M0 = n \ 10000
fac1.M1 = (n Mod 10000) * 100000
fac1.exponent = fac1.exponent + 4
ElseIf n < 100000000000000 Then
fac1.M0 = n \ 100000
fac1.M1 = (n Mod 100000) * 10000
fac1.exponent = fac1.exponent + 5
ElseIf n < 1000000000000000 Then
fac1.M0 = n \ 1000000
fac1.M1 = (n Mod 1000000) * 1000
fac1.exponent = fac1.exponent + 6
ElseIf n < 10000000000000000 Then
fac1.M0 = n \ 10000000
fac1.M1 = (n Mod 10000000) * 100
fac1.exponent = fac1.exponent + 7
ElseIf n < 100000000000000000 Then
fac1.M0 = n \ 100000000
fac1.M1 = (n Mod 100000000) * 10
fac1.exponent = fac1.exponent + 8
ElseIf n < 1000000000000000000 Then
fac1.M0 = n \ 1000000000
fac1.M1 = n Mod 1000000000
fac1.exponent = fac1.exponent + 9
End If
End If
If m < 0 Then
fac1.sign = &H8000
Else
fac1.sign = 0
End If
result = fac1
End Sub
Sub RSHIFT_n (mantissa As decfloat, n As Long)
If n = 9 Then
mantissa.M6 = mantissa.M5
mantissa.M5 = mantissa.M4
mantissa.M4 = mantissa.M3
mantissa.M3 = mantissa.M2
mantissa.M2 = mantissa.M1
mantissa.M1 = mantissa.M0
mantissa.M0 = 0
Exit Sub
Else
Dim As _Unsigned Long v1, v2, c1, c2
c1 = shift_constants(n - 1)
c2 = shift_constants(8 - n)
v1 = mantissa.M6 \ c1
v2 = mantissa.M5 Mod c1
v2 = v2 * c2 + v1
mantissa.M6 = v2
v1 = mantissa.M5 \ c1
v2 = mantissa.M4 Mod c1
v2 = v2 * c2 + v1
mantissa.M5 = v2
v1 = mantissa.M4 \ c1
v2 = mantissa.M3 Mod c1
v2 = v2 * c2 + v1
mantissa.M4 = v2
v1 = mantissa.M3 \ c1
v2 = mantissa.M2 Mod c1
v2 = v2 * c2 + v1
mantissa.M3 = v2
v1 = mantissa.M2 \ c1
v2 = mantissa.M1 Mod c1
v2 = v2 * c2 + v1
mantissa.M2 = v2
v1 = mantissa.M1 \ c1
v2 = mantissa.M0 Mod c1
v2 = v2 * c2 + v1
mantissa.M1 = v2
mantissa.M0 = mantissa.M0 \ c1
End If
End Sub
Sub LSHIFT_n (mantissa As decfloat, n As Long)
If n = 9 Then
mantissa.M0 = mantissa.M1
mantissa.M1 = mantissa.M2
mantissa.M2 = mantissa.M3
mantissa.M3 = mantissa.M4
mantissa.M4 = mantissa.M5
mantissa.M5 = mantissa.M6
mantissa.M6 = 0
Exit Sub
Else
Dim As _Unsigned Long v1, v2, c1, c2
c1 = shift_constants(n - 1)
c2 = shift_constants(8 - n)
v1 = mantissa.M0 Mod c2
v2 = mantissa.M1 \ c2
mantissa.M0 = v1 * c1 + v2
mantissa.M1 = mantissa.M1 Mod c2
v1 = mantissa.M1 Mod c2
v2 = mantissa.M2 \ c2
mantissa.M1 = v1 * c1 + v2
mantissa.M2 = mantissa.M2 Mod c2
v1 = mantissa.M2 Mod c2
v2 = mantissa.M3 \ c2
mantissa.M2 = v1 * c1 + v2
mantissa.M3 = mantissa.M3 Mod c2
v1 = mantissa.M3 Mod c2
v2 = mantissa.M4 \ c2
mantissa.M3 = v1 * c1 + v2
mantissa.M4 = mantissa.M4 Mod c2
v1 = mantissa.M4 Mod c2
v2 = mantissa.M5 \ c2
mantissa.M4 = v1 * c1 + v2
mantissa.M5 = mantissa.M5 Mod c2
v1 = mantissa.M5 Mod c2
v2 = mantissa.M6 \ c2
mantissa.M5 = v1 * c1 + v2
mantissa.M6 = mantissa.M6 Mod c2
mantissa.M6 = c1 * (mantissa.M6 Mod c2)
End If
End Sub
Sub LSHIFT_a (mantissa() As _Unsigned Long, n As Long, k As Long)
Dim As _Unsigned Long v1, v2, c1, c2
Dim As Long i
c1 = shift_constants(k - 1)
c2 = shift_constants(8 - k)
For i = 0 To n
v1 = mantissa(i) Mod c2
v2 = mantissa(i + 1) \ c2
mantissa(i) = v1 * c1 + v2
mantissa(i + 1) = mantissa(i + 1) Mod c2
Next
mantissa(n) = c1 * (mantissa(n) Mod c2)
End Sub
Sub RSHIFT_a (mantissa() As _Unsigned Long, n As Long, k As Long)
Dim As Long i
If k = 9 Then
For i = n To 1 Step -1
mantissa(i) = mantissa(i - 1)
Next
mantissa(0) = 0
Exit Sub
Else
Dim As _Unsigned Long v1, v2, c1, c2
c1 = shift_constants(k - 1)
c2 = shift_constants(8 - k)
For i = n To 1 Step -1
v1 = mantissa(i) \ c1
v2 = mantissa(i - 1) Mod c1
v2 = v2 * c2 + v1
mantissa(i) = v2
Next
mantissa(0) = mantissa(0) \ c1
End If
End Sub
Sub LSHIFT_da (mantissa() As Double, n As Long, k As Long)
Dim As _Unsigned Long v1, v2, c1, c2
Dim As Long i
c1 = shift_constants(k - 1)
c2 = shift_constants(3 - k)
For i = 2 To n - 1
v1 = mantissa(i) Mod c2
v2 = mantissa(i + 1) \ c2
mantissa(i) = v1 * c1 + v2
mantissa(i + 1) = Fix(mantissa(i + 1)) Mod c2
Next
mantissa(n) = c1 * (mantissa(n) Mod c2)
End Sub
Function fpcmp& (x As decfloat, y As decfloat)
Dim As Long c
If x.sign < y.sign Then
fpcmp& = -1: Exit Function
End If
If x.sign > y.sign Then
fpcmp& = 1: Exit Function
End If
If x.exponent < y.exponent Then
If x.sign = 0 Then
fpcmp& = -1: Exit Function
Else
fpcmp& = 1: Exit Function
End If
End If
If x.exponent > y.exponent Then
If x.sign = 0 Then
fpcmp& = 1: Exit Function
Else
fpcmp& = -1: Exit Function
End If
End If
c = x.M0 - y.M0
If c <> 0 Then GoTo fpcmpcompare
c = x.M1 - y.M1
If c <> 0 Then GoTo fpcmpcompare
c = x.M2 - y.M2
If c <> 0 Then GoTo fpcmpcompare
c = x.M3 - y.M3
If c <> 0 Then GoTo fpcmpcompare
c = x.M4 - y.M4
If c <> 0 Then GoTo fpcmpcompare
c = x.M5 - y.M5
If c <> 0 Then GoTo fpcmpcompare
c = x.M6 - y.M6
fpcmpcompare:
If c = 0 Then
fpcmp& = 0: Exit Function
End If
If c < 0 Then
If x.sign = 0 Then
fpcmp& = -1: Exit Function
Else
fpcmp& = 1: Exit Function
End If
End If
If c > 0 Then
If x.sign = 0 Then
fpcmp& = 1: Exit Function
Else
fpcmp& = -1: Exit Function
End If
End If
End Function
Sub NORM_FAC1 (fac1 As decfloat)
Dim As Long er, f
' 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.
er = 0: f = 0
If fac1.M0 > 0 Then f = 1
If fac1.M1 > 0 Then f = 1
If fac1.M2 > 0 Then f = 1
If fac1.M3 > 0 Then f = 1
If fac1.M4 > 0 Then f = 1
If fac1.M5 > 0 Then f = 1
If fac1.M6 > 0 Then f = 1
If f = 0 Then
fac1.exponent = 0
Exit Sub
'if the highmost Digit in fac1_man is nonzero,
'shift the mantissa right 1 Digit and
'increment the exponent
ElseIf fac1.M0 > 999999999 Then
RSHIFT_n fac1, 1
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 fac1.M0 = 0
LSHIFT_n fac1, 9
fac1.exponent = fac1.exponent - 9
If fac1.exponent <= 0 Then
Print "NORM_FAC1=EXPU_ERR"
Exit Sub
End If
Wend
If fac1.M0 < 10 Then
LSHIFT_n fac1, 8
fac1.exponent = fac1.exponent - 8
ElseIf fac1.M0 < 100 Then
LSHIFT_n fac1, 7
fac1.exponent = fac1.exponent - 7
ElseIf fac1.M0 < 1000 Then
LSHIFT_n fac1, 6
fac1.exponent = fac1.exponent - 6
ElseIf fac1.M0 < 10000 Then
LSHIFT_n fac1, 5
fac1.exponent = fac1.exponent - 5
ElseIf fac1.M0 < 100000 Then
LSHIFT_n fac1, 4
fac1.exponent = fac1.exponent - 4
ElseIf fac1.M0 < 1000000 Then
LSHIFT_n fac1, 3
fac1.exponent = fac1.exponent - 3
ElseIf fac1.M0 < 10000000 Then
LSHIFT_n fac1, 2
fac1.exponent = fac1.exponent - 2
ElseIf fac1.M0 < 100000000 Then
LSHIFT_n fac1, 1
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, carryover As Long)
Dim As Long v, c
c = carryover
v = fac2.M6 + fac1.M6 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M6 = v
v = fac2.M5 + fac1.M5 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M5 = v
v = fac2.M4 + fac1.M4 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M4 = v
v = fac2.M3 + fac1.M3 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M3 = v
v = fac2.M2 + fac1.M2 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M2 = v
v = fac2.M1 + fac1.M1 + c
If v > 999999999 Then
v = v - 1000000000
c = 1
Else
c = 0
End If
fac1.M1 = v
v = fac1.M0 + fac2.M0 + c
fac1.M0 = v
NORM_FAC1 fac1
End Sub
Sub fpsub_aux (fac1 As decfloat, fac2 As decfloat, carryover As Long)
Dim As Long v, c
c = carryover
v = fac1.M6 - fac2.M6 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M6 = v
v = fac1.M5 - fac2.M5 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M5 = v
v = fac1.M4 - fac2.M4 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M4 = v
v = fac1.M3 - fac2.M3 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M3 = v
v = fac1.M2 - fac2.M2 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M2 = v
v = fac1.M1 - fac2.M1 - c
If v < 0 Then
v = v + 1000000000
c = 1
Else
c = 0
End If
fac1.M1 = v
v = fac1.M0 - fac2.M0 - c
fac1.M0 = v
NORM_FAC1 fac1
End Sub
Sub fpadd (result As decfloat, x As decfloat, y As decfloat)
Dim As decfloat fac1, fac2
Dim As Long i, t, c, xsign, ysign, m
Dim As Unsigned Long fac3(DWORDS)
xsign = x.sign: x.sign = 0
ysign = y.sign: y.sign = 0
c = fpcmp(x, y)
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 &H7FFF) - BIAS - 1) - ((fac2.exponent And &H7FFF) - BIAS - 1)
If t < 63 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
fac3(0) = fac2.M0
fac3(1) = fac2.M1
fac3(2) = fac2.M2
fac3(3) = fac2.M3
fac3(4) = fac2.M4
fac3(5) = fac2.M5
fac3(6) = fac2.M6
If t > 0 And t < 63 Then 'shift
i = t \ 9
While i > 0
RSHIFT_a fac3(), DWORDS, 9
'RSHIFT_n fac2, 9
t = t - 9
i = i - 1
Wend
If t = 8 Then
RSHIFT_a fac3(), DWORDS, 8
'RSHIFT_n fac2, 8
ElseIf t = 7 Then
RSHIFT_a fac3(), DWORDS, 7
'RSHIFT_n fac2, 7
ElseIf t = 6 Then
RSHIFT_a fac3(), DWORDS, 6
'RSHIFT_n fac2, 6
ElseIf t = 5 Then
RSHIFT_a fac3(), DWORDS, 5
'RSHIFT_n fac2, 5
ElseIf t = 4 Then
RSHIFT_a fac3(), DWORDS, 4
'RSHIFT_n fac2, 4
ElseIf t = 3 Then
RSHIFT_a fac3(), DWORDS, 3
'RSHIFT_n fac2, 3
ElseIf t = 2 Then
RSHIFT_a fac3(), DWORDS, 2
'RSHIFT_n fac2, 2
ElseIf t = 1 Then
RSHIFT_a fac3(), DWORDS, 1
'RSHIFT_n fac2, 1
End If
End If
m = 0
If fac3(DWORDS) >= 500000000 Then m = 1
fac2.M0 = fac3(0)
fac2.M1 = fac3(1)
fac2.M2 = fac3(2)
fac2.M3 = fac3(3)
fac2.M4 = fac3(4)
fac2.M5 = fac3(5)
fac2.M6 = fac3(6)
'See if the signs of the two numbers
'are the same. If so, add; if not, subtract.
If fac1.sign = fac2.sign Then 'add
fpadd_aux fac1, fac2, m
Else
fpsub_aux fac1, fac2, m
End If
End If
result = fac1
'now check exponent of result, watch for overflow.
If result.exponent > 32768 Then
'er = EXPO_ERR
Print "add: Exponent over/under flow"
End If
End Sub
Sub fpsub (result As decfloat, x As decfloat, y As decfloat)
Dim As decfloat fac1, fac2
fac1 = x
fac2 = y
fac2.sign = fac2.sign Xor &H8000
fpadd fac1, fac1, fac2
result = fac1
'now check exponent of result, watch for overflow.
If result.exponent > 32768 Then
'er = EXPO_ERR
Print "sub: Exponent over/under flow"
End If
End Sub
Sub fpmul (result As decfloat, x As decfloat, y As decfloat)
'Dim As decfloat fac1,fac2
Dim As Integer i, j, ex, den, num
Dim As _Integer64 digit, carry, prod
Dim As _Unsigned Long fac3(0 To DWORDS * 2 + 1)
Dim As Long fac1(0 To DWORDS - 1), fac2(0 To DWORDS - 1)
Dim As _Unsigned Integer fac1exponent, fac2exponent
Dim As Integer fac1sign, fac2sign
' fac1=x
' fac2=y
'check exponents. if either is zero,
'the result is zero
If x.exponent = 0 Or y.exponent = 0 Then 'result is zero...clear fac1.
result.exponent = 0
result.sign = 0
result.M0 = 0
result.M1 = 0
result.M2 = 0
result.M3 = 0
result.M4 = 0
result.M5 = 0
result.M6 = 0
Exit Sub
Else
fac1(0) = x.M0
fac1(1) = x.M1
fac1(2) = x.M2
fac1(3) = x.M3
fac1(4) = x.M4
fac1(5) = x.M5
fac1(6) = x.M6
fac1exponent = x.exponent
fac1sign = x.sign
fac2(0) = y.M0
fac2(1) = y.M1
fac2(2) = y.M2
fac2(3) = y.M3
fac2(4) = y.M4
fac2(5) = y.M5
fac2(6) = y.M6
fac2exponent = y.exponent
fac2sign = y.sign
'clear fac3 mantissa
For i = 0 To DWORDS * 2 + 1
fac3(i) = 0
Next
den = DWORDS - 1
While fac2(den) = 0
den = den - 1
Wend
num = DWORDS - 1
While fac1(num) = 0
num = num - 1
Wend
If num < den Then
'Swap fac1, fac2
fac2(0) = x.M0
fac2(1) = x.M1
fac2(2) = x.M2
fac2(3) = x.M3
fac2(4) = x.M4
fac2(5) = x.M5
fac2(6) = x.M6
fac2exponent = x.exponent
fac2sign = x.sign
fac1(0) = y.M0
fac1(1) = y.M1
fac1(2) = y.M2
fac1(3) = y.M3
fac1(4) = y.M4
fac1(5) = y.M5
fac1(6) = y.M6
fac1exponent = y.exponent
fac1sign = y.sign
Swap den, num
End If
For j = den To 0 Step -1
carry = 0
digit = fac2(j)
For i = num To 0 Step -1
prod = fac3(i + j + 1) + digit * fac1(i) + carry
carry = prod \ 1000000000
fac3(i + j + 1) = (prod Mod 1000000000)
Next
fac3(j) = carry
Next
j = 0
If fac3(0) < 10~& Then
j = 8
ElseIf fac3(0) < 100~& Then
j = 7
ElseIf fac3(0) < 1000~& Then
j = 6
ElseIf fac3(0) < 10000~& Then
j = 5
ElseIf fac3(0) < 100000~& Then
j = 4
ElseIf fac3(0) < 1000000~& Then
j = 3
ElseIf fac3(0) < 10000000~& Then
j = 2
ElseIf fac3(0) < 100000000~& Then
j = 1
End If
If j > 0 Then
LSHIFT_a fac3(), 7, j
End If
result.M0 = fac3(0)
result.M1 = fac3(1)
result.M2 = fac3(2)
result.M3 = fac3(3)
result.M4 = fac3(4)
result.M5 = fac3(5)
result.M6 = fac3(6)
End If
'now determine exponent of result.
'as you do...watch for overflow.
ex = x.exponent - BIAS + y.exponent
result.exponent = ex - j
If multrndup Then
If fac3(DWORDS) >= 500000000~& Then
Dim As decfloat fac
fac.exponent = ex - j
fac.sign = 0
fac.M0 = 0
fac.M1 = 0
fac.M2 = 0
fac.M3 = 0
fac.M4 = 0
fac.M5 = 0
fac.M6 = 1
fpadd_aux result, fac, 0
End If
End If
'determine the sign of the product
result.sign = x.sign Xor y.sign
NORM_FAC1 result
'now check exponent of result, watch for overflow.
If result.exponent > 32768 Then
'er = EXPO_ERR
Print "mul: Exponent over/under flow"
End If
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_out As decfloat, x As decfloat, y As decfloat)
Dim As Long fac1(6), fac2(6), r
Dim As Long i, er, is_power_of_ten, rndup
Dim As _Unsigned Integer fac1exponent, fac2exponent
Dim As Integer fac1sign, fac2sign
fac1(0) = x.M0
fac1(1) = x.M1
fac1(2) = x.M2
fac1(3) = x.M3
fac1(4) = x.M4
fac1(5) = x.M5
fac1(6) = x.M6
fac1exponent = x.exponent
fac1sign = x.sign
fac2(0) = y.M0
fac2(1) = y.M1
fac2(2) = y.M2
fac2(3) = y.M3
fac2(4) = y.M4
fac2(5) = y.M5
fac2(6) = y.M6
fac2exponent = y.exponent
fac2sign = y.sign
If fac2exponent = 0 Then ' if fac2 = 0, return
' a divide-by-zero error and
' bail out.
result_out.M0 = 999999999
result_out.M1 = 999999999
result_out.M2 = 999999999
result_out.M3 = 999999999
result_out.M4 = 999999999
result_out.M5 = 999999999
result_out.M6 = 999999999
result_out.exponent = 9999 + BIAS + 1
'er = DIVZ_ERR
Print "division by zero"
Exit Sub
ElseIf fac1exponent = 0 Then 'fact1=0, just return
er = 0
result_out.M0 = 0
result_out.M1 = 0
result_out.M2 = 0
result_out.M3 = 0
result_out.M4 = 0
result_out.M5 = 0
result_out.M6 = 0
result_out.exponent = 0
result_out.sign = 0
Exit Sub
Else
'check to see if fac2 is a power of ten
is_power_of_ten = 0
If fac2(0) = 100000000 Then
is_power_of_ten = 1
For i = 1 To 6
If fac2(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
result_out.sign = fac1sign Xor fac2sign
result_out.exponent = fac1exponent - fac2exponent + BIAS + 1
result_out.M0 = fac1(0)
result_out.M1 = fac1(1)
result_out.M2 = fac1(2)
result_out.M3 = fac1(3)
result_out.M4 = fac1(4)
result_out.M5 = fac1(5)
result_out.M6 = fac1(6)
Exit Sub
End If
Const dm = 18
Dim As Double result(1 To dm), n(1 To dm), d(1 To dm)
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 6
' n(2 * j + 2) = fac1(j) \ 10000
' n(2 * j + 3) = fac1(j) Mod 10000
' d(2 * j + 2) = fac2(j) \ 10000
' d(2 * j + 3) = fac2(j) Mod 10000
' Next
n(2) = x.M0 \ 100000: r = x.M0 Mod 100000
n(3) = r \ 10: r = r Mod 10
n(4) = r * 1000 + x.M1 \ 1000000: r = x.M1 Mod 1000000
n(5) = r \ 100: r = r Mod 100
n(6) = r * 100 + x.M2 \ 10000000: r = x.M2 Mod 10000000
n(7) = r \ 1000: r = r Mod 1000
n(8) = r * 10 + x.M3 \ 100000000: r = x.M3 Mod 100000000
n(9) = r \ 10000: r = r Mod 10000
n(10) = r
n(11) = x.M4 \ 100000: r = x.M4 Mod 100000
n(12) = r \ 10: r = r Mod 10
n(13) = r * 1000 + x.M5 \ 1000000: r = x.M5 Mod 1000000
n(14) = r \ 100: r = r Mod 100
n(15) = r * 100 + x.M6 \ 10000000: r = x.M6 Mod 10000000
n(16) = r \ 1000: r = r Mod 1000
n(17) = r * 10
d(2) = y.M0 \ 100000: r = y.M0 Mod 100000
d(3) = r \ 10: r = r Mod 10
d(4) = r * 1000 + y.M1 \ 1000000: r = y.M1 Mod 1000000
d(5) = r \ 100: r = r Mod 100
d(6) = r * 100 + y.M2 \ 10000000: r = y.M2 Mod 10000000
d(7) = r \ 1000: r = r Mod 1000
d(8) = r * 10 + y.M3 \ 100000000: r = y.M3 Mod 100000000
d(9) = r \ 10000: r = r Mod 10000
d(10) = r
d(11) = y.M4 \ 100000: r = y.M4 Mod 100000
d(12) = r \ 10: r = r Mod 10
d(13) = r * 1000 + y.M5 \ 1000000: r = y.M5 Mod 1000000
d(14) = r \ 100: r = r Mod 100
d(15) = r * 100 + y.M6 \ 10000000: r = y.M6 Mod 10000000
d(16) = r \ 1000: r = r Mod 1000
d(17) = r * 10
n(1) = (fac1exponent And &H7FFF) - BIAS - 1
d(1) = (fac2exponent And &H7FFF) - 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))
subtract w(), q, d(), (stp + 2), last
normalize w(), (stp + 2), q
Next
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)
j = 0
If result(2) < 10~& Then
j = 3
ElseIf result(2) < 100~& Then
j = 2
ElseIf result(2) < 1000~& Then
j = 1
End If
If j > 0 Then
LSHIFT_da result(), dm, j
End If
j = 2
result_out.M0 = (result(j) * 10000 + result(j + 1)) * 10 + result(j + 2) \ 1000: r = result(j + 2) Mod 1000: j = 5
result_out.M1 = (r * 10000 + result(j)) * 100 + result(j + 1) \ 100: r = result(j + 1) Mod 100: j = 7
result_out.M2 = (r * 10000 + result(j)) * 1000 + result(j + 1) \ 10: r = result(j + 1) Mod 10: j = 9
result_out.M3 = (r * 10000 + result(j)) * 10000 + result(j + 1): j = 11
result_out.M4 = (result(j) * 10000 + result(j + 1)) * 10 + result(j + 2) \ 1000: r = result(j + 2) Mod 1000: j = 14
result_out.M5 = (r * 10000 + result(j)) * 100 + result(j + 1) \ 100: r = result(j + 1) Mod 100: j = 16
result_out.M6 = (r * 10000 + result(j)) * 1000 + result(j + 1) \ 10
rndup = 0
If (result(j + 1) Mod 10) > 4 Then
Dim As decfloat ru
ru.M0 = 100000000
ru.exponent = 1 + BIAS
result_out.exponent = 63 + BIAS
fpadd result_out, result_out, ru
If result_out.exponent > (63 + BIAS) Then
rndup = 1
End If
End If
NORM_FAC1 result_out
result_out.exponent = (result(1) + BIAS) + rndup
End If
result_out.sign = fac1sign Xor fac2sign
End Sub
Sub fpdiv_si (result As decfloat, num As decfloat, den As Long)
Dim As Long fac1(6)
Dim As _Unsigned Integer fac1exponent
Dim As Integer fac1sign
Dim As _Unsigned _Integer64 carry, remder
Dim As _Integer64 i, divisor
Dim As _Integer64 quotient
remder = 0
divisor = Abs(den)
fac1(0) = num.M0
fac1(1) = num.M1
fac1(2) = num.M2
fac1(3) = num.M3
fac1(4) = num.M4
fac1(5) = num.M5
fac1(6) = num.M6
fac1exponent = num.exponent
fac1sign = num.sign
result.M0 = num.M0
result.M1 = num.M1
result.M2 = num.M2
result.M3 = num.M3
result.M4 = num.M4
result.M5 = num.M5
result.M6 = num.M6
result.exponent = num.exponent
result.sign = num.sign
If divisor = 0 Then
Print "error: divisor = 0"
Exit Sub
End If
If divisor > 2147483647 Then
Print "error: divisor too large"
Exit Sub
End If
For i = 0 To 6
quotient = fac1(i) + remder * 1000000000
remder = quotient Mod divisor
fac1(i) = quotient \ divisor
Next
quotient = remder * 1000000000
quotient = quotient \ divisor
carry = fac1(0)
result.M0 = fac1(0)
result.M1 = fac1(1)
result.M2 = fac1(2)
result.M3 = fac1(3)
result.M4 = fac1(4)
result.M5 = fac1(5)
result.M6 = fac1(6)
result.exponent = fac1exponent
result.sign = fac1sign
If carry = 0 Then
LSHIFT_n result, 9
result.exponent = result.exponent - 9
result.M6 = result.M6 + quotient
ElseIf carry < 10 Then
LSHIFT_n result, 8
result.exponent = result.exponent - 8
result.M6 = result.M6 + quotient \ 10
ElseIf carry < 100 Then
LSHIFT_n result, 7
result.exponent = result.exponent - 7
result.M6 = result.M6 + quotient \ 100
ElseIf carry < 1000 Then
LSHIFT_n result, 6
result.exponent = result.exponent - 6
result.M6 = result.M6 + quotient \ 1000
ElseIf carry < 10000 Then
LSHIFT_n result, 5
result.exponent = result.exponent - 5
result.M6 = result.M6 + quotient \ 10000
ElseIf carry < 100000 Then
LSHIFT_n result, 4
result.exponent = result.exponent - 4
result.M6 = result.M6 + quotient \ 100000
ElseIf carry < 1000000 Then
LSHIFT_n result, 3
result.exponent = result.exponent - 3
result.M6 = result.M6 + quotient \ 1000000
ElseIf carry < 10000000 Then
LSHIFT_n result, 2
result.exponent = result.exponent - 2
result.M6 = result.M6 + quotient \ 10000000
ElseIf carry < 100000000 Then
LSHIFT_n result, 1
result.exponent = result.exponent - 1
result.M6 = result.M6 + quotient \ 100000000
End If
'NORM_FAC1(fac1)
result.sign = fac1sign
If den < 0 Then
result.sign = fac1sign Xor &H8000
End If
'now check exponent of result, watch for overflow.
If result.exponent > 32768 Then
'er = EXPO_ERR
Print "div_si: Exponent over/under flow"
End If
End Sub
Sub fpmul_si (result As decfloat, x As decfloat, y As _Integer64)
Dim As decfloat fac1, fac2
Dim As Long i
Dim As _Integer64 carry, digit, prod, value
Dim As _Unsigned Long fac3(7)
fac1 = x
digit = Abs(y)
If digit > 999999999 Then
Si2fp fac2, y
fpmul fac1, fac1, fac2
result = fac1: Exit Sub
End If
If fac1.exponent = 0 Or y = 0 Then 'result is zero...clear fac1.
fac1.sign = 0
fac1.exponent = 0
fac1.M0 = 0
fac1.M1 = 0
fac1.M2 = 0
fac1.M3 = 0
fac1.M4 = 0
fac1.M5 = 0
fac1.M6 = 0
NORM_FAC1 fac1
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
fac3(0) = fac1.M0
fac3(1) = fac1.M1
fac3(2) = fac1.M2
fac3(3) = fac1.M3
fac3(4) = fac1.M4
fac3(5) = fac1.M5
fac3(6) = fac1.M6
fac3(7) = 0
carry = 0
For i = 6 To 0 Step -1
prod = digit * fac3(i) + carry
value = prod Mod 1000000000
fac3(i) = value
carry = prod \ 1000000000
Next
If carry < 10 Then
RSHIFT_a fac3(), 7, 1
fac1.exponent = fac1.exponent + 1
fac3(0) = fac3(0) + carry * 100000000
ElseIf carry < 100 Then
RSHIFT_a fac3(), 7, 2
fac1.exponent = fac1.exponent + 2
fac3(0) = fac3(0) + carry * 10000000
ElseIf carry < 1000 Then
RSHIFT_a fac3(), 7, 3
fac1.exponent = fac1.exponent + 3
fac3(0) = fac3(0) + carry * 1000000
ElseIf carry < 10000 Then
RSHIFT_a fac3(), 7, 4
fac1.exponent = fac1.exponent + 4
fac3(0) = fac3(0) + carry * 100000
ElseIf carry < 100000 Then
RSHIFT_a fac3(), 7, 5
fac1.exponent = fac1.exponent + 5
fac3(0) = fac3(0) + carry * 10000
ElseIf carry < 1000000 Then
RSHIFT_a fac3(), 7, 6
fac1.exponent = fac1.exponent + 6
fac3(0) = fac3(0) + carry * 1000
ElseIf carry < 10000000 Then
RSHIFT_a fac3(), 7, 7
fac1.exponent = fac1.exponent + 7
fac3(0) = fac3(0) + carry * 100
ElseIf carry < 100000000 Then
RSHIFT_a fac3(), 7, 8
fac1.exponent = fac1.exponent + 8
fac3(0) = fac3(0) + carry * 10
ElseIf carry < 1000000000 Then
RSHIFT_a fac3(), 7, 9
fac1.exponent = fac1.exponent + 9
fac3(0) = fac3(0) + carry
End If
End If
fac1.M0 = fac3(0)
fac1.M1 = fac3(1)
fac1.M2 = fac3(2)
fac1.M3 = fac3(3)
fac1.M4 = fac3(4)
fac1.M5 = fac3(5)
fac1.M6 = fac3(6)
If fac3(7) >= 500000000~& Then
Dim As decfloat fac
fac.exponent = fac1.exponent
fac.sign = 0
fac.M0 = 0
fac.M1 = 0
fac.M2 = 0
fac.M3 = 0
fac.M4 = 0
fac.M5 = 0
fac.M6 = 1
fpadd_aux fac1, fac, 0
End If
NORM_FAC1 fac1
If y < 0 Then
fac1.sign = fac1.sign Xor &H8000
End If
result = fac1
'now check exponent of result, watch for overflow.
If result.exponent > 32768 Then
'er = EXPO_ERR
Print "mul_si: Exponent over/under flow"
End If
End Sub
'integer part of num
Sub fpfix (result As decfloat, num As decfloat)
Dim As decfloat ips
Dim As Long ex, ex2, j, k
ex = (num.exponent And &H7FFF) - BIAS
If ex < 1 Then
result = ips
Exit Sub
End If
If ex >= 63 Then
result = num
Exit Sub
End If
Dim As _Unsigned Long ip(6), numa(6)
numa(0) = num.M0
numa(1) = num.M1
numa(2) = num.M2
numa(3) = num.M3
numa(4) = num.M4
numa(5) = num.M5
numa(6) = num.M6
ex2 = ex \ 9
k = ex2
j = ex Mod 9
While ex2 > 0
ex2 = ex2 - 1
ip(ex2) = numa(ex2)
Wend
If j = 1 Then
ip(k) = 100000000 * (numa(k) \ 100000000)
ElseIf j = 2 Then
ip(k) = 10000000 * (numa(k) \ 10000000)
ElseIf j = 3 Then
ip(k) = 1000000 * (numa(k) \ 1000000)
ElseIf j = 4 Then
ip(k) = 100000 * (numa(k) \ 100000)
ElseIf j = 5 Then
ip(k) = 10000 * (numa(k) \ 10000)
ElseIf j = 6 Then
ip(k) = 1000 * (numa(k) \ 1000)
ElseIf j = 7 Then
ip(k) = 100 * (numa(k) \ 100)
ElseIf j = 8 Then
ip(k) = 10 * (numa(k) \ 10)
ElseIf j = 9 Then
ip(k) = numa(k)
End If
result.exponent = ex + BIAS
result.sign = num.sign
result.M0 = ip(0)
result.M1 = ip(1)
result.M2 = ip(2)
result.M3 = ip(3)
result.M4 = ip(4)
result.M5 = ip(5)
result.M6 = ip(6)
End Sub
Sub fpfrac (result As decfloat, num As decfloat)
Dim As decfloat ip, fp
fpfix ip, num
fpsub fp, num, ip
result = fp
End Sub
Function fpfix_is_odd& (num As decfloat)
Dim As Long ex, ex2, j, k
ex = (num.exponent And &H7FFF) - BIAS
If ex < 0 Then
fpfix_is_odd = 0
Exit Function
End If
If ex >= 63 Then
fpfix_is_odd = 0
Exit Function
End If
Dim As _Unsigned Long numa(6)
numa(0) = num.M0
numa(1) = num.M1
numa(2) = num.M2
numa(3) = num.M3
numa(4) = num.M4
numa(5) = num.M5
numa(6) = num.M6
ex = ex - 1
ex2 = ex \ 9
k = ex2
j = ex Mod 9
If j = 0 Then
fpfix_is_odd = (numa(k) \ 100000000) And 1: Exit Function
ElseIf j = 1 Then
fpfix_is_odd = (numa(k) \ 10000000) And 1: Exit Function
ElseIf j = 2 Then
fpfix_is_odd = (numa(k) \ 1000000) And 1: Exit Function
ElseIf j = 3 Then
fpfix_is_odd = (numa(k) \ 100000) And 1: Exit Function
ElseIf j = 4 Then
fpfix_is_odd = (numa(k) \ 10000) And 1: Exit Function
ElseIf j = 5 Then
fpfix_is_odd = (numa(k) \ 1000) And 1: Exit Function
ElseIf j = 6 Then
fpfix_is_odd = (numa(k) \ 100) And 1: Exit Function
ElseIf j = 7 Then
fpfix_is_odd = (numa(k) \ 10) And 1: Exit Function
ElseIf j = 8 Then
fpfix_is_odd = numa(k) And 1: Exit Function
End If
fpfix_is_odd = 0
End Function
Posts: 2,186
Threads: 222
Joined: Apr 2022
Reputation:
104
That's the first program I tried on the new QB64PE v3.3.
Nice!
Question. I can do SQR() to get square roots. Any notation for nth roots like cube roots 4th roots, etc.?
It looks like it does decimal roots. If you don't have a method set up for nth roots yet, all you would need is your own notation like 3r() and a function to manipulate the equation as a decimal root. You know, the cube root of 27 is also 27 ^ .333..., etc. (Provided the rounding formula handles it correctly, which it looks like would require running the repeating portion out a certain number of decimals past the reporting decimal limit.)
Pete
|
