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