decfloat -- again

[Jack]

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