Posts: 733
Threads: 30
Joined: Apr 2022
Reputation:
43
(09-15-2022, 03:37 AM)Pete Wrote: Beats "Pride" month. Rainbows are just half-ascii circles.
Pete
Technically, a rainbow is a circle. We just don't usually get to see the other half because of the horizon or other obstructions.
Tread on those who tread on you
Posts: 422
Threads: 27
Joined: Apr 2022
Reputation:
26
09-16-2022, 09:03 PM
(This post was last modified: 09-16-2022, 09:04 PM by Jack.)
just for the fun of it here's a version with 56-digit and maximum exponent of 16384
Code: (Select All) $NoPrefix
$Console:Only
Dest Console
Dim Shared As Long shift_constants(6)
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
Const BIAS = 16384 '2 ^ 14
Type decfloat
As Unsigned Integer exponent
As Integer sign
As Long M0
As Long M1
As Long M2
As Long M3
As Long M4
As Long M5
As Long M6
End Type
' Error definitions
Const DIVZ_ERR = 1 'Divide by zero
Const EXPO_ERR = 2 'Exponent overflow error
Const EXPU_ERR = 3 'Exponent underflow error
Dim As decfloat x, y, z, pi1
Dim As Integer64 i, j, k
'si2fp x, -9223372036854775808
'str2fp y, "2.7182818284590452353602874713526624977572470936999595750"
str2fp pi1, "3.1415926535897932384626433832795028841971693993751058210"
str2fp x, "9.9999999999999999999999999999999999999999999999999999999"
str2fp y, "8.8888888888888888888888888888888888888888888888888888888"
Print "pi1 = "; fp2str(pi1)
Print "x = "; fp2str(x)
Print "y = "; fp2str(y)
fpadd z, x, y
Print "x + y = "; fp2str(z)
fpsub z, x, y
Print "x - y = "; fp2str(z)
fpmul z, x, y
Print "x * y = "; fp2str(z)
fpdiv z, x, y
Print "x / y = "; fp2str(z)
fpdiv_si z, pi1, 2
Print "Pi / 2 = "; fp2str(z)
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) > (56 + 1 + 8) Then
f1 = Mid$(f1, 1, (56 + 1 + 8))
End If
While Len(f1) < (56 + 1 + 8)
f1 = f1 + "0"
Wend
j = 1
ts = Mid$(f1, j, 8)
ulng = Val(ts)
n.M0 = ulng
If ulng <> 0 Then fp = 1
j = j + 8
ts = Mid$(f1, j, 8)
ulng = Val(ts)
n.M1 = ulng
If ulng <> 0 Then fp = 1
j = j + 8
ts = Mid$(f1, j, 8)
ulng = Val(ts)
n.M2 = ulng
If ulng <> 0 Then fp = 1
j = j + 8
ts = Mid$(f1, j, 8)
ulng = Val(ts)
n.M3 = ulng
If ulng <> 0 Then fp = 1
j = j + 8
ts = Mid$(f1, j, 8)
ulng = Val(ts)
n.M4 = ulng
If ulng <> 0 Then fp = 1
j = j + 8
ts = Mid$(f1, j, 8)
ulng = Val(ts)
n.M5 = ulng
If ulng <> 0 Then fp = 1
j = j + 8
ts = Mid$(f1, j, 8)
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
Function fp2str$ (n As decfloat)
Dim As Long ex, ex1, i, ln
Dim As String s, sd, sdl, sdr, se, sign
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 < 8 Then
s = String$(8 - ln, "0") + s
End If
sd = sd + s
s = Trim$(Str$(n.M2))
ln = Len(s)
If ln < 8 Then
s = String$(8 - ln, "0") + s
End If
sd = sd + s
s = Trim$(Str$(n.M3))
ln = Len(s)
If ln < 8 Then
s = String$(8 - ln, "0") + s
End If
sd = sd + s
s = Trim$(Str$(n.M4))
ln = Len(s)
If ln < 8 Then
s = String$(8 - ln, "0") + s
End If
sd = sd + s
s = Trim$(Str$(n.M5))
ln = Len(s)
If ln < 8 Then
s = String$(8 - ln, "0") + s
End If
sd = sd + s
s = Trim$(Str$(n.M6))
ln = Len(s)
If ln < 8 Then
s = String$(8 - ln, "0") + s
End If
sd = sd + s
ln = Len(sd)
If ex >= 0 Then
If ex < 55 Then
sd = Left$(sd, ex + 1) + "." + Mid$(sd, ex + 2)
ElseIf ex > 55 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
Sub si2fp (result As decfloat, m As Integer64)
Dim As decfloat fac1
Dim As Long i
Dim As _Unsigned Integer64 n
n = Abs(m)
If n > 9999999999999999 Then
Call str2fp(fac1, Str$(m))
result = fac1: 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
result = fac1: Exit Sub
End If
fac1.exponent = BIAS
If n < 100000000 Then
If n < 10 Then
fac1.M0 = n * 10000000
fac1.exponent = fac1.exponent + 1
ElseIf n < 100 Then
fac1.M0 = n * 1000000
fac1.exponent = fac1.exponent + 2
ElseIf n < 1000 Then
fac1.M0 = n * 100000
fac1.exponent = fac1.exponent + 3
ElseIf n < 10000 Then
fac1.M0 = n * 10000
fac1.exponent = fac1.exponent + 4
ElseIf n < 100000 Then
fac1.M0 = n * 1000
fac1.exponent = fac1.exponent + 5
ElseIf n < 1000000 Then
fac1.M0 = n * 100
fac1.exponent = fac1.exponent + 6
ElseIf n < 10000000 Then
fac1.M0 = n * 10
fac1.exponent = fac1.exponent + 7
ElseIf n < 100000000 Then
fac1.M0 = n
fac1.exponent = fac1.exponent + 8
End If
End If
If n > 99999999 Then
fac1.exponent = fac1.exponent + 8
If n < 1000000000 Then
fac1.M0 = n \ 10
fac1.M1 = (n Mod 10) * 10000000
fac1.exponent = fac1.exponent + 1
ElseIf n < 100000000000 Then
fac1.M0 = n \ 100
fac1.M1 = (n Mod 100) * 1000000
fac1.exponent = fac1.exponent + 2
ElseIf n < 1000000000000 Then
fac1.M0 = n \ 1000
fac1.M1 = (n Mod 1000) * 100000
fac1.exponent = fac1.exponent + 3
ElseIf n < 10000000000000 Then
fac1.M0 = n \ 10000
fac1.M1 = (n Mod 10000) * 10000
fac1.exponent = fac1.exponent + 4
ElseIf n < 100000000000000 Then
fac1.M0 = n \ 100000
fac1.M1 = (n Mod 100000) * 1000
fac1.exponent = fac1.exponent + 5
ElseIf n < 1000000000000000 Then
fac1.M0 = n \ 1000000
fac1.M1 = (n Mod 1000000) * 100
fac1.exponent = fac1.exponent + 6
ElseIf n < 10000000000000000 Then
fac1.M0 = n \ 10000000
fac1.M1 = (n Mod 10000000) * 10
fac1.exponent = fac1.exponent + 7
ElseIf n < 100000000000000000 Then
fac1.M0 = n \ 100000000
fac1.M1 = n Mod 100000000
fac1.exponent = fac1.exponent + 8
End If
End If
If m < 0 Then
fac1.M0 = fac1.M0 Or &H80000000
End If
result = fac1
End Sub
Sub RSHIFT_n (mantissa As decfloat, n As Long)
If n = 8 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(7 - 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 = 8 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(7 - 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
Function fpcmp& (x As decfloat, y As decfloat)
Dim As Long i
Dim As Integer64 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 i, 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 > 99999999 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, 8
fac1.exponent = fac1.exponent - 8
If fac1.exponent = 0 Then
Print "NORM_FAC1=EXPU_ERR"
Exit Sub
End If
Wend
If fac1.M0 < 10 Then
LSHIFT_n fac1, 7
fac1.exponent = fac1.exponent - 7
ElseIf fac1.M0 < 100 Then
LSHIFT_n fac1, 6
fac1.exponent = fac1.exponent - 6
ElseIf fac1.M0 < 1000 Then
LSHIFT_n fac1, 5
fac1.exponent = fac1.exponent - 5
ElseIf fac1.M0 < 10000 Then
LSHIFT_n fac1, 4
fac1.exponent = fac1.exponent - 4
ElseIf fac1.M0 < 100000 Then
LSHIFT_n fac1, 3
fac1.exponent = fac1.exponent - 3
ElseIf fac1.M0 < 1000000 Then
LSHIFT_n fac1, 2
fac1.exponent = fac1.exponent - 2
ElseIf fac1.M0 < 10000000 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, i
c = 0
v = fac2.M6 + fac1.M6 + c
If v > 99999999 Then
v = v - 100000000
c = 1
Else
c = 0
End If
fac1.M6 = v
v = fac2.M5 + fac1.M5 + c
If v > 99999999 Then
v = v - 100000000
c = 1
Else
c = 0
End If
fac1.M5 = v
v = fac2.M4 + fac1.M4 + c
If v > 99999999 Then
v = v - 100000000
c = 1
Else
c = 0
End If
fac1.M4 = v
v = fac2.M3 + fac1.M3 + c
If v > 99999999 Then
v = v - 100000000
c = 1
Else
c = 0
End If
fac1.M3 = v
v = fac2.M2 + fac1.M2 + c
If v > 99999999 Then
v = v - 100000000
c = 1
Else
c = 0
End If
fac1.M2 = v
v = fac2.M1 + fac1.M1 + c
If v > 99999999 Then
v = v - 100000000
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, i
c = 0
v = fac1.M6 - fac2.M6 - c
If v < 0 Then
v = v + 100000000
c = 1
Else
c = 0
End If
fac1.M6 = v
v = fac1.M5 - fac2.M5 - c
If v < 0 Then
v = v + 100000000
c = 1
Else
c = 0
End If
fac1.M5 = v
v = fac1.M4 - fac2.M4 - c
If v < 0 Then
v = v + 100000000
c = 1
Else
c = 0
End If
fac1.M4 = v
v = fac1.M3 - fac2.M3 - c
If v < 0 Then
v = v + 100000000
c = 1
Else
c = 0
End If
fac1.M3 = v
v = fac1.M2 - fac2.M2 - c
If v < 0 Then
v = v + 100000000
c = 1
Else
c = 0
End If
fac1.M2 = v
v = fac1.M1 - fac2.M1 - c
If v < 0 Then
v = v + 100000000
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 < 56 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 < 56 Then 'shift
i = t \ 8
While i > 0
RSHIFT_n fac2, 8
t = t - 8
i = i - 1
Wend
If 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
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
End Sub
Sub fpmul (result As decfloat, x As decfloat, y As decfloat)
'Dim As decfloat fac1,fac2
Dim As Integer i, j, ex, er, den, num
Dim As Integer64 digit, carry, prod
Dim As Unsigned Long fac3(0 To 7 * 2 + 1)
Dim As Long fac1(0 To 6), fac2(0 To 6)
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 7 * 2 + 1
fac3(i) = 0
Next
den = 6
While fac2(den) = 0
den = den - 1
Wend
num = 6
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 \ 100000000
fac3(i + j + 1) = (prod Mod 100000000)
Next
fac3(j) = carry
Next
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
'determine the sign of the product
result.sign = x.sign Xor y.sign
NORM_FAC1 result
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)
Dim As Long i, er, is_power_of_ten
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 = 99999999
result_out.M1 = 99999999
result_out.M2 = 99999999
result_out.M3 = 99999999
result_out.M4 = 99999999
result_out.M5 = 99999999
result_out.M6 = 99999999
result_out.exponent = 9999 + BIAS + 1
er = DIVZ_ERR
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) = 10000000 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
Dim As Double result(1 To 15), n(1 To 15), d(1 To 15)
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 15 + 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(1) = (fac1exponent And &H7FFF) - BIAS - 1
d(1) = (fac2exponent And &H7FFF) - BIAS - 1
For j = 15 To 19
w(j) = 0
Next
t = 14
w(1) = n(1) - d(1) + 1
w(2) = 0
For j = 2 To 15
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, 19)
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
result_out.M0 = result(2 * j + 2) * 10000 + result(2 * j + 3): j = j + 1
result_out.M1 = result(2 * j + 2) * 10000 + result(2 * j + 3): j = j + 1
result_out.M2 = result(2 * j + 2) * 10000 + result(2 * j + 3): j = j + 1
result_out.M3 = result(2 * j + 2) * 10000 + result(2 * j + 3): j = j + 1
result_out.M4 = result(2 * j + 2) * 10000 + result(2 * j + 3): j = j + 1
result_out.M5 = result(2 * j + 2) * 10000 + result(2 * j + 3): j = j + 1
result_out.M6 = result(2 * j + 2) * 10000 + result(2 * j + 3)
result_out.exponent = fac1exponent
result_out.sign = fac1sign
NORM_FAC1 result_out
fac1exponent = (result(1) + BIAS)
End If
result_out.sign = fac1sign Xor fac2sign
result_out.exponent = fac1exponent
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 * 100000000
remder = quotient Mod divisor
fac1(i) = quotient \ divisor
Next
quotient = remder * 100000000
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, 8
result.exponent = result.exponent - 8
result.M6 = result.M6 + quotient
ElseIf carry < 10 Then
LSHIFT_n result, 7
result.exponent = result.exponent - 7
result.M6 = result.M6 + quotient \ 10
ElseIf carry < 100 Then
LSHIFT_n result, 6
result.exponent = result.exponent - 6
result.M6 = result.M6 + quotient \ 100
ElseIf carry < 1000 Then
LSHIFT_n result, 5
result.exponent = result.exponent - 5
result.M6 = result.M6 + quotient \ 1000
ElseIf carry < 10000 Then
LSHIFT_n result, 4
result.exponent = result.exponent - 4
result.M6 = result.M6 + quotient \ 10000
ElseIf carry < 100000 Then
LSHIFT_n result, 3
result.exponent = result.exponent - 3
result.M6 = result.M6 + quotient \ 100000
ElseIf carry < 1000000 Then
LSHIFT_n result, 2
result.exponent = result.exponent - 2
result.M6 = result.M6 + quotient \ 1000000
ElseIf carry < 10000000 Then
LSHIFT_n result, 1
result.exponent = result.exponent - 1
result.M6 = result.M6 + quotient \ 10000000
End If
'NORM_FAC1(fac1)
result.sign = fac1sign
If den < 0 Then
result.sign = fac1sign Xor &H8000
End If
End Sub
sample run
Code: (Select All) pi1 = 3.1415926535897932384626433832795028841971693993751058210
x = 9.9999999999999999999999999999999999999999999999999999999
y = 8.8888888888888888888888888888888888888888888888888888888
x + y = 18.888888888888888888888888888888888888888888888888888888
x - y = 1.1111111111111111111111111111111111111111111111111111111
x * y = 88.888888888888888888888888888888888888888888888888888887
x / y = 1.1250000000000000000000000000000000000000000000000000000
Pi / 2 = 1.5707963267948966192313216916397514420985846996875529105
Posts: 2,171
Threads: 222
Joined: Apr 2022
Reputation:
103
More involved, but this looks similar to something I was trying out before giving up on numeric computations.
Pete
Posts: 422
Threads: 27
Joined: Apr 2022
Reputation:
26
there are many thing that need to be added and tweaked, for example I can add 4 guard digits in the sign integer
Posts: 2,171
Threads: 222
Joined: Apr 2022
Reputation:
103
Yep. More like opening a drum of worms, vs a can.
Pete
Posts: 422
Threads: 27
Joined: Apr 2022
Reputation:
26
09-19-2022, 01:48 AM
(This post was last modified: 09-19-2022, 01:49 AM by Jack.)
Pete, why did you give up on numeric computations?
it took me a long time, but dogged determination brought me success, writing math routines from scratch can be very hard to debug and to get working because of many dependent routines.
just for fun I am continuing to develop this 32-byte decfloat package but I decided to change from 8 digits per element to 9, which took me a while to get working right, this way you get 63 digits but some operations tend to loose accuracy so the the print routine is restricted to 55 digits with round up, I got the basic 4 arithmetic and some helper functions done
Posts: 733
Threads: 30
Joined: Apr 2022
Reputation:
43
I want to say there is even a library for a similar number format. Can not remember what it is called, though.
Tread on those who tread on you
Posts: 2,171
Threads: 222
Joined: Apr 2022
Reputation:
103
09-20-2022, 02:30 AM
(This post was last modified: 09-20-2022, 02:31 AM by Pete.)
(09-19-2022, 01:48 AM)Jack Wrote: Pete, why did you give up on numeric computations?
it took me a long time, but dogged determination brought me success, writing math routines from scratch can be very hard to debug and to get working because of many dependent routines.
just for fun I am continuing to develop this 32-byte decfloat package but I decided to change from 8 digits per element to 9, which took me a while to get working right, this way you get 63 digits but some operations tend to loose accuracy so the the print routine is restricted to 55 digits with round up, I got the basic 4 arithmetic and some helper functions done
@Jack
Keep in mind I started migrating to string math before some of the better underscore math types were built for QB64. When I took a second look, and tried the larger data types, I got some stuff working with numeric computation but sure enough, when I ramped up the size of the the computations enough, fail.
With string math very large computations involving addition, subtraction, multiplication, non-decimal powers, square roots, and percentages are done remarkably well and now, fast.
What I'm finding a problem with is general roots (includes decimal roots) and decimal powers. With solutions involving approximation routines, especially those which rely on decimal factors in the equation like .333, .666, etc. well, errors just keep getting compounded. There is way too much irony involved with having a math system that can print out 10,000 digits in a decent amount of time but, in some instances, every one of those digits could be off!
So what is needed with any of the mat routines is fraction manipulations with long division when needed. That seems quite difficult to develop for general roots and decimal powers.
Pete
Posts: 422
Threads: 27
Joined: Apr 2022
Reputation:
26
rational arithmetic has it's place, but the fractions quite often become really huge
Posts: 2,171
Threads: 222
Joined: Apr 2022
Reputation:
103
(09-20-2022, 12:18 PM)Jack Wrote: rational arithmetic has it's place, but the fractions quite often become really huge
They do, indeed. Of course you can always waste a lot of time doing GCF iterations to reduce those fractions, but I found it is faster with string math to just let those huge fractions grow. That is especially true when you get complex fractions that hardly reduce at all, but a lot of cpu time is spent trying to see if they can be reduced. At the end of a routine, the numerator is simply divided by the denominator to display a decimal representation, so it really doesn't matter if even the end resulting fraction gets reduced. For ongoing calculations, that might slow things down to a point where a GFC reduction should be utilized.
Pete
Shoot first and shoot people who ask questions, later.
|
