Author Topic: decfloat 2.0 (Read 10363 times)

jack · « **on:** October 10, 2021, 08:00:41 am »

FreeBASIC supports arrays in UDT's and I had written this using arrays in the decfloat type
so in order to test how much slower this routines would perform using a string * some_number
I rewrote it using the same string scheme as you see here, for 1000 digit precision it was 178 times slower

Code: QB64: [Select]

$NoPrefix
$Console:Only
Dest Console
 
Option Explicit
 
Const NUMBER_OF_DIGITS = 512
Const NUM_DIGITS = NUMBER_OF_DIGITS
Const NUM_DWORDS = NUM_DIGITS \ 8
Const BIAS = 1073741824 '2 ^ 30
Const MANTISSA_BYTES = (NUM_DWORDS + 1) * 8 * 4
 
Type decfloat
    As Long sign
    As Unsigned Long exponent
    As String * Mantissa_bytes mantissa
End Type
 
' Error definitions
 
Const DIVZ_ERR = 1 'Divide by zero
Const EXPO_ERR = 2 'Exponent overflow error
Const EXPU_ERR = 3 'Exponent underflow error
Dim Shared As decfloat pi_dec
Call pi_brent_salamin(pi_dec, NUMBER_OF_DIGITS)
 
Dim As decfloat n, m, z
 
Dim As Double t
t = Timer
Call fpfactorial(z, 10000, NUMBER_OF_DIGITS)
t = Timer - t
Print "factorial(10000) to "; NUMBER_OF_DIGITS; " digits"
Print fp2str(z, NUMBER_OF_DIGITS)
Print
Print "in "; t; " seconds"
Print
Call str2fp(n, "3")
Call str2fp(m, "9")
t = Timer
Call fpdiv(z, n, m, NUMBER_OF_DIGITS)
t = Timer - t
Print "3/9 to "; NUMBER_OF_DIGITS; " digits"
Print fp2str(z, NUMBER_OF_DIGITS)
Print
Print "in "; t; " seconds"
Print
t = Timer
Call fpmul(z, z, z, NUMBER_OF_DIGITS)
t = Timer - t
Print "and square it"
Print fp2str(z, NUMBER_OF_DIGITS)
Print
Call str2fp(n, "2")
t = Timer
Call fplog(z, n, NUMBER_OF_DIGITS)
t = Timer - t
Print "ln(2) to "; NUMBER_OF_DIGITS; " digits"
Print fp2str(z, NUMBER_OF_DIGITS)
Print
Print "in "; t; " seconds"
Print
t = Timer
Call fpexp(z, z, NUMBER_OF_DIGITS)
t = Timer - t
Print "and the inverse"
Print fp2str(z, NUMBER_OF_DIGITS)
Print
Print "in "; t; " seconds"
Call si2fp(n, 3, NUMBER_OF_DIGITS)
Print
Print "reciprocal of 3"
t = Timer
Call fpinv(z, n, NUMBER_OF_DIGITS)
t = Timer - t
Print fp2str(z, NUMBER_OF_DIGITS)
Print
Print "in "; t; " seconds"
 
Sub set (s$, i&, v~&)
    Mid$(s$, 4 * i& + 1, 4) = MKL$(v~&)
End Sub
 
Function git~& (s$, i&)
    git = CVL(Mid$(s$, 4 * i& + 1, 4))
End Function
 
Function fp2str_exp$ (n As decfloat, places As Long)
    Dim As Long i, ex
    Dim As String v, f, ts
    If n.exponent > 0 Then
        ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
    Else
        ex = 0
    End If
    If n.sign Then v = "-" Else v = " "
    ts = Str$(git(n.mantissa, 0))
    ts = Trim$(ts)
    If Len(ts) < 8 Then
        ts = ts + String$(8 - Len(ts), "0")
    End If
    v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
    For i = 1 To NUM_DWORDS - 1
        ts = Str$(git(n.mantissa, i))
        ts = Trim$(ts)
        If Len(ts) < 8 Then
            ts = String$(8 - Len(ts), "0") + ts
        End If
        v = v + ts
    Next
    v = Left$(v, places + 3)
    f = Trim$(Str$(Abs(ex)))
    f = String$(5 - Len(f), "0") + f
    If ex < 0 Then v = v + "E-" Else v = v + "E+"
    v = v + f
    fp2str_exp$ = v
End Function
 
'long part of num
Sub fpfix (result As decfloat, num As decfloat)
    Dim As decfloat ip
    Dim As Long ex, ex2, j, k
 
    ex = (num.exponent And &H7FFFFFFF) - BIAS
    If ex < 1 Then
        result = ip: Exit Sub
    End If
    If ex >= (NUM_DIGITS) Then
        result = num: Exit Sub
    End If
    ex2 = ex \ 8
    k = ex2
    j = ex Mod 8
    While ex2 > 0
        ex2 = ex2 - 1
        Call set(ip.mantissa, ex2, git(num.mantissa, ex2))
    Wend
    If j = 1 Then
        Call set(ip.mantissa, k, 10000000 * (git(num.mantissa, k) \ 10000000))
    ElseIf j = 2 Then
        Call set(ip.mantissa, k, 1000000 * (git(num.mantissa, k) \ 1000000))
    ElseIf j = 3 Then
        Call set(ip.mantissa, k, 100000 * (git(num.mantissa, k) \ 100000))
    ElseIf j = 4 Then
        Call set(ip.mantissa, k, 10000 * (git(num.mantissa, k) \ 10000))
    ElseIf j = 5 Then
        Call set(ip.mantissa, k, 1000 * (git(num.mantissa, k) \ 1000))
    ElseIf j = 6 Then
        Call set(ip.mantissa, k, 100 * (git(num.mantissa, k) \ 100))
    ElseIf j = 7 Then
        Call set(ip.mantissa, k, 10 * (git(num.mantissa, k) \ 10))
    ElseIf j = 8 Then
        Call set(ip.mantissa, k, git(num.mantissa, k))
    End If
    ip.exponent = ex + BIAS
    ip.sign = num.sign
    result = ip
End Sub
 
Function fp2dbl# (n As decfloat)
    Dim As Long ex
    Dim As String v, f, ts
    If n.exponent > 0 Then
        ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
    Else
        ex = 0
    End If
    If n.sign Then v = "-" Else v = " "
    ts = Trim$(Str$(git(n.mantissa, 0)))
    If Len(ts) < 8 Then
        ts = ts + String$(8 - Len(ts), "0")
    End If
    v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
    ts = Trim$(Str$(git(n.mantissa, 1)))
    If Len(ts) < 8 Then
        ts = String$(8 - Len(ts), "0") + ts
    End If
    v = v + ts
 
    f = Str$(Abs(ex))
    f = String$(5 - Len(f), "0") + f
    If ex < 0 Then v = v + "E-" Else v = v + "E+"
    v = v + f
    fp2dbl# = Val(v)
End Function
 
Function fp2str_fix$ (n As decfloat, places As Long)
    Dim As Long i, ex
    Dim As String v, ts, s
    If n.exponent > 0 Then
        ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
    Else
        ex = 0
    End If
    If n.sign Then s = "-" Else s = " "
    ts = Trim$(Str$(git(n.mantissa, 0)))
    If Len(ts) < 8 Then
        ts = ts + String$(8 - Len(ts), "0")
    End If
    v = ts
    For i = 1 To NUM_DWORDS - 1
        ts = Trim$(Str$(git(n.mantissa, i)))
        If Len(ts) < 8 Then
            ts = String$(8 - Len(ts), "0") + ts
        End If
        v = v + ts
    Next
    If places < NUM_DIGITS Then
        v = Left$(v, places)
    End If
    If ex = 0 Then
        v = Left$(v, 1) + "." + Mid$(v, 2)
    ElseIf ex < 0 Then
        v = "0." + String$(Abs(ex) - 1, "0") + v
    ElseIf ex > 0 Then
        v = Left$(v, ex + 1) + "." + Mid$(v, ex + 2)
    End If
    fp2str_fix$ = s + v
End Function
 
Function fp2str$ (n As decfloat, places As Long)
    Dim As Long ex
    If n.exponent > 0 Then
        ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
    Else
        ex = 0
    End If
    If Abs(ex) < places Then
        fp2str = fp2str_fix(n, places)
    Else
        fp2str = fp2str_exp(n, places)
    End If
End Function
 
Sub str2fp (result As decfloat, value As String)
    Dim As Long j, s, d, e, ep, ex, es, i, f, fp, fln
    Dim As String c, f1, f2, f3, ts
    Dim As Unsigned Long ulng
    Dim n As decfloat
    j = 1
    s = 1
    d = 0
    e = 0
    ep = 0
    ex = 0
    es = 1
    i = 0
    f = 0
    fp = 0
    f1 = ""
    f2 = ""
    f3 = ""
    value = UCase$(value)
    fln = Len(value)
 
    While j <= fln
        c = Mid$(value, j, 1)
        If ep = 1 Then
            If c = " " Then
                j = j + 1
                GoTo skip_while
            End If
            If c = "-" Then
                es = -es
                c = ""
            End If
            If c = "+" Then
                j = j + 1
                GoTo skip_while
            End If
            If (c = "0") And (f3 = "") Then
                j = j + 1
                GoTo skip_while
            End If
            If (c > "/") And (c < ":") Then 'c is digit between 0 and 9
                f3 = f3 + c
                ex = 10 * ex + (Asc(c) - 48)
                j = j + 1
                GoTo skip_while
            End If
        End If
 
        If c = " " Then
            j = j + 1
            GoTo skip_while
        End If
        If c = "-" Then
            s = -s
            j = j + 1
            GoTo skip_while
        End If
        If c = "+" Then
            j = j + 1
            GoTo skip_while
        End If
        If c = "." Then
            If d = 1 Then
                j = j + 1
                GoTo skip_while
            End If
            d = 1
        End If
        If (c > "/") And (c < ":") Then 'c is digit between 0 and 9
            If ((c = "0") And (i = 0)) Then
                If d = 0 Then
                    j = j + 1
                    GoTo skip_while
                End If
                If (d = 1) And (f = 0) Then
                    e = e - 1
                    j = j + 1
                    GoTo skip_while
                End If
            End If
            If d = 0 Then
                f1 = f1 + c
                i = i + 1
            Else
                If (c > "0") Then
                    fp = 1
                End If
                f2 = f2 + c
                f = f + 1
            End If
        End If
        If c = "E" Or c = "D" Then
            ep = 1
        End If
        j = j + 1
        skip_while:
    Wend
    If fp = 0 Then
        f = 0
        f2 = ""
    End If
 
    If s = -1 Then s = &H8000 Else s = 0
    n.sign = s
    ex = es * ex - 1 + i + e
    f1 = f1 + f2
    f1 = Mid$(f1, 1, 1) + Right$(f1, Len(f1) - 1)
    fln = Len(f1)
    If Len(f1) > (NUM_DIGITS + 1 + 8) Then
        f1 = Mid$(f1, 1, (NUM_DIGITS + 1 + 8))
    End If
    While Len(f1) < (NUM_DIGITS + 1 + 8)
        f1 = f1 + "0"
    Wend
    j = 1
    For i = 0 To NUM_DWORDS
        ts = Mid$(f1, j, 8)
        ulng = Val(ts)
        Call set(n.mantissa, i, ulng)
        If ulng <> 0 Then fp = 1
        j = j + 8
    Next
    If fp Then n.exponent = (ex + BIAS + 1) Else n.exponent = 0
    result = n
End Sub
 
Sub si2fp (result As decfloat, m As Integer64, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat fac1
    Dim As Long i
    Dim As Integer64 n
    n = Abs(m)
 
    If n > 9999999999999999 Then
        Call str2fp(fac1, Str$(m))
        result = fac1: Exit Sub
    End If
 
    For i = 1 To digits
        Call set(fac1.mantissa, i, 0)
    Next
 
    If m = 0 Then
        fac1.exponent = 0
        fac1.sign = 0
        result = fac1: Exit Sub
    End If
 
    fac1.exponent = BIAS
    If n < 100000000 Then
        If n < 10 Then
            Call set(fac1.mantissa, 0, n * 10000000)
            fac1.exponent = fac1.exponent + 1
        ElseIf n < 100 Then
            Call set(fac1.mantissa, 0, n * 1000000)
            fac1.exponent = fac1.exponent + 2
        ElseIf n < 1000 Then
            Call set(fac1.mantissa, 0, n * 100000)
            fac1.exponent = fac1.exponent + 3
        ElseIf n < 10000 Then
            Call set(fac1.mantissa, 0, n * 10000)
            fac1.exponent = fac1.exponent + 4
        ElseIf n < 100000 Then
            Call set(fac1.mantissa, 0, n * 1000)
            fac1.exponent = fac1.exponent + 5
        ElseIf n < 1000000 Then
            Call set(fac1.mantissa, 0, n * 100)
            fac1.exponent = fac1.exponent + 6
        ElseIf n < 10000000 Then
            Call set(fac1.mantissa, 0, n * 10)
            fac1.exponent = fac1.exponent + 7
        ElseIf n < 100000000 Then
            Call set(fac1.mantissa, 0, n)
            fac1.exponent = fac1.exponent + 8
        End If
    End If
    If n > 99999999 Then
        fac1.exponent = fac1.exponent + 8
        If n < 1000000000 Then
            Call set(fac1.mantissa, 0, n \ 10)
            Call set(fac1.mantissa, 1, (n Mod 10) * 10000000)
            fac1.exponent = fac1.exponent + 1
        ElseIf n < 100000000000 Then
            Call set(fac1.mantissa, 0, n \ 100)
            Call set(fac1.mantissa, 1, (n Mod 100) * 1000000)
            fac1.exponent = fac1.exponent + 2
        ElseIf n < 1000000000000 Then
            Call set(fac1.mantissa, 0, n \ 1000)
            Call set(fac1.mantissa, 1, (n Mod 1000) * 100000)
            fac1.exponent = fac1.exponent + 3
        ElseIf n < 10000000000000 Then
            Call set(fac1.mantissa, 0, n \ 10000)
            Call set(fac1.mantissa, 1, (n Mod 10000) * 10000)
            fac1.exponent = fac1.exponent + 4
        ElseIf n < 100000000000000 Then
            Call set(fac1.mantissa, 0, n \ 100000)
            Call set(fac1.mantissa, 1, (n Mod 100000) * 1000)
            fac1.exponent = fac1.exponent + 5
        ElseIf n < 1000000000000000 Then
            Call set(fac1.mantissa, 0, n \ 1000000)
            Call set(fac1.mantissa, 1, (n Mod 1000000) * 100)
            fac1.exponent = fac1.exponent + 6
        ElseIf n < 10000000000000000 Then
            Call set(fac1.mantissa, 0, n \ 10000000)
            Call set(fac1.mantissa, 1, (n Mod 10000000) * 10)
            fac1.exponent = fac1.exponent + 7
        ElseIf n < 100000000000000000 Then
            Call set(fac1.mantissa, 0, n \ 100000000)
            Call set(fac1.mantissa, 1, n Mod 100000000)
            fac1.exponent = fac1.exponent + 8
        End If
    End If
    If m < 0 Then
        fac1.sign = &H8000
    Else
        fac1.sign = 0
    End If
    result = fac1
End Sub
 
Sub RSHIFT_1 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = digits To 1 Step -1
        v1 = git(mantissa.mantissa, i) \ 10
        v2 = git(mantissa.mantissa, i - 1) Mod 10
        v2 = v2 * 10000000 + v1
        Call set(mantissa.mantissa, i, v2)
    Next
    Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 10)
End Sub
 
Sub LSHIFT_1 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = 0 To digits - 1
        v1 = git(mantissa.mantissa, i) Mod 10000000
        v2 = git(mantissa.mantissa, i + 1) \ 10000000
        Call set(mantissa.mantissa, i, v1 * 10 + v2)
        Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 10000000)
    Next
    Call set(mantissa.mantissa, digits, 10 * (git(mantissa.mantissa, digits) Mod 10000000))
End Sub
 
Sub RSHIFT_2 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = digits To 1 Step -1
        v1 = git(mantissa.mantissa, i) \ 100
        v2 = git(mantissa.mantissa, i - 1) Mod 100
        v2 = v2 * 1000000 + v1
        Call set(mantissa.mantissa, i, v2)
    Next
    Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 100)
End Sub
 
Sub LSHIFT_2 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = 0 To digits - 1
        v1 = git(mantissa.mantissa, i) Mod 1000000
        v2 = git(mantissa.mantissa, i + 1) \ 1000000
        Call set(mantissa.mantissa, i, v1 * 100 + v2)
        Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 1000000)
    Next
    Call set(mantissa.mantissa, digits, 100 * (git(mantissa.mantissa, digits) Mod 1000000))
End Sub
 
Sub RSHIFT_3 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = digits To 1 Step -1
        v1 = git(mantissa.mantissa, i) \ 1000
        v2 = git(mantissa.mantissa, i - 1) Mod 1000
        v2 = v2 * 100000 + v1
        Call set(mantissa.mantissa, i, v2)
    Next
    Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 1000)
End Sub
 
Sub LSHIFT_3 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = 0 To digits - 1
        v1 = git(mantissa.mantissa, i) Mod 100000
        v2 = git(mantissa.mantissa, i + 1) \ 100000
        Call set(mantissa.mantissa, i, v1 * 1000 + v2)
        Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 100000)
    Next
    Call set(mantissa.mantissa, digits, 1000 * (git(mantissa.mantissa, digits) Mod 100000))
End Sub
 
Sub RSHIFT_4 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = digits To 1 Step -1
        v1 = git(mantissa.mantissa, i) \ 10000
        v2 = git(mantissa.mantissa, i - 1) Mod 10000
        v2 = v2 * 10000 + v1
        Call set(mantissa.mantissa, i, v2)
    Next
    Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 10000)
End Sub
 
Sub LSHIFT_4 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = 0 To digits - 1
        v1 = git(mantissa.mantissa, i) Mod 10000
        v2 = git(mantissa.mantissa, i + 1) \ 10000
        Call set(mantissa.mantissa, i, v1 * 10000 + v2)
        Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 10000)
    Next
    Call set(mantissa.mantissa, digits, 10000 * (git(mantissa.mantissa, digits) Mod 10000))
End Sub
 
Sub RSHIFT_5 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = digits To 1 Step -1
        v1 = git(mantissa.mantissa, i) \ 100000
        v2 = git(mantissa.mantissa, i - 1) Mod 100000
        v2 = v2 * 1000 + v1
        Call set(mantissa.mantissa, i, v2)
    Next
    Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 100000)
End Sub
 
Sub LSHIFT_5 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = 0 To digits - 1
        v1 = git(mantissa.mantissa, i) Mod 1000
        v2 = git(mantissa.mantissa, i + 1) \ 1000
        Call set(mantissa.mantissa, i, v1 * 100000 + v2)
        Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 1000)
    Next
    Call set(mantissa.mantissa, digits, 100000 * (git(mantissa.mantissa, digits) Mod 1000))
End Sub
 
Sub RSHIFT_6 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = digits To 1 Step -1
        v1 = git(mantissa.mantissa, i) \ 1000000
        v2 = git(mantissa.mantissa, i - 1) Mod 1000000
        v2 = v2 * 100 + v1
        Call set(mantissa.mantissa, i, v2)
    Next
    Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 1000000)
End Sub
 
Sub LSHIFT_6 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = 0 To digits - 1
        v1 = git(mantissa.mantissa, i) Mod 100
        v2 = git(mantissa.mantissa, i + 1) \ 100
        Call set(mantissa.mantissa, i, v1 * 1000000 + v2)
        Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 100)
    Next
    Call set(mantissa.mantissa, digits, 1000000 * (git(mantissa.mantissa, digits) Mod 100))
End Sub
 
Sub RSHIFT_7 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = digits To 1 Step -1
        v1 = git(mantissa.mantissa, i) \ 10000000
        v2 = git(mantissa.mantissa, i - 1) Mod 10000000
        v2 = v2 * 10 + v1
        Call set(mantissa.mantissa, i, v2)
    Next
    Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 10000000)
End Sub
 
Sub LSHIFT_7 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = 0 To digits - 1
        v1 = git(mantissa.mantissa, i) Mod 10
        v2 = git(mantissa.mantissa, i + 1) \ 10
        Call set(mantissa.mantissa, i, v1 * 10000000 + v2)
        Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 10)
    Next
    Call set(mantissa.mantissa, digits, 10000000 * (git(mantissa.mantissa, digits) Mod 10))
End Sub
 
Sub RSHIFT_8 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As Long i
    For i = digits To 1 Step -1
        Call set(mantissa.mantissa, i, git(mantissa.mantissa, i - 1))
    Next
    Call set(mantissa.mantissa, 0, 0)
End Sub
 
Sub LSHIFT_8 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As Long i
    For i = 0 To digits - 1
        Call set(mantissa.mantissa, i, git(mantissa.mantissa, i + 1))
    Next
    Call set(mantissa.mantissa, digits, 0)
End Sub
 
Function fpcmp& (x As decfloat, y As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As Long c, i
    If x.sign < y.sign Then fpcmp = -1: Exit Function
    If x.sign > y.sign Then fpcmp = 1: Exit Function
    If x.sign = y.sign Then
        If x.exponent = y.exponent Then
            For i = 0 To digits
                c = git(x.mantissa, i) - git(y.mantissa, i)
                If c <> 0 Then Exit For
            Next
            If c < 0 Then fpcmp = -1: Exit Function
            If c = 0 Then fpcmp = 0: Exit Function
            If c > 0 Then fpcmp = 1: Exit Function
        End If
        If x.exponent < y.exponent Then fpcmp = -1: Exit Function
        If x.exponent > y.exponent Then fpcmp = 1: Exit Function
    End If
    ' if we reach this point it means that the signs are different
    ' and if the sign of x is set meaning that x is negative then x < y
 
End Function
 
Sub NORM_FAC1 (fac1 As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    ' normalize the number in fac1
    ' all routines exit through this one.
 
    'see if the mantissa is all zeros.
    'if so, set the exponent and sign equal to 0.
    Dim As Long i, er, f
    er = 0: f = 0
    For i = 0 To digits
        If git(fac1.mantissa, i) > 0 Then f = 1
    Next
    If f = 0 Then
        fac1.exponent = 0
        fac1.sign = 0
        Exit Sub
        'if the highmost Digit in fac1_man is nonzero,
        'shift the mantissa right 1 Digit and
        'increment the exponent
    ElseIf git(fac1.mantissa, 0) > 99999999 Then
        Call RSHIFT_1(fac1, digits)
        fac1.exponent = fac1.exponent + 1
    Else
        'now shift fac1_man 1 to the left until a
        'nonzero digit appears in the next-to-highest
        'Digit of fac1_man.  decrement exponent for
        'each shift.
        While git(fac1.mantissa, 0) = 0
            Call LSHIFT_8(fac1, digits)
            fac1.exponent = fac1.exponent - 8
            If fac1.exponent = 0 Then
                Print "NORM_FAC1=EXPU_ERR"
                Exit Sub
            End If
        Wend
        If git(fac1.mantissa, 0) < 10 Then
            Call LSHIFT_7(fac1, digits)
            fac1.exponent = fac1.exponent - 7
        ElseIf git(fac1.mantissa, 0) < 100 Then
            Call LSHIFT_6(fac1, digits)
            fac1.exponent = fac1.exponent - 6
        ElseIf git(fac1.mantissa, 0) < 1000 Then
            Call LSHIFT_5(fac1, digits)
            fac1.exponent = fac1.exponent - 5
        ElseIf git(fac1.mantissa, 0) < 10000 Then
            Call LSHIFT_4(fac1, digits)
            fac1.exponent = fac1.exponent - 4
        ElseIf git(fac1.mantissa, 0) < 100000 Then
            Call LSHIFT_3(fac1, digits)
            fac1.exponent = fac1.exponent - 3
        ElseIf git(fac1.mantissa, 0) < 1000000 Then
            Call LSHIFT_2(fac1, digits)
            fac1.exponent = fac1.exponent - 2
        ElseIf git(fac1.mantissa, 0) < 10000000 Then
            Call LSHIFT_1(fac1, digits)
            fac1.exponent = fac1.exponent - 1
        End If
    End If
    'check for overflow/underflow
    If fac1.exponent < 0 Then
        Print "NORM_FAC1=EXPO_ERR"
    End If
End Sub
 
Sub fpadd_aux (fac1 As decfloat, fac2 As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As Unsigned Long v, c, i
    c = 0
    For i = digits To 1 Step -1
        v = git(fac2.mantissa, i) + git(fac1.mantissa, i) + c
        If v > 99999999 Then
            v = v - 100000000
            c = 1
        Else
            c = 0
        End If
        Call set(fac1.mantissa, i, v)
    Next
    v = git(fac1.mantissa, 0) + git(fac2.mantissa, 0) + c
    Call set(fac1.mantissa, 0, v)
 
    Call NORM_FAC1(fac1, digits)
 
End Sub
 
Sub fpsub_aux (fac1 As decfloat, fac2 As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As Long v, c, i
    c = 0
    For i = digits To 1 Step -1
        v = git(fac1.mantissa, i) - git(fac2.mantissa, i) - c
        If v < 0 Then
            v = v + 100000000
            c = 1
        Else
            c = 0
        End If
        Call set(fac1.mantissa, i, v)
    Next
    v = git(fac1.mantissa, 0) - git(fac2.mantissa, 0) - c
    Call set(fac1.mantissa, 0, v)
 
    Call NORM_FAC1(fac1, digits)
End Sub
 
Sub fpadd (result As decfloat, x As decfloat, y As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As decfloat fac1, fac2
    Dim As Long i, t, c, xsign, ysign
 
    xsign = x.sign: x.sign = 0
    ysign = y.sign: y.sign = 0
    c = fpcmp(x, y, digits)
 
    x.sign = xsign
    y.sign = ysign
    If c < 0 Then
        fac1 = y
        fac2 = x
    Else
        fac1 = x
        fac2 = y
    End If
    t = fac1.exponent - fac2.exponent
    t = ((fac1.exponent And &H7FFFFFFF) - BIAS - 1) - ((fac2.exponent And &H7FFFFFFF) - BIAS - 1)
 
    If t < (NUM_DIGITS + 8) Then
        'The difference between the two
        'exponents indicate how many times
        'we have to multiply the mantissa
        'of FAC2 by 10 (i.e., shift it right 1 place).
        'If we have to shift more times than
        'we have digits, the result is already in FAC1.
        t = fac1.exponent - fac2.exponent
        If t > 0 And t < (NUM_DIGITS + 8) Then 'shift
 
            i = t \ 8
            While i > 0
                Call RSHIFT_8(fac2, digits)
                t = t - 8
                i = i - 1
            Wend
            If t = 7 Then
                Call RSHIFT_7(fac2, digits)
            ElseIf t = 6 Then
                Call RSHIFT_6(fac2, digits)
            ElseIf t = 5 Then
                Call RSHIFT_5(fac2, digits)
            ElseIf t = 4 Then
                Call RSHIFT_4(fac2, digits)
            ElseIf t = 3 Then
                Call RSHIFT_3(fac2, digits)
            ElseIf t = 2 Then
                Call RSHIFT_2(fac2, digits)
            ElseIf t = 1 Then
                Call RSHIFT_1(fac2, digits)
            End If
        End If
        'See if the signs of the two numbers
        'are the same.  If so, add; if not, subtract.
        If fac1.sign = fac2.sign Then 'add
            Call fpadd_aux(fac1, fac2, digits)
        Else
            Call fpsub_aux(fac1, fac2, digits)
        End If
    End If
    result = fac1
End Sub
 
Sub fpsub (result As decfloat, x As decfloat, y As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat fac1, fac2
    fac1 = x
    fac2 = y
    fac2.sign = fac2.sign Xor &H8000
    Call fpadd(fac1, fac1, fac2, digits)
    result = fac1
End Sub
 
Sub fpmul (result As decfloat, x As decfloat, y As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat fac1, fac2
    Dim As Long i, j, ex, er, den, num
    Dim As Integer64 digit, carry
    Dim As Unsigned Integer64 fac3(0 To digits)
 
    fac1 = x
    fac2 = y
    'check exponents.  if either is zero,
    'the result is zero
    If fac1.exponent = 0 Or fac2.exponent = 0 Then 'result is zero...clear fac1.
        fac1.sign = 0
        fac1.exponent = 0
        For i = 0 To digits
            Call set(fac1.mantissa, i, 0)
        Next
        'NORM_FAC1(fac1)
        result = fac1: Exit Sub
    Else
 
        If ex < 0 Then
            er = EXPO_ERR
            result = fac1: Exit Sub
        End If
 
        'clear fac3 mantissa
        For i = 0 To digits
            fac3(i) = 0
        Next
 
        den = digits
        While git(fac2.mantissa, den) = 0
            den = den - 1
        Wend
        num = digits
        While git(fac1.mantissa, num) = 0
            num = num - 1
        Wend
 
        If num < den Then
            Swap fac1, fac2
            'fac1=y
            'fac2=x
            Swap den, num
        End If
 
        For j = den To 0 Step -1
            carry = 0
            digit = git(fac2.mantissa, j)
            For i = num To 0 Step -1
                fac3(i) = fac3(i) + digit * git(fac1.mantissa, i)
            Next
            For i = num To 0 Step -1
                fac3(i) = fac3(i) + carry
                carry = fac3(i) \ 100000000
                fac3(i) = (fac3(i) Mod 100000000)
            Next
 
            For i = digits To 1 Step -1
                fac3(i) = fac3(i - 1)
            Next
            fac3(0) = carry
        Next
 
        For i = 0 To digits
            Call set(fac1.mantissa, i, fac3(i))
        Next
    End If
    'now determine exponent of result.
    'as you do...watch for overflow.
    ex = fac2.exponent - BIAS + fac1.exponent
    fac1.exponent = ex
    'determine the sign of the product
    fac1.sign = fac1.sign Xor fac2.sign
    Call NORM_FAC1(fac1, digits)
    result = fac1
End Sub
 
Sub fpmul_si (result As decfloat, x As decfloat, y As Integer64, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat fac1, fac2
    Dim As Long count, ex, er, i
    Dim As Integer64 carry, digit, prod, value
    fac1 = x
    digit = Abs(y)
    If digit > 99999999 Then
        Call si2fp(fac2, y, digits)
        Call fpmul(fac1, fac1, fac2, digits)
        result = fac1: Exit Sub
    End If
    'check exponents.  if either is zero,
    'the result is zero
    If fac1.exponent = 0 Or y = 0 Then 'result is zero...clear fac1.
        fac1.sign = 0
        fac1.exponent = 0
        For count = 0 To digits
            Call set(fac1.mantissa, count, 0)
        Next
        Call NORM_FAC1(fac1, digits)
        result = fac1: Exit Sub
    Else
        If digit = 1 Then
            If y < 0 Then
                fac1.sign = fac1.sign Xor &H8000
            End If
            result = fac1: Exit Sub
        End If
        'now determine exponent of result.
        'as you do...watch for overflow.
 
        If ex < 0 Then
            er = EXPO_ERR
            result = fac1: Exit Sub
        End If
 
        carry = 0
 
        For i = digits To 0 Step -1
            prod = digit * git(fac1.mantissa, i) + carry
            value = (prod Mod 100000000)
            Call set(fac1.mantissa, i, value)
            carry = prod \ 100000000
        Next
 
        If carry < 10 Then
            Call RSHIFT_1(fac1, digits)
            fac1.exponent = fac1.exponent + 1
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10000000)
        ElseIf carry < 100 Then
            Call RSHIFT_2(fac1, digits)
            fac1.exponent = fac1.exponent + 2
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 1000000)
        ElseIf carry < 1000 Then
            Call RSHIFT_3(fac1, digits)
            fac1.exponent = fac1.exponent + 3
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 100000)
        ElseIf carry < 10000 Then
            Call RSHIFT_4(fac1, digits)
            fac1.exponent = fac1.exponent + 4
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10000)
        ElseIf carry < 100000 Then
            Call RSHIFT_5(fac1, digits)
            fac1.exponent = fac1.exponent + 5
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 1000)
        ElseIf carry < 1000000 Then
            Call RSHIFT_6(fac1, digits)
            fac1.exponent = fac1.exponent + 6
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 100)
        ElseIf carry < 10000000 Then
            Call RSHIFT_7(fac1, digits)
            fac1.exponent = fac1.exponent + 7
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10)
        ElseIf carry < 100000000 Then
            Call RSHIFT_8(fac1, digits)
            fac1.exponent = fac1.exponent + 8
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry)
        End If
 
    End If
 
    Call NORM_FAC1(fac1, digits)
 
    If y < 0 Then
        fac1.sign = fac1.sign Xor &H8000
    End If
    result = fac1
End Sub
 
Sub fpmul_ll (result As decfloat, x As decfloat, y As Integer64, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat fac1, fac2
    Dim As Long count, ex, er, i
    Dim As Integer64 carry, digit, prod, value, n0, n1
    fac1 = x
    digit = Abs(y)
    If digit > 99999999 And digit < 10000000000000000 Then
        n0 = digit \ 100000000
        n1 = digit Mod 100000000
        Call fpmul_si(fac2, fac1, n1, digits)
        fac1.exponent = fac1.exponent + 8
        Call fpmul_si(fac1, fac1, n0, digits)
        Call fpadd(fac1, fac1, fac2, digits)
        If y < 0 Then fac1.sign = fac1.sign Xor &H8000
        result = fac1: Exit Sub
    End If
    If digit > 9999999999999999 Then
        Call str2fp(fac2, Str$(y))
        Call fpmul(fac1, fac1, fac2, digits)
        result = fac1: Exit Sub
    End If
    Call si2fp(fac2, y, digits)
 
    'check exponents.  if either is zero,
    'the result is zero
    If fac1.exponent = 0 Or y = 0 Then 'result is zero...clear fac1.
        fac1.sign = 0
        fac1.exponent = 0
        For count = 0 To digits
            Call set(fac1.mantissa, count, 0)
        Next
        Call NORM_FAC1(fac1, digits)
        result = fac1: Exit Sub
    Else
        If digit = 1 Then
            If y < 0 Then
                fac1.sign = fac1.sign Xor &H8000
            End If
            result = fac1: Exit Sub
        End If
        'now determine exponent of result.
        'as you do...watch for overflow.
 
        If ex < 0 Then
            er = EXPO_ERR
            result = fac1: Exit Sub
        End If
 
        carry = 0
 
        For i = digits To 0 Step -1
            prod = digit * git(fac1.mantissa, i) + carry
            value = (prod Mod 100000000)
            Call set(fac1.mantissa, i, value)
            carry = prod \ 100000000
        Next
 
        If carry < 10 Then
            Call RSHIFT_1(fac1, digits)
            fac1.exponent = fac1.exponent + 1
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10000000)
        ElseIf carry < 100 Then
            Call RSHIFT_2(fac1, digits)
            fac1.exponent = fac1.exponent + 2
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 1000000)
        ElseIf carry < 1000 Then
            Call RSHIFT_3(fac1, digits)
            fac1.exponent = fac1.exponent + 3
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 100000)
        ElseIf carry < 10000 Then
            Call RSHIFT_4(fac1, digits)
            fac1.exponent = fac1.exponent + 4
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10000)
        ElseIf carry < 100000 Then
            Call RSHIFT_5(fac1, digits)
            fac1.exponent = fac1.exponent + 5
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 1000)
        ElseIf carry < 1000000 Then
            Call RSHIFT_6(fac1, digits)
            fac1.exponent = fac1.exponent + 6
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 100)
        ElseIf carry < 10000000 Then
            Call RSHIFT_7(fac1, digits)
            fac1.exponent = fac1.exponent + 7
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10)
        ElseIf carry < 100000000 Then
            Call RSHIFT_8(fac1, digits)
            fac1.exponent = fac1.exponent + 8
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry)
        End If
 
    End If
 
    Call NORM_FAC1(fac1, digits)
 
    If y < 0 Then
        fac1.sign = fac1.sign Xor &H8000
    End If
    result = fac1
End Sub
 
Sub fpinv (result As decfloat, m As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As Double x
    Dim As Long k, l, ex
    Dim As Long prec
    Dim As decfloat r, r2, one, n
    n = m
    l = Log((NUM_DIGITS + 8) * 0.0625) * 1.5
    Dim As String v, ts
    If n.exponent > 0 Then
        ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
    Else
        ex = 0
    End If
 
    If n.sign Then v = "-" Else v = " "
    ts = Str$(git(n.mantissa, 0))
    If Len(ts) < 8 Then
        ts = ts + String$(8 - Len(ts), "0")
    End If
    v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
    ts = Str$(git(n.mantissa, 1))
    If Len(ts) < 8 Then
        ts = String$(8 - Len(ts), "0") + ts
    End If
    v = v + ts
    x = Val(v)
    If x = 0 Then Print "Div 0": result = r: Exit Sub
    If x = 1 And ex = 0 Then
        Call str2fp(r, "1")
        result = r: Exit Sub
    ElseIf x = 1 Then
        x = 10
        ex = ex - 1
    End If
    ex = (-1) - ex
    x = 1 / x
    Call str2fp(r, Str$(x))
    r.exponent = ex + BIAS + 1
    Call set(one.mantissa, 0, 10000000)
    one.exponent = BIAS + 1
    r2 = r
    prec = 3
    For k = 1 To l
        prec = 2 * prec - 1
        Call fpmul(r2, n, r, prec)
        Call fpsub(r2, one, r2, prec)
        Call fpmul(r2, r, r2, prec)
        Call fpadd(r, r, r2, prec)
    Next
    result = r
End Sub
 
Function min& (a As Long, b As Long)
    If a < b Then min& = a Else min& = b
End Function
 
Function RealW# (w() As Double, j As Long)
    Dim wx As Double
    wx = ((w(j - 1) * 10000 + w(j)) * 10000 + w(j + 1)) * 10000
    If UBound(w) >= (j + 2) Then wx = wx + w(j + 2)
    RealW# = wx
End Function
 
Sub subtract (w() As Double, q As Long, d() As Double, ka As Long, kb As Long)
    Dim As Long j
    For j = ka To kb
        w(j) = w(j) - q * d(j - ka + 2)
    Next
End Sub
 
Sub normalize (w() As Double, ka As Long, q As Long)
    w(ka) = w(ka) + w(ka - 1) * 10000
    w(ka - 1) = q
End Sub
 
Sub finalnorm (w() As Double, kb As Long)
    Dim As Long carry, j
    For j = kb To 3 Step -1
        If w(j) < 0 Then
            carry = ((-w(j) - 1) \ 10000) + 1
        Else
            If w(j) >= 10000 Then
                carry = -(w(j) \ 10000)
            Else
                carry = 0
            End If
        End If
        w(j) = w(j) + carry * 10000
        w(j - 1) = w(j - 1) - carry
    Next
End Sub
 
Sub fpdiv (result As decfloat, x As decfloat, y As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat fac1, fac2
    Dim As Long i, er, is_power_of_ten
 
    fac1 = x
    fac2 = y
 
    If fac2.exponent = 0 Then ' if fac2 = 0, return
        ' a divide-by-zero error and
        ' bail out.
        For i = 0 To digits
            Call set(fac1.mantissa, i, 99999999)
        Next
        fac1.exponent = 99999 + BIAS + 1
        er = DIVZ_ERR
        result = fac1
        Exit Sub
    ElseIf fac1.exponent = 0 Then 'fact1=0, just return
        er = 0
        result = fac1
        Exit Sub
    Else
        'check to see if fac2 is a power of ten
        is_power_of_ten = 0
        If git(fac2.mantissa, 0) = 10000000 Then
            is_power_of_ten = 1
            For i = 1 To digits
                If git(fac2.mantissa, i) <> 0 Then
                    is_power_of_ten = 0
                    Exit For
                End If
            Next
        End If
        'if fac2 is a power of ten then all we need to do is to adjust the sign and exponent and we are finished
        If is_power_of_ten = 1 Then
            fac1.sign = fac1.sign Xor fac2.sign
            fac1.exponent = fac1.exponent - fac2.exponent + BIAS + 1
            result = fac1
            Exit Sub
        End If
 
        Dim As Double result(1 To 2 * digits + 3), n(1 To 2 * digits + 3), d(1 To 2 * digits + 3)
        Const b = 10000
        Dim As Long j, last, laststep, q, t
        Dim As Long stp
        Dim As Double xd, xn, rund
        Dim As Double w(1 To UBound(n) + 4)
 
        For j = 0 To digits
            n(2 * j + 2) = git(fac1.mantissa, j) \ 10000
            n(2 * j + 3) = git(fac1.mantissa, j) Mod 10000
            d(2 * j + 2) = git(fac2.mantissa, j) \ 10000
            d(2 * j + 3) = git(fac2.mantissa, j) Mod 10000
        Next
        n(1) = (fac1.exponent And &H7FFFFFFF) - BIAS - 1
        d(1) = (fac2.exponent And &H7FFFFFFF) - BIAS - 1
        For j = UBound(n) To UBound(w)
            w(j) = 0
        Next
        t = UBound(n) - 1
        w(1) = n(1) - d(1) + 1
        w(2) = 0
        For j = 2 To UBound(n)
            w(j + 1) = n(j)
        Next
        xd = (d(2) * b + d(3)) * b + d(4) + d(5) / b
        laststep = t + 2
        For stp = 1 To laststep
            xn = RealW(w(), (stp + 2))
            q = Int(xn / xd)
            last = min(stp + t + 1, UBound(w))
            Call subtract(w(), q, d(), (stp + 2), last)
            Call normalize(w(), (stp + 2), q)
        Next
        Call finalnorm(w(), (laststep + 1))
        If w(2) <> 0 Then laststep = laststep - 1
        rund = w(laststep + 1) / b
        If rund >= 0.5 Then w(laststep) = w(laststep) + 1
        If w(2) = 0 Then
            For j = 1 To t + 1
                result(j) = w(j + 1)
            Next
        Else
            For j = 1 To t + 1
                result(j) = w(j)
            Next
        End If
        If w(2) = 0 Then result(1) = w(1) - 1 Else result(1) = w(1)
 
        For j = 0 To digits
            Call set(fac1.mantissa, j, result(2 * j + 2) * 10000 + result(2 * j + 3))
        Next
        Call NORM_FAC1(fac1, digits)
        fac1.exponent = (result(1) + BIAS)
    End If
    fac1.sign = fac1.sign Xor fac2.sign
    result = fac1
End Sub
 
' sqrt(num)
Sub fpsqr (result As decfloat, num As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat r, r2, tmp, n, half
    Dim As Long ex, k, l, prec
    Dim As String ts, v
    Dim As Double x
    l = Log((NUM_DIGITS + 8) * 0.0625) * 1.5
    'l=estimated number of iterations needed
    'first estimate is accurate to about 16 digits
    'l is approximatly = to log2((NUM_DIGITS+9)/16)
    'NUM_DIGITS+9 because decfloat has an extra 9 guard digits
    n = num
    Call si2fp(tmp, 0, digits)
    If fpcmp(n, tmp, digits) = 0 Then
        Call si2fp(r, 0, digits)
        result = r
        Exit Sub
    End If
    Call si2fp(tmp, 1, digits)
    If fpcmp(n, tmp, digits) = 0 Then
        Call si2fp(r, 1, digits)
        result = r
        Exit Sub
    End If
    Call si2fp(tmp, 0, digits)
    If fpcmp(n, tmp, digits) < 0 Then
        Call si2fp(r, 0, digits)
        result = r
        Exit Sub
    End If
    '=====================================================================
    'hack to bypass the limitation of double exponent range
    'in case the number is larger than what a double can handle
    'for example, if the number is 2e500
    'we separate the exponent and mantissa in this case 2
    'if the exponent is odd then multiply the mantissa by 10
    'take the square root and assign it to decfloat
    'divide the exponent in half for square root
    'in this case 1.414213562373095e250
    If n.exponent > 0 Then
        ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
    Else
        ex = 0
    End If
    ts = Trim$(Str$(git(n.mantissa, 0)))
    If Len(ts) < 8 Then
        ts = ts + String$(8 - Len(ts), "0")
    End If
    v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
    ts = Trim$(Str$(git(n.mantissa, 1)))
    If Len(ts) < 8 Then
        ts = String$(8 - Len(ts), "0") + ts
    End If
    v = v + ts
    x = Val(v)
    If x = 0 Then Print "Div 0": result = r: Exit Sub
    If x = 1 And ex = 0 Then
        Call si2fp(r, 1, digits)
        result = r
        Exit Sub
    End If
    If Abs(ex) And 1 Then
        x = x * 10
        ex = ex - 1
    End If
    x = Sqr(x) 'approximation
    v = Trim$(Str$(x))
    k = InStr(v, ".")
    Call str2fp(r, v)
    r.exponent = ex \ 2 + BIAS + 1
    If Len(v) > 1 And k = 0 Then r.exponent = r.exponent + 1
    For k = 1 To digits
        Call set(half.mantissa, k, 0)
    Next
    Call set(half.mantissa, 0, 50000000)
    half.exponent = BIAS
    half.sign = 0
    '=====================================================================
    'Newton-Raphson method
    prec = 3
    For k = 1 To l + 1
        prec = 2 * prec - 1
        Call fpdiv(tmp, n, r, prec)
        Call fpadd(r2, r, tmp, prec)
        Call fpmul(r, r2, half, prec)
    Next
    result = r
End Sub
 
Sub fpdiv_si (result As decfloat, num As decfloat, den As Long, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat fac1
    Dim As Unsigned Integer64 carry, remder
    Dim As Integer64 i, divisor
    Dim As Integer64 quotient
    remder = 0
    divisor = Abs(den)
    fac1 = num
    If divisor = 0 Then
        Print "error: divisor = 0"
        result = fac1: Exit Sub
    End If
    If divisor > 99999999 Then
        Print "error: divisor too large"
        result = fac1: Exit Sub
    End If
 
    For i = 0 To digits
        quotient = git(fac1.mantissa, i) + remder * 100000000
        remder = quotient Mod divisor
        Call set(fac1.mantissa, i, quotient \ divisor)
    Next
    quotient = remder * 100000000
    quotient = quotient \ divisor
    carry = git(fac1.mantissa, 0)
 
    If carry = 0 Then
        Call LSHIFT_8(fac1, digits)
        fac1.exponent = fac1.exponent - 8
        Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient)
    ElseIf carry < 10 Then
        Call LSHIFT_7(fac1, digits)
        fac1.exponent = fac1.exponent - 7
        Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 10)
    ElseIf carry < 100 Then
        Call LSHIFT_6(fac1, digits)
        fac1.exponent = fac1.exponent - 6
        Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 100)
    ElseIf carry < 1000 Then
        Call LSHIFT_5(fac1, digits)
        fac1.exponent = fac1.exponent - 5
        Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 1000)
    ElseIf carry < 10000 Then
        Call LSHIFT_4(fac1, digits)
        fac1.exponent = fac1.exponent - 4
        Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 10000)
    ElseIf carry < 100000 Then
        Call LSHIFT_3(fac1, digits)
        fac1.exponent = fac1.exponent - 3
        Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 100000)
    ElseIf carry < 1000000 Then
        Call LSHIFT_2(fac1, digits)
        fac1.exponent = fac1.exponent - 2
        Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 1000000)
    ElseIf carry < 10000000 Then
        Call LSHIFT_1(fac1, digits)
        fac1.exponent = fac1.exponent - 1
        Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 10000000)
    End If
 
    'NORM_FAC1(fac1)
    If den < 0 Then
        fac1.sign = fac1.sign Xor &H8000
    End If
    result = fac1
End Sub
 
'fractional part of num
Sub fpfrac (result As decfloat, num As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat n
    Call fpfix(n, num)
    Call fpsub(n, num, n, digits)
    result = n
End Sub
 
'returns the positive of n
Sub fpabs (result As decfloat, n As decfloat)
    Dim As decfloat x
    x = n
    x.sign = 0
    result = x
End Sub
 
'changes the sign of n, if n is positive then n will be negative & vice versa
Sub fpneg (result As decfloat, n As decfloat)
    Dim As decfloat x
    x = n
    x.sign = x.sign Xor &H8000
    result = x
End Sub
 
'returns the negative of n regardless of the sign of n
Sub fpnegative (result As decfloat, n As decfloat)
    Dim As decfloat x
    x = n
    x.sign = &H8000
    result = x
End Sub
 
Sub fpfmod (quotient As decfloat, f_mod As decfloat, num As decfloat, denom As decfloat, digits As Long)
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat q, fm 'make copy in case the destination and source are the same
    Call fpdiv(fm, num, denom, digits): Call fpfix(q, fm)
    Call fpsub(fm, fm, q, digits): Call fpmul(fm, fm, denom, digits)
    quotient = q
    f_mod = fm
End Sub
 
Sub fpeps (result As decfloat, digits As Long)
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat ep
    Call si2fp(ep, 1, digits)
    ep.exponent = (-(NUM_DIGITS) + BIAS + 1)
    result = ep
End Sub
 
Sub fpipow (result As decfloat, x As decfloat, e As Integer64, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    'take x to an Long power
    Dim As decfloat y
    Dim As decfloat z
    Dim As Integer64 n, c
    Dim As Long i
    c = 0
    y = x
    n = Abs(e)
    z.sign = 0
    z.exponent = (BIAS + 1)
    Call set(z.mantissa, 0, 10000000)
    For i = 1 To NUM_DWORDS
        Call set(z.mantissa, i, 0)
    Next
    While n > 0
        While (n And 1) = 0
            n = n \ 2
            Call fpmul(y, y, y, digits)
            c = c + 1
        Wend
        n = n - 1
        Call fpmul(z, y, z, digits)
        c = c + 1
    Wend
    If e < 0 Then
        Call fpinv(z, z, digits)
    End If
    result = z
End Sub
 
Sub fpfactorial (result As decfloat, n As Long, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    Dim As Long i
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat f
    If n < 0 Then Print "inf": result = f: Exit Sub
    If n = 0 Or n = 1 Then
        Call si2fp(f, 1, digits)
        result = f: Exit Sub
    End If
    Call si2fp(f, 2, digits)
    If n = 2 Then result = f: Exit Sub
    For i = 3 To n
        Call fpmul_si(f, f, i, digits)
    Next
    result = f
End Sub
 
Sub fplogTaylor (result As decfloat, x As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    'taylor series
    '====================Log Guard==================
    Dim As decfloat g, zero
    Dim As Long i
    Call fpabs(g, x)
    zero.sign = 0
    zero.exponent = 0
    For i = 0 To NUM_DWORDS
        Call set(zero.mantissa, i, 0)
    Next
    If fpcmp(g, x, digits) <> 0 Then result = zero: Exit Sub
    If fpcmp(x, zero, digits) = 0 Then result = zero: Exit Sub
    '=============================================
    Dim As Long invflag
    Dim As decfloat XX, Term, Accum, x9, tmp, tmp2
    Dim As decfloat T, B, one, Q, two
 
    one.sign = 0
    one.exponent = (BIAS + 1)
    Call set(one.mantissa, 0, 10000000)
    two.sign = 0
    two.exponent = (BIAS + 1)
    Call set(two.mantissa, 0, 20000000)
    For i = 1 To NUM_DWORDS
        Call set(one.mantissa, i, 0)
        Call set(two.mantissa, i, 0)
    Next
    x9 = x
    If fpcmp(x, one, digits) < 0 Then
        invflag = 1
        Call fpdiv(x9, one, x9, digits)
    End If
    Call fpsub(T, x9, one, digits)
    Call fpadd(B, x9, one, digits)
    Call fpdiv(Accum, T, B, digits)
    Call fpdiv(Q, T, B, digits)
    tmp = Q
    Call fpmul(XX, Q, Q, digits)
    Dim As Long c
    c = 1
    Do
        c = c + 2
        tmp2 = tmp
        Call fpmul(Q, Q, XX, digits)
        Call fpdiv_si(Term, Q, c, digits)
        Call fpadd(Accum, tmp, Term, digits)
        Swap tmp, Accum
    Loop Until fpcmp(tmp, tmp2, digits) = 0
    Call fpmul_si(Accum, Accum, 2, digits)
    If invflag Then
        Call fpneg(Accum, Accum)
        result = Accum: Exit Sub
    End If
    result = Accum
End Sub
 
Sub fplog (result As decfloat, x As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    '====================Log Guard==================
    Dim As decfloat g, one, zero
    Dim As Long i, factor
    zero.sign = 0
    zero.exponent = 0
    For i = 0 To NUM_DWORDS
        Call set(zero.mantissa, i, 0)
    Next
    Call si2fp(one, 1, digits)
    Call fpabs(g, x)
    If fpcmp(g, x, digits) <> 0 Then result = zero: Exit Sub
    If fpcmp(x, zero, digits) = 0 Then result = zero: Exit Sub
    If fpcmp(x, one, digits) = 0 Then result = zero: Exit Sub
    '=============================================
    Dim As decfloat approx, ans, logx
    approx = x
    factor = 4096
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fplogTaylor(logx, approx, digits)
    Call fpmul_si(ans, logx, factor, digits)
    result = ans
End Sub
 
Sub fpexp (result As decfloat, x As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    'taylor series
    Dim As decfloat fac, x9, temp, accum, p, term
    Dim As Long i, c
    Call si2fp(temp, 0, digits)
 
    fac.sign = 0
    fac.exponent = (BIAS + 1)
    Call set(fac.mantissa, 0, 10000000)
    For i = 1 To NUM_DWORDS
        Call set(fac.mantissa, i, 0)
    Next
    If fpcmp(x, temp, digits) = 0 Then result = fac: Exit Sub
    c = 1
    Call fpdiv_si(x9, x, 8192, digits) 'fpdiv_si(x, 67108864) '
    p = x9
    Call fpadd(accum, fac, x9, digits) '1 + x
 
    Do
        c = c + 1
        temp = accum
        Call fpdiv_si(fac, fac, c, digits)
        Call fpmul(p, p, x9, digits)
        Call fpmul(term, p, fac, digits)
        Call fpadd(accum, temp, term, digits)
    Loop Until fpcmp(accum, temp, digits) = 0
    For i = 1 To 13
        Call fpmul(accum, accum, accum, digits)
    Next
    result = accum
End Sub
 
Sub fpexp_cf (result As decfloat, x As decfloat, n As Long, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    '                       x
    '    cf =   ------------------------------
    '                         x^2
    '           2 + -------------------------
    '                           x^2
    '              6 + ----------------------
    '                             x^2
    '                 10 + ------------------
    '                               x^2
    '                     14 + --------------
    '                                 x^2
    '                          18 + ---------
    '                                   x^2
    '                              22 + -----
    '                                   26 + ...
    '           (1 + cf)
    ' exp(x) = ----------
    '           (1 - cf)
 
    Dim i As Long
    Dim As decfloat s, x1, x2, tmp, tmp2
    Call fpdiv_si(x1, x, 8192, digits) ' 2^13
    Call fpmul(x2, x1, x1, digits)
    Call si2fp(s, 4 * n + 6, digits)
    For i = n To 0 Step -1
        Call si2fp(tmp, 4 * i + 2, digits)
        Call fpdiv(s, x2, s, digits)
        Call fpadd(s, s, tmp, digits)
    Next
    Call fpdiv(s, x1, s, digits)
    Call si2fp(tmp, 1, digits)
    tmp2 = tmp
    Call fpadd(tmp, tmp, s, digits)
    Call fpsub(tmp2, tmp2, s, digits)
    Call fpdiv(s, tmp, tmp2, digits)
    For i = 1 To 13 ' s^13
        Call fpmul(s, s, s, digits)
    Next
    result = s
End Sub
 
Sub fplog_r (result As decfloat, x As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    '====================Log Guard==================
    Dim As decfloat x2, y, tmp
    Dim As Long i
    Dim As Long ex, terms, prec, n
    Dim As Double t, t2
 
    tmp.sign = 0
    tmp.exponent = (BIAS + 1)
    Call set(tmp.mantissa, 0, 10000000)
    For i = 1 To NUM_DWORDS
        Call set(tmp.mantissa, i, 0)
    Next
    If fpcmp(x, tmp, digits) = 0 Then result = x2: Exit Sub
    ex = (x.exponent And &H7FFFFFFF) - BIAS - 1
    t = git(x.mantissa, 0) + git(x.mantissa, 1) / 100000000
    t = t / 10000000
    t2 = Log(t) + Log(10) * ex
 
    n = Log((digits + 1) * 0.5) * 1.5
 
    Call str2fp(x2, Str$(t2))
    prec = 3
    terms = 2
    prec = 2 * prec - 1
    Call fpexp_cf(y, x2, terms, prec)
    Call fpsub(tmp, y, x, prec)
    Call fpdiv(tmp, tmp, y, prec)
    Call fpsub(x2, x2, tmp, prec)
    If digits > 4 And digits <= 9 Then
        prec = 2 * prec - 1
        terms = 2 * terms + 1
        Call fpexp_cf(y, x2, terms, prec)
        Call fpsub(tmp, y, x, prec)
        Call fpdiv(tmp, tmp, y, prec)
        Call fpsub(x2, x2, tmp, prec)
    ElseIf digits > 9 And digits < 18 Then
        For i = 1 To 3
            prec = 2 * prec - 1
            terms = 2 * terms + 1
            Call fpexp_cf(y, x2, terms, prec)
            Call fpsub(tmp, y, x, prec)
            Call fpdiv(tmp, tmp, y, prec)
            Call fpsub(x2, x2, tmp, prec)
        Next
    ElseIf digits > 17 And digits < 34 Then
        For i = 1 To 4
            prec = 2 * prec - 1
            terms = 2 * terms + 1
            Call fpexp_cf(y, x2, terms, prec)
            Call fpsub(tmp, y, x, prec)
            Call fpdiv(tmp, tmp, y, prec)
            Call fpsub(x2, x2, tmp, prec)
        Next
    ElseIf digits > 33 And digits < 65 Then
        For i = 1 To 5
            prec = 2 * prec - 1
            terms = 2 * terms + 1
            Call fpexp_cf(y, x2, terms, prec)
            Call fpsub(tmp, y, x, prec)
            Call fpdiv(tmp, tmp, y, prec)
            Call fpsub(x2, x2, tmp, prec)
        Next
    ElseIf digits > 64 Then
        For i = 1 To 6
            prec = 2 * prec - 1
            terms = 2 * terms + 1
            Call fpexp_cf(y, x2, terms, prec)
            Call fpsub(tmp, y, x, prec)
            Call fpdiv(tmp, tmp, y, prec)
            Call fpsub(x2, x2, tmp, prec)
        Next
    End If
    result = x2
End Sub
 
Sub fppow (result As decfloat, lhs As decfloat, rhs As decfloat)
    Dim As decfloat lhs2
    Call fplog(lhs2, lhs, NUM_DWORDS)
    Call fpmul(lhs2, lhs2, rhs, NUM_DWORDS)
    Call fpexp(lhs2, lhs2, NUM_DWORDS)
    result = lhs2
End Sub
 
Sub fpnroot (result As decfloat, x As decfloat, p As Long, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat ry, tmp, tmp2
    Dim As Double t, t2
    Dim As Long i, ex, l, prec
 
    l = Log((NUM_DIGITS + 8) * 0.0625) * 1.5
    ex = (x.exponent And &H7FFFFFFF) - BIAS - 1
    t = git(x.mantissa, 0) + git(x.mantissa, 1) / 100000000
    t = t / 10000000
    t2 = Log(t) / p + Log(10) * ex / p
    t2 = Exp(t2)
    Call str2fp(ry, Str$(t2))
    prec = 3
 
    Call fpipow(tmp, ry, p - 1, prec)
    Call fpdiv(tmp, x, tmp, prec)
    Call fpmul_si(tmp2, ry, p - 1, prec)
    Call fpadd(tmp2, tmp2, tmp, prec)
    Call fpdiv_si(ry, tmp2, p, prec)
    For i = 1 To l
        prec = 2 * prec - 1
        Call fpipow(tmp, ry, p - 1, prec)
        Call fpdiv(tmp, x, tmp, prec)
        Call fpmul_si(tmp2, ry, p - 1, prec)
        Call fpadd(tmp2, tmp2, tmp, prec)
        Call fpdiv_si(ry, tmp2, p, prec)
    Next
    result = ry
End Sub
 
Sub pi_brent_salamin (pi_bs As decfloat, digits_in As Unsigned Long)
    Dim As Unsigned Long digits
    digits = digits_in
    If digits > NUMBER_OF_DIGITS Then digits = NUMBER_OF_DIGITS
 
    Dim As Long limit2
    Dim As decfloat c0, c1, c2, c05
    Dim As decfloat a, b, sum
    Dim As decfloat ak, bk, ck
    Dim As decfloat ab, asq
    Dim As decfloat pow2, tmp
 
    limit2 = -digits + BIAS + 1
    Call si2fp(c0, 0, NUMBER_OF_DIGITS): ak = c0: bk = c0: ab = c0: asq = c0
    Call si2fp(c1, 1, NUMBER_OF_DIGITS): a = c1: ck = c1: pow2 = c1
    Call si2fp(c2, 2, NUMBER_OF_DIGITS): b = c2
    Call str2fp(c05, ".5"): sum = c05
    Call si2fp(pi_bs, 3, NUMBER_OF_DIGITS)
 
    Call fpsqr(b, b, NUMBER_OF_DIGITS)
    Call fpdiv(b, c1, b, NUMBER_OF_DIGITS)
    While fpcmp(ck, c0, NUMBER_OF_DIGITS) <> 0 And ck.exponent > limit2
        Call fpadd(ak, a, b, NUMBER_OF_DIGITS)
        Call fpmul(ak, c05, ak, NUMBER_OF_DIGITS)
        Call fpmul(ab, a, b, NUMBER_OF_DIGITS)
        Call fpsqr(bk, ab, NUMBER_OF_DIGITS)
        Call fpmul(asq, ak, ak, NUMBER_OF_DIGITS)
        Call fpsub(ck, asq, ab, NUMBER_OF_DIGITS)
        Call fpmul(pow2, pow2, c2, NUMBER_OF_DIGITS)
        Call fpmul(tmp, pow2, ck, NUMBER_OF_DIGITS)
        Call fpsub(sum, sum, tmp, NUMBER_OF_DIGITS)
        a = ak: b = bk
    Wend
    Call fpdiv(tmp, asq, sum, NUMBER_OF_DIGITS)
    Call fpmul(pi_bs, c2, tmp, NUMBER_OF_DIGITS)
End Sub
 
Sub fpsin (result As decfloat, x As decfloat, digits_in As Unsigned Long)
    Dim As decfloat XX, Term, Accum, p, temp2, fac, x_2
    Dim As decfloat pi2, circ, Ab
    Dim As Long precision
    precision = digits_in
 
    x_2 = x
    pi2 = pi_dec
    Call fpmul_si(circ, pi2, 2, precision)
    Call fpabs(Ab, x)
    If fpcmp(Ab, circ, precision) > 0 Then
        '======== CENTRALIZE ==============
        'floor/ceil to centralize
        Dim As decfloat tmp, tmp2
        If precision > 20 Then
            pi2 = pi_dec
        End If
        Call fpmul_si(pi2, pi2, 2, precision) 'got 2*pi
        Call fpdiv(tmp, x_2, pi2, precision)
        tmp2 = tmp
        Call fpfix(tmp, tmp) 'int part
        Call fpsub(tmp, tmp2, tmp, precision) 'frac part
        Call fpmul(tmp, tmp, pi2, precision)
        x_2 = tmp
    End If
 
    Dim As Long lm, limit2, i
    Dim As decfloat factor
    lm = precision
    limit2 = Int(-0.45344993886092585968 + 0.022333002852398072433 * lm + 5.0461814408333079844E-7 * lm * lm - 4.2338453039804235772E-11 * lm * lm * lm)
 
    Call si2fp(factor, 5, precision)
    Call fpipow(factor, factor, limit2, precision)
    Call fpdiv(x_2, x_2, factor, precision) 'x_=x_/5^limit2
    '==================================
    Dim As Long sign(3): sign(3) = 1
    Dim As Long c: c = 1
 
    Accum = x_2
    Call si2fp(fac, 1, precision)
    p = x_2
    Call fpmul(XX, x_2, x_2, precision)
    Do
        c = c + 2
        temp2 = Accum
        Call fpmul_si(fac, fac, c * (c - 1), precision)
        Call fpmul(p, p, XX, precision)
        Call fpdiv(Term, p, fac, precision)
        If sign(c And 3) Then
            Call fpsub(Accum, temp2, Term, precision)
        Else
            Call fpadd(Accum, temp2, Term, precision)
        End If
    Loop Until fpcmp(Accum, temp2, precision) = 0
    'multiply the result by 5^limit2
 
    For i = 1 To limit2
        Call fpmul(p, Accum, Accum, precision)
        Call fpmul(temp2, Accum, p, precision)
        '*** sin(5*x) = 5 * sin(x) - 20 * sin(x)^3 + 16 * sin(x)^5
        Call fpmul_si(Accum, Accum, 5, precision)
        Call fpmul_si(Term, temp2, 20, precision)
        Call fpmul_si(XX, temp2, 16, precision)
        Call fpmul(XX, XX, p, precision)
        Call fpsub(Accum, Accum, Term, precision)
        Call fpadd(Accum, Accum, XX, precision)
    Next i
    result = Accum
End Sub
 
Sub fpcos (result As decfloat, z As decfloat, digits_in As Unsigned Long)
    Dim As decfloat x_2, pi2
    Dim As Long precision
    precision = digits_in
    Call fpdiv_si(pi2, pi_dec, 2, precision)
    Call fpsub(x_2, pi2, z, precision)
 
    Call fpsin(result, x_2, precision)
End Sub
 
Sub fptan (result As decfloat, z As decfloat, digits_in As Unsigned Long)
    Dim As decfloat x_2, s, c
    Dim As Long precision
    precision = digits_in
    x_2 = z
    Call fpsin(s, x_2, precision)
    x_2 = z
    Call fpcos(c, x_2, precision)
    Call fpdiv(result, s, c, precision)
End Sub
 
Sub fpatn (result As decfloat, x As decfloat, digits_in As Unsigned Long)
    Dim As Long precision, z
    precision = digits_in
    Dim As Long sign(3): sign(3) = 1
    Dim As Unsigned Long c: c = 1
    Dim As decfloat XX, Term, Accum, strC, x_2, mt, mt2, p
    Dim As decfloat decnum, one, decnum2, factor
 
    decnum2 = x
    decnum2.sign = 0
    Call si2fp(one, 1, precision)
    If fpcmp(decnum2, one, precision) = 0 Then
        Call fpdiv_si(result, pi_dec, 4, precision)
        result.sign = x.sign
        Exit Sub
    End If
    decnum2.sign = x.sign
    Dim As Long limit2: limit2 = 16
    Call si2fp(factor, SHL(2, limit2 - 1), precision)
    For z = 1 To limit2
        Call fpmul(decnum, decnum2, decnum2, precision)
        Call fpadd(decnum, decnum, one, precision)
        Call fpsqr(decnum, decnum, precision)
        Call fpadd(decnum, decnum, one, precision)
        Call fpdiv(decnum, decnum2, decnum, precision)
        decnum2 = decnum
    Next z
 
    mt = decnum
    x_2 = decnum
    p = decnum
    Call fpmul(XX, x_2, x_2, precision)
    Do
        c = c + 2
        mt2 = mt
        Call si2fp(strC, c, precision)
        Call fpmul(p, p, XX, precision)
        Call fpdiv(Term, p, strC, precision)
        If sign(c And 3) Then
            Call fpsub(Accum, mt, Term, precision)
        Else
            Call fpadd(Accum, mt, Term, precision)
        End If
        Swap mt, Accum
    Loop Until fpcmp(mt, mt2, precision) = 0
    Call fpmul(result, factor, mt, precision)
End Sub
 
Sub fpasin (result As decfloat, x As decfloat, digits_in As Unsigned Long)
    Dim As Long precision
    precision = digits_in
    Dim As Double num
    ' ARCSIN = ATN(x / SQR(-x * x + 1))
    '============= ARCSIN GUARD =========
    num = fp2dbl(x)
    If num > 1 Then Exit Sub
    If num < -1 Then Exit Sub
    '========================
    Dim As decfloat one, T, B, term1, minusone
    Call si2fp(one, 1, precision)
    Call si2fp(minusone, -1, precision)
    T = x
    Call fpmul(B, x, x, precision) 'x*x
    'for 1 and -1
    If fpcmp(B, one, precision) = 0 Then
        Dim As decfloat two, atn1
        Call si2fp(two, 2, precision)
        Call fpatn(atn1, one, precision)
        If fpcmp(x, minusone, precision) = 0 Then
            Call fpmul(two, two, atn1, precision)
            Call fpmul(result, two, minusone, precision)
            Exit Sub
        Else
            Call fpmul(result, two, atn1, precision)
            Exit Sub
        End If
    End If
    Call fpsub(B, one, B, precision) '1-x*x
    Call fpsqr(B, B, precision) 'sqr(1-x*x)
    Call fpdiv(term1, T, B, precision)
    Call fpatn(result, term1, precision)
End Sub
 
Sub fpacos (result As decfloat, x As decfloat, digits_in As Unsigned Long)
    Dim As Long precision
    precision = digits_in
    Dim As decfloat one, minusone, two, atn1, tail, T, B, term1, atnterm1 ',_x,temp
    Dim As Double num
    'ARCCOS = ATN(-x / SQR(-x * x + 1)) + 2 * ATN(1)
    '============= ARCCOS GUARD =========
    num = fp2dbl(x)
    If num > 1 Then Exit Sub
    If num < -1 Then Exit Sub
    '========================
    Call si2fp(one, 1, precision)
    Call si2fp(minusone, -1, precision)
    Call si2fp(two, 2, precision)
    Call fpatn(atn1, one, precision)
    Call fpmul(tail, two, atn1, precision) '2*atn(1)
    Call fpmul(T, minusone, x, precision) '-x
    Call fpmul(B, x, x, precision) 'x*x
    If fpcmp(B, one, precision) = 0 Then
        'for 1 and -1
        If fpcmp(x, minusone, precision) = 0 Then
            Dim As decfloat four
            Call si2fp(four, 4, precision)
            Call fpmul(result, four, atn1, precision)
            Exit Sub
        Else
            Call si2fp(result, 0, precision)
            Exit Sub
        End If
    End If
    Call fpsub(B, one, B, precision) '1-x*x
    Call fpsqr(B, B, precision) 'sqr(1-x*x)
    Call fpdiv(term1, T, B, precision)
    Call fpatn(atnterm1, term1, precision)
    Call fpadd(result, atnterm1, tail, precision)
End Sub
 

output

Code: [Select]

factorial(10000) to  512  digits
 2.8462596809170545189064132121198688901480514017027992307941799942744113400037644437729907867577847758158840621423175288300423399401535187390524211613827161748198241998275924182892597878981242531205946599625986706560161572036032397926328736717055741975962099479720346153698119897092611277500484198845410475544642442136573303076703628825803548967461117097369578603670191071512730587281041158640561281165385325968425825995584688146430425589836649317059251717204276597407446133400054194052462303436869154059404066227E+35659

in  .494140625  seconds

3/9 to  512  digits
 0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333

in  0  seconds

and square it
 0.11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111

ln(2) to  512  digits
 0.69314718055994530941723212145817656807550013436025525412068000949339362196969471560586332699641868754200148102057068573368552023575813055703267075163507596193072757082837143519030703862389167347112335011536449795523912047517268157493206515552473413952588295045300709532636664265410423915781495204374043038550080194417064167151864471283996817178454695702627163106454615025720740248163777338963855069526066834113727387372292895649354702576265209885969320196505855476470330679365443254763274495125040606943814710468

in  .111328125  seconds

and the inverse
 1.9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999

in  .21875  seconds

reciprocal of 3
 0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333

in  0  seconds

Press any key to continue

bplus · « **Reply #1 on:** October 10, 2021, 09:19:30 am »

Might be slower but I bet it runs rings around my string math. Looks like you might be using binary math?

Very impressive set of functions Power and Log... (I'm just scanning code).

I went to run the code and the Option Explicit line is red-lined in IDE, not first statement in code. Was that a last minute addition? Seems fixed when I moved it just under $NoPrefix which also demand to be top line. ;-)) Certainly $noprefix wins top spot specially if want "Option Explicit" without _prefix.

Ha! didn't take but a second or 2! Fast enough for me! Very nice work @jack

FellippeHeitor · « **Reply #2 on:** October 10, 2021, 09:31:45 am »

Option _Explicit, _ExplicitArray, and $NoPrefix no longer have to be the first line in a program in v2.0 😉😉

jack · « **Reply #3 on:** October 10, 2021, 09:33:20 am »

thanks bplus for testing
I am using base 100000000, that is 8 digits per element

jack · « **Reply #4 on:** October 10, 2021, 10:04:13 am »

bplus
the reason I added a precision parameter was so that you could vary the precision in Newton–Raphson type calculations, I use it in the square-root and reciprocal routines
it's a pity that QB64 doesn't support array in type and optional parameters

bplus · « **Reply #5 on:** October 10, 2021, 10:45:01 am »

Quote

QB64 doesn't support array in type and optional parameters.

Actually might be an advantage if you are really, really after speed and not so much pretty object like coding. It will discourage you from using Type. What was I working on that I was trying to get done much much faster? Oh yeah, the thing with word mods, trying to get from one word to another by only changing a letter at a time. This also worked great for me in tracking worm segments (x, y) positions, DIM x(1 to nWorms, 1 to MaxSegmentsAllowed), y(1 to nWorms, 1 to MaxSegmentsAllowed), ItemLimit(1 to nWorms)

Anyway the lesson was to dump Type stuff and use @2D dimension arrays (1 to nItems, 1 to MaxNumberYouWillEverNeedPlusOneForSafety). You will see significant jump in speed performance, no strings, no ReDim _Preserves, just track highest element for each item but I suspect you will use same amount for each "String Array" part probably the number of digits you want your precision plus a safety margin for rounding errors that creep in after thousands of iterations.

jack · « **Reply #6 on:** October 10, 2021, 11:11:01 am »

bplus
I find that type's are essential and operator and function overloading are very handy and the slowdown is minimum, my guess less than 5% slowdown but you get huge benefits.

bplus · « **Reply #7 on:** October 10, 2021, 11:50:38 am »

Quote from: jack on October 10, 2021, 11:11:01 am

bplus
I find that type's are essential and operator and function overloading are very handy and the slowdown is minimum, my guess less than 5% slowdown but you get huge benefits.

jack,

I'd sure like to pick your brain on what is "essential" about Types, this thread probably not appropriate for such discussion and besides all you might mean by that is they offer code organization and clarity, to which I would agree whole heartedly.

But to me essential means you couldn't do something unless in Type form, I am pretty sure in QB64 that is not the case, that's why I want to pick your brain a bit more. I have learned to watch out when I am pretty sure about something ;-))

PS I quoted you in case you decided to change out " essential" in edit. I seem to be picking arguments with people's wording today, it is one way to learn not necessarily essential to building friendships. :)

jack · « **Reply #8 on:** October 10, 2021, 12:09:08 pm »

where are the emoticons?
sure could use laughing emoticon now, but you basically stated it
just to be a little absurd, imagine we didn't have long's or double's just loose bytes all over the place, what a nightmare it would be to write readable and maintainable code

SMcNeill · « **Reply #9 on:** October 10, 2021, 12:11:03 pm »

Quote from: jack on October 10, 2021, 12:09:08 pm

where are the emoticons?
sure could use laughing emoticon now, but you basically stated it
just to be a little absurd, imagine we didn't have long's or double's just loose bytes all over the place, what a nightmare it would be to write readable and maintainable code

Sounds like coding everything using just _MEMNEW to hold the data.

And we has emoticons! 😂🤣🥲😅😇🙃😊☺️🙂😘😗😙🥰😍😋😚😛🤪😜😝🤨🧐🤓😗😁😆.... And MANY more!

bplus · « **Reply #10 on:** October 10, 2021, 12:58:11 pm »

Quote

And we has emoticons! 😂🤣🥲😅😇🙃😊☺️🙂😘😗😙🥰😍😋😚😛🤪😜😝🤨🧐🤓😗😁😆.... And MANY more!

Must be an QB64 Developer privilege, closet I am finding is this font.

Wait!

Quote

Emoticons are punctuation marks, letters, and numbers used to create pictorial icons that generally display an emotion or sentiment.

Quote

Emoji were invented in 1999 by Shigetaka Kurita and were intended for a Japanese user base. The first emoji were very simple—only 12 pixels by 12 pixels—and were inspired by manga art and kanji characters. In order to attract Japanese customers, Apple hid an emoji keyboard in the first iPhone back in 2007, but North American users quickly became aware of the keyboard. Now, emoji are available in almost all messaging apps, and while different apps have distinct emoji styles, emoji can translate across platforms, thanks to Unicode. This is why an iPhone user is able to receive the smiling pile of poo emoji from someone using a Samsung Galaxy.

So, if you come across a smiley face that contains a character you can find on your computer keyboard, it’s an emoticon. If it’s a little cartoon figure that is free from the binds of punctuation, numbers, and letters, it’s an emoji.

Emoji must be one of those things know to QB64 Developers.

FellippeHeitor · « **Reply #11 on:** October 10, 2021, 12:59:47 pm »

Lol, nothing exclusive. Emoji are available on all mobile devices and all modern operating systems have a way to allow you to insert them (virtual keyboards, look for those).
🤯⌨️

jack · « **Reply #12 on:** October 10, 2021, 01:59:02 pm »

https://fsymbols.com/emoticons/

Cobalt · « **Reply #13 on:** October 12, 2021, 03:51:12 am »

Quote from: FellippeHeitor on October 10, 2021, 12:59:47 pm

Lol, nothing exclusive. Emoji are available on all mobile devices and all modern operating systems have a way to allow you to insert them (virtual keyboards, look for those).

Woo Hoo my Vista machine is "Modern" again!

🥳

Long standing advocate for Arrays in UDTs.

jack · « **Reply #14 on:** October 12, 2021, 10:24:52 pm »

if someone were to actually use this it would be better to have a bi, bm
decfloat.bi

Code: QB64: [Select]

'Const NUMBER_OF_DIGITS = 512
Const NUM_DIGITS = NUMBER_OF_DIGITS
Const NUM_DWORDS = NUM_DIGITS \ 8
Const BIAS = 1073741824 '2 ^ 30
Const MANTISSA_BYTES = (NUM_DWORDS + 1) * 8 * 4
 
Type decfloat
    As Long sign
    As Unsigned Long exponent
    As String * Mantissa_bytes mantissa
End Type
 
' Error definitions
 
Const DIVZ_ERR = 1 'Divide by zero
Const EXPO_ERR = 2 'Exponent overflow error
Const EXPU_ERR = 3 'Exponent underflow error
Dim Shared As decfloat pi_dec
Call pi_brent_salamin(pi_dec, NUMBER_OF_DIGITS)
 

decfloat.bm

Code: QB64: [Select]

Sub set (s$, i&, v~&)
    Mid$(s$, 4 * i& + 1, 4) = MKL$(v~&)
End Sub
 
Function git~& (s$, i&)
    git = CVL(Mid$(s$, 4 * i& + 1, 4))
End Function
 
Function fp2str_exp$ (n As decfloat, places As Long)
    Dim As Long i, ex
    Dim As String v, f, ts
    If n.exponent > 0 Then
        ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
    Else
        ex = 0
    End If
    If n.sign Then v = "-" Else v = " "
    ts = Str$(git(n.mantissa, 0))
    ts = Trim$(ts)
    If Len(ts) < 8 Then
        ts = ts + String$(8 - Len(ts), "0")
    End If
    v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
    For i = 1 To NUM_DWORDS - 1
        ts = Str$(git(n.mantissa, i))
        ts = Trim$(ts)
        If Len(ts) < 8 Then
            ts = String$(8 - Len(ts), "0") + ts
        End If
        v = v + ts
    Next
    v = Left$(v, places + 3)
    f = Trim$(Str$(Abs(ex)))
    f = String$(5 - Len(f), "0") + f
    If ex < 0 Then v = v + "E-" Else v = v + "E+"
    v = v + f
    fp2str_exp$ = v
End Function
 
'long part of num
Sub fpfix (result As decfloat, num As decfloat)
    Dim As decfloat ip
    Dim As Long ex, ex2, j, k
 
    ex = (num.exponent And &H7FFFFFFF) - BIAS
    If ex < 1 Then
        result = ip: Exit Sub
    End If
    If ex >= (NUM_DIGITS) Then
        result = num: Exit Sub
    End If
    ex2 = ex \ 8
    k = ex2
    j = ex Mod 8
    While ex2 > 0
        ex2 = ex2 - 1
        Call set(ip.mantissa, ex2, git(num.mantissa, ex2))
    Wend
    If j = 1 Then
        Call set(ip.mantissa, k, 10000000 * (git(num.mantissa, k) \ 10000000))
    ElseIf j = 2 Then
        Call set(ip.mantissa, k, 1000000 * (git(num.mantissa, k) \ 1000000))
    ElseIf j = 3 Then
        Call set(ip.mantissa, k, 100000 * (git(num.mantissa, k) \ 100000))
    ElseIf j = 4 Then
        Call set(ip.mantissa, k, 10000 * (git(num.mantissa, k) \ 10000))
    ElseIf j = 5 Then
        Call set(ip.mantissa, k, 1000 * (git(num.mantissa, k) \ 1000))
    ElseIf j = 6 Then
        Call set(ip.mantissa, k, 100 * (git(num.mantissa, k) \ 100))
    ElseIf j = 7 Then
        Call set(ip.mantissa, k, 10 * (git(num.mantissa, k) \ 10))
    ElseIf j = 8 Then
        Call set(ip.mantissa, k, git(num.mantissa, k))
    End If
    ip.exponent = ex + BIAS
    ip.sign = num.sign
    result = ip
End Sub
 
Function fp2dbl# (n As decfloat)
    Dim As Long ex
    Dim As String v, f, ts
    If n.exponent > 0 Then
        ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
    Else
        ex = 0
    End If
    If n.sign Then v = "-" Else v = " "
    ts = Trim$(Str$(git(n.mantissa, 0)))
    If Len(ts) < 8 Then
        ts = ts + String$(8 - Len(ts), "0")
    End If
    v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
    ts = Trim$(Str$(git(n.mantissa, 1)))
    If Len(ts) < 8 Then
        ts = String$(8 - Len(ts), "0") + ts
    End If
    v = v + ts
 
    f = Str$(Abs(ex))
    f = String$(5 - Len(f), "0") + f
    If ex < 0 Then v = v + "E-" Else v = v + "E+"
    v = v + f
    fp2dbl# = Val(v)
End Function
 
Function fp2str_fix$ (n As decfloat, places As Long)
    Dim As Long i, ex
    Dim As String v, ts, s
    If n.exponent > 0 Then
        ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
    Else
        ex = 0
    End If
    If n.sign Then s = "-" Else s = " "
    ts = Trim$(Str$(git(n.mantissa, 0)))
    If Len(ts) < 8 Then
        ts = ts + String$(8 - Len(ts), "0")
    End If
    v = ts
    For i = 1 To NUM_DWORDS - 1
        ts = Trim$(Str$(git(n.mantissa, i)))
        If Len(ts) < 8 Then
            ts = String$(8 - Len(ts), "0") + ts
        End If
        v = v + ts
    Next
    If places < NUM_DIGITS Then
        v = Left$(v, places)
    End If
    If ex = 0 Then
        v = Left$(v, 1) + "." + Mid$(v, 2)
    ElseIf ex < 0 Then
        v = "0." + String$(Abs(ex) - 1, "0") + v
    ElseIf ex > 0 Then
        v = Left$(v, ex + 1) + "." + Mid$(v, ex + 2)
    End If
    fp2str_fix$ = s + v
End Function
 
Function fp2str$ (n As decfloat, places As Long)
    Dim As Long ex
    If n.exponent > 0 Then
        ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
    Else
        ex = 0
    End If
    If Abs(ex) < places Then
        fp2str = fp2str_fix(n, places)
    Else
        fp2str = fp2str_exp(n, places)
    End If
End Function
 
Sub str2fp (result As decfloat, value As String)
    Dim As Long j, s, d, e, ep, ex, es, i, f, fp, fln
    Dim As String c, f1, f2, f3, ts
    Dim As Unsigned Long ulng
    Dim n As decfloat
    j = 1
    s = 1
    d = 0
    e = 0
    ep = 0
    ex = 0
    es = 1
    i = 0
    f = 0
    fp = 0
    f1 = ""
    f2 = ""
    f3 = ""
    value = UCase$(value)
    fln = Len(value)
 
    While j <= fln
        c = Mid$(value, j, 1)
        If ep = 1 Then
            If c = " " Then
                j = j + 1
                GoTo skip_while
            End If
            If c = "-" Then
                es = -es
                c = ""
            End If
            If c = "+" Then
                j = j + 1
                GoTo skip_while
            End If
            If (c = "0") And (f3 = "") Then
                j = j + 1
                GoTo skip_while
            End If
            If (c > "/") And (c < ":") Then 'c is digit between 0 and 9
                f3 = f3 + c
                ex = 10 * ex + (Asc(c) - 48)
                j = j + 1
                GoTo skip_while
            End If
        End If
 
        If c = " " Then
            j = j + 1
            GoTo skip_while
        End If
        If c = "-" Then
            s = -s
            j = j + 1
            GoTo skip_while
        End If
        If c = "+" Then
            j = j + 1
            GoTo skip_while
        End If
        If c = "." Then
            If d = 1 Then
                j = j + 1
                GoTo skip_while
            End If
            d = 1
        End If
        If (c > "/") And (c < ":") Then 'c is digit between 0 and 9
            If ((c = "0") And (i = 0)) Then
                If d = 0 Then
                    j = j + 1
                    GoTo skip_while
                End If
                If (d = 1) And (f = 0) Then
                    e = e - 1
                    j = j + 1
                    GoTo skip_while
                End If
            End If
            If d = 0 Then
                f1 = f1 + c
                i = i + 1
            Else
                If (c > "0") Then
                    fp = 1
                End If
                f2 = f2 + c
                f = f + 1
            End If
        End If
        If c = "E" Or c = "D" Then
            ep = 1
        End If
        j = j + 1
        skip_while:
    Wend
    If fp = 0 Then
        f = 0
        f2 = ""
    End If
 
    If s = -1 Then s = &H8000 Else s = 0
    n.sign = s
    ex = es * ex - 1 + i + e
    f1 = f1 + f2
    f1 = Mid$(f1, 1, 1) + Right$(f1, Len(f1) - 1)
    fln = Len(f1)
    If Len(f1) > (NUM_DIGITS + 1 + 8) Then
        f1 = Mid$(f1, 1, (NUM_DIGITS + 1 + 8))
    End If
    While Len(f1) < (NUM_DIGITS + 1 + 8)
        f1 = f1 + "0"
    Wend
    j = 1
    For i = 0 To NUM_DWORDS
        ts = Mid$(f1, j, 8)
        ulng = Val(ts)
        Call set(n.mantissa, i, ulng)
        If ulng <> 0 Then fp = 1
        j = j + 8
    Next
    If fp Then n.exponent = (ex + BIAS + 1) Else n.exponent = 0
    result = n
End Sub
 
Sub si2fp (result As decfloat, m As Integer64, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat fac1
    Dim As Long i
    Dim As Integer64 n
    n = Abs(m)
 
    If n > 9999999999999999 Then
        Call str2fp(fac1, Str$(m))
        result = fac1: Exit Sub
    End If
 
    For i = 1 To digits
        Call set(fac1.mantissa, i, 0)
    Next
 
    If m = 0 Then
        fac1.exponent = 0
        fac1.sign = 0
        result = fac1: Exit Sub
    End If
 
    fac1.exponent = BIAS
    If n < 100000000 Then
        If n < 10 Then
            Call set(fac1.mantissa, 0, n * 10000000)
            fac1.exponent = fac1.exponent + 1
        ElseIf n < 100 Then
            Call set(fac1.mantissa, 0, n * 1000000)
            fac1.exponent = fac1.exponent + 2
        ElseIf n < 1000 Then
            Call set(fac1.mantissa, 0, n * 100000)
            fac1.exponent = fac1.exponent + 3
        ElseIf n < 10000 Then
            Call set(fac1.mantissa, 0, n * 10000)
            fac1.exponent = fac1.exponent + 4
        ElseIf n < 100000 Then
            Call set(fac1.mantissa, 0, n * 1000)
            fac1.exponent = fac1.exponent + 5
        ElseIf n < 1000000 Then
            Call set(fac1.mantissa, 0, n * 100)
            fac1.exponent = fac1.exponent + 6
        ElseIf n < 10000000 Then
            Call set(fac1.mantissa, 0, n * 10)
            fac1.exponent = fac1.exponent + 7
        ElseIf n < 100000000 Then
            Call set(fac1.mantissa, 0, n)
            fac1.exponent = fac1.exponent + 8
        End If
    End If
    If n > 99999999 Then
        fac1.exponent = fac1.exponent + 8
        If n < 1000000000 Then
            Call set(fac1.mantissa, 0, n \ 10)
            Call set(fac1.mantissa, 1, (n Mod 10) * 10000000)
            fac1.exponent = fac1.exponent + 1
        ElseIf n < 100000000000 Then
            Call set(fac1.mantissa, 0, n \ 100)
            Call set(fac1.mantissa, 1, (n Mod 100) * 1000000)
            fac1.exponent = fac1.exponent + 2
        ElseIf n < 1000000000000 Then
            Call set(fac1.mantissa, 0, n \ 1000)
            Call set(fac1.mantissa, 1, (n Mod 1000) * 100000)
            fac1.exponent = fac1.exponent + 3
        ElseIf n < 10000000000000 Then
            Call set(fac1.mantissa, 0, n \ 10000)
            Call set(fac1.mantissa, 1, (n Mod 10000) * 10000)
            fac1.exponent = fac1.exponent + 4
        ElseIf n < 100000000000000 Then
            Call set(fac1.mantissa, 0, n \ 100000)
            Call set(fac1.mantissa, 1, (n Mod 100000) * 1000)
            fac1.exponent = fac1.exponent + 5
        ElseIf n < 1000000000000000 Then
            Call set(fac1.mantissa, 0, n \ 1000000)
            Call set(fac1.mantissa, 1, (n Mod 1000000) * 100)
            fac1.exponent = fac1.exponent + 6
        ElseIf n < 10000000000000000 Then
            Call set(fac1.mantissa, 0, n \ 10000000)
            Call set(fac1.mantissa, 1, (n Mod 10000000) * 10)
            fac1.exponent = fac1.exponent + 7
        ElseIf n < 100000000000000000 Then
            Call set(fac1.mantissa, 0, n \ 100000000)
            Call set(fac1.mantissa, 1, n Mod 100000000)
            fac1.exponent = fac1.exponent + 8
        End If
    End If
    If m < 0 Then
        fac1.sign = &H8000
    Else
        fac1.sign = 0
    End If
    result = fac1
End Sub
 
Sub RSHIFT_1 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = digits To 1 Step -1
        v1 = git(mantissa.mantissa, i) \ 10
        v2 = git(mantissa.mantissa, i - 1) Mod 10
        v2 = v2 * 10000000 + v1
        Call set(mantissa.mantissa, i, v2)
    Next
    Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 10)
End Sub
 
Sub LSHIFT_1 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = 0 To digits - 1
        v1 = git(mantissa.mantissa, i) Mod 10000000
        v2 = git(mantissa.mantissa, i + 1) \ 10000000
        Call set(mantissa.mantissa, i, v1 * 10 + v2)
        Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 10000000)
    Next
    Call set(mantissa.mantissa, digits, 10 * (git(mantissa.mantissa, digits) Mod 10000000))
End Sub
 
Sub RSHIFT_2 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = digits To 1 Step -1
        v1 = git(mantissa.mantissa, i) \ 100
        v2 = git(mantissa.mantissa, i - 1) Mod 100
        v2 = v2 * 1000000 + v1
        Call set(mantissa.mantissa, i, v2)
    Next
    Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 100)
End Sub
 
Sub LSHIFT_2 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = 0 To digits - 1
        v1 = git(mantissa.mantissa, i) Mod 1000000
        v2 = git(mantissa.mantissa, i + 1) \ 1000000
        Call set(mantissa.mantissa, i, v1 * 100 + v2)
        Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 1000000)
    Next
    Call set(mantissa.mantissa, digits, 100 * (git(mantissa.mantissa, digits) Mod 1000000))
End Sub
 
Sub RSHIFT_3 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = digits To 1 Step -1
        v1 = git(mantissa.mantissa, i) \ 1000
        v2 = git(mantissa.mantissa, i - 1) Mod 1000
        v2 = v2 * 100000 + v1
        Call set(mantissa.mantissa, i, v2)
    Next
    Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 1000)
End Sub
 
Sub LSHIFT_3 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = 0 To digits - 1
        v1 = git(mantissa.mantissa, i) Mod 100000
        v2 = git(mantissa.mantissa, i + 1) \ 100000
        Call set(mantissa.mantissa, i, v1 * 1000 + v2)
        Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 100000)
    Next
    Call set(mantissa.mantissa, digits, 1000 * (git(mantissa.mantissa, digits) Mod 100000))
End Sub
 
Sub RSHIFT_4 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = digits To 1 Step -1
        v1 = git(mantissa.mantissa, i) \ 10000
        v2 = git(mantissa.mantissa, i - 1) Mod 10000
        v2 = v2 * 10000 + v1
        Call set(mantissa.mantissa, i, v2)
    Next
    Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 10000)
End Sub
 
Sub LSHIFT_4 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = 0 To digits - 1
        v1 = git(mantissa.mantissa, i) Mod 10000
        v2 = git(mantissa.mantissa, i + 1) \ 10000
        Call set(mantissa.mantissa, i, v1 * 10000 + v2)
        Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 10000)
    Next
    Call set(mantissa.mantissa, digits, 10000 * (git(mantissa.mantissa, digits) Mod 10000))
End Sub
 
Sub RSHIFT_5 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = digits To 1 Step -1
        v1 = git(mantissa.mantissa, i) \ 100000
        v2 = git(mantissa.mantissa, i - 1) Mod 100000
        v2 = v2 * 1000 + v1
        Call set(mantissa.mantissa, i, v2)
    Next
    Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 100000)
End Sub
 
Sub LSHIFT_5 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = 0 To digits - 1
        v1 = git(mantissa.mantissa, i) Mod 1000
        v2 = git(mantissa.mantissa, i + 1) \ 1000
        Call set(mantissa.mantissa, i, v1 * 100000 + v2)
        Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 1000)
    Next
    Call set(mantissa.mantissa, digits, 100000 * (git(mantissa.mantissa, digits) Mod 1000))
End Sub
 
Sub RSHIFT_6 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = digits To 1 Step -1
        v1 = git(mantissa.mantissa, i) \ 1000000
        v2 = git(mantissa.mantissa, i - 1) Mod 1000000
        v2 = v2 * 100 + v1
        Call set(mantissa.mantissa, i, v2)
    Next
    Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 1000000)
End Sub
 
Sub LSHIFT_6 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = 0 To digits - 1
        v1 = git(mantissa.mantissa, i) Mod 100
        v2 = git(mantissa.mantissa, i + 1) \ 100
        Call set(mantissa.mantissa, i, v1 * 1000000 + v2)
        Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 100)
    Next
    Call set(mantissa.mantissa, digits, 1000000 * (git(mantissa.mantissa, digits) Mod 100))
End Sub
 
Sub RSHIFT_7 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = digits To 1 Step -1
        v1 = git(mantissa.mantissa, i) \ 10000000
        v2 = git(mantissa.mantissa, i - 1) Mod 10000000
        v2 = v2 * 10 + v1
        Call set(mantissa.mantissa, i, v2)
    Next
    Call set(mantissa.mantissa, 0, git(mantissa.mantissa, 0) \ 10000000)
End Sub
 
Sub LSHIFT_7 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As Unsigned Long v1, v2
    Dim As Long i
    For i = 0 To digits - 1
        v1 = git(mantissa.mantissa, i) Mod 10
        v2 = git(mantissa.mantissa, i + 1) \ 10
        Call set(mantissa.mantissa, i, v1 * 10000000 + v2)
        Call set(mantissa.mantissa, i + 1, git(mantissa.mantissa, i + 1) Mod 10)
    Next
    Call set(mantissa.mantissa, digits, 10000000 * (git(mantissa.mantissa, digits) Mod 10))
End Sub
 
Sub RSHIFT_8 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As Long i
    For i = digits To 1 Step -1
        Call set(mantissa.mantissa, i, git(mantissa.mantissa, i - 1))
    Next
    Call set(mantissa.mantissa, 0, 0)
End Sub
 
Sub LSHIFT_8 (mantissa As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As Long i
    For i = 0 To digits - 1
        Call set(mantissa.mantissa, i, git(mantissa.mantissa, i + 1))
    Next
    Call set(mantissa.mantissa, digits, 0)
End Sub
 
Function fpcmp& (x As decfloat, y As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As Long c, i
    If x.sign < y.sign Then fpcmp = -1: Exit Function
    If x.sign > y.sign Then fpcmp = 1: Exit Function
    If x.sign = y.sign Then
        If x.exponent = y.exponent Then
            For i = 0 To digits
                c = git(x.mantissa, i) - git(y.mantissa, i)
                If c <> 0 Then Exit For
            Next
            If c < 0 Then fpcmp = -1: Exit Function
            If c = 0 Then fpcmp = 0: Exit Function
            If c > 0 Then fpcmp = 1: Exit Function
        End If
        If x.exponent < y.exponent Then fpcmp = -1: Exit Function
        If x.exponent > y.exponent Then fpcmp = 1: Exit Function
    End If
    ' if we reach this point it means that the signs are different
    ' and if the sign of x is set meaning that x is negative then x < y
 
End Function
 
Sub NORM_FAC1 (fac1 As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    ' normalize the number in fac1
    ' all routines exit through this one.
 
    'see if the mantissa is all zeros.
    'if so, set the exponent and sign equal to 0.
    Dim As Long i, er, f
    er = 0: f = 0
    For i = 0 To digits
        If git(fac1.mantissa, i) > 0 Then f = 1
    Next
    If f = 0 Then
        fac1.exponent = 0
        fac1.sign = 0
        Exit Sub
        'if the highmost Digit in fac1_man is nonzero,
        'shift the mantissa right 1 Digit and
        'increment the exponent
    ElseIf git(fac1.mantissa, 0) > 99999999 Then
        Call RSHIFT_1(fac1, digits)
        fac1.exponent = fac1.exponent + 1
    Else
        'now shift fac1_man 1 to the left until a
        'nonzero digit appears in the next-to-highest
        'Digit of fac1_man.  decrement exponent for
        'each shift.
        While git(fac1.mantissa, 0) = 0
            Call LSHIFT_8(fac1, digits)
            fac1.exponent = fac1.exponent - 8
            If fac1.exponent = 0 Then
                Print "NORM_FAC1=EXPU_ERR"
                Exit Sub
            End If
        Wend
        If git(fac1.mantissa, 0) < 10 Then
            Call LSHIFT_7(fac1, digits)
            fac1.exponent = fac1.exponent - 7
        ElseIf git(fac1.mantissa, 0) < 100 Then
            Call LSHIFT_6(fac1, digits)
            fac1.exponent = fac1.exponent - 6
        ElseIf git(fac1.mantissa, 0) < 1000 Then
            Call LSHIFT_5(fac1, digits)
            fac1.exponent = fac1.exponent - 5
        ElseIf git(fac1.mantissa, 0) < 10000 Then
            Call LSHIFT_4(fac1, digits)
            fac1.exponent = fac1.exponent - 4
        ElseIf git(fac1.mantissa, 0) < 100000 Then
            Call LSHIFT_3(fac1, digits)
            fac1.exponent = fac1.exponent - 3
        ElseIf git(fac1.mantissa, 0) < 1000000 Then
            Call LSHIFT_2(fac1, digits)
            fac1.exponent = fac1.exponent - 2
        ElseIf git(fac1.mantissa, 0) < 10000000 Then
            Call LSHIFT_1(fac1, digits)
            fac1.exponent = fac1.exponent - 1
        End If
    End If
    'check for overflow/underflow
    If fac1.exponent < 0 Then
        Print "NORM_FAC1=EXPO_ERR"
    End If
End Sub
 
Sub fpadd_aux (fac1 As decfloat, fac2 As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As Unsigned Long v, c, i
    c = 0
    For i = digits To 1 Step -1
        v = git(fac2.mantissa, i) + git(fac1.mantissa, i) + c
        If v > 99999999 Then
            v = v - 100000000
            c = 1
        Else
            c = 0
        End If
        Call set(fac1.mantissa, i, v)
    Next
    v = git(fac1.mantissa, 0) + git(fac2.mantissa, 0) + c
    Call set(fac1.mantissa, 0, v)
 
    Call NORM_FAC1(fac1, digits)
 
End Sub
 
Sub fpsub_aux (fac1 As decfloat, fac2 As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As Long v, c, i
    c = 0
    For i = digits To 1 Step -1
        v = git(fac1.mantissa, i) - git(fac2.mantissa, i) - c
        If v < 0 Then
            v = v + 100000000
            c = 1
        Else
            c = 0
        End If
        Call set(fac1.mantissa, i, v)
    Next
    v = git(fac1.mantissa, 0) - git(fac2.mantissa, 0) - c
    Call set(fac1.mantissa, 0, v)
 
    Call NORM_FAC1(fac1, digits)
End Sub
 
Sub fpadd (result As decfloat, x As decfloat, y As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    Dim As decfloat fac1, fac2
    Dim As Long i, t, c, xsign, ysign
 
    xsign = x.sign: x.sign = 0
    ysign = y.sign: y.sign = 0
    c = fpcmp(x, y, digits)
 
    x.sign = xsign
    y.sign = ysign
    If c < 0 Then
        fac1 = y
        fac2 = x
    Else
        fac1 = x
        fac2 = y
    End If
    t = fac1.exponent - fac2.exponent
    t = ((fac1.exponent And &H7FFFFFFF) - BIAS - 1) - ((fac2.exponent And &H7FFFFFFF) - BIAS - 1)
 
    If t < (NUM_DIGITS + 8) Then
        'The difference between the two
        'exponents indicate how many times
        'we have to multiply the mantissa
        'of FAC2 by 10 (i.e., shift it right 1 place).
        'If we have to shift more times than
        'we have digits, the result is already in FAC1.
        t = fac1.exponent - fac2.exponent
        If t > 0 And t < (NUM_DIGITS + 8) Then 'shift
 
            i = t \ 8
            While i > 0
                Call RSHIFT_8(fac2, digits)
                t = t - 8
                i = i - 1
            Wend
            If t = 7 Then
                Call RSHIFT_7(fac2, digits)
            ElseIf t = 6 Then
                Call RSHIFT_6(fac2, digits)
            ElseIf t = 5 Then
                Call RSHIFT_5(fac2, digits)
            ElseIf t = 4 Then
                Call RSHIFT_4(fac2, digits)
            ElseIf t = 3 Then
                Call RSHIFT_3(fac2, digits)
            ElseIf t = 2 Then
                Call RSHIFT_2(fac2, digits)
            ElseIf t = 1 Then
                Call RSHIFT_1(fac2, digits)
            End If
        End If
        'See if the signs of the two numbers
        'are the same.  If so, add; if not, subtract.
        If fac1.sign = fac2.sign Then 'add
            Call fpadd_aux(fac1, fac2, digits)
        Else
            Call fpsub_aux(fac1, fac2, digits)
        End If
    End If
    result = fac1
End Sub
 
Sub fpsub (result As decfloat, x As decfloat, y As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat fac1, fac2
    fac1 = x
    fac2 = y
    fac2.sign = fac2.sign Xor &H8000
    Call fpadd(fac1, fac1, fac2, digits)
    result = fac1
End Sub
 
Sub fpmul (result As decfloat, x As decfloat, y As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat fac1, fac2
    Dim As Long i, j, ex, er, den, num
    Dim As Integer64 digit, carry
    Dim As Unsigned Integer64 fac3(0 To digits)
 
    fac1 = x
    fac2 = y
    'check exponents.  if either is zero,
    'the result is zero
    If fac1.exponent = 0 Or fac2.exponent = 0 Then 'result is zero...clear fac1.
        fac1.sign = 0
        fac1.exponent = 0
        For i = 0 To digits
            Call set(fac1.mantissa, i, 0)
        Next
        'NORM_FAC1(fac1)
        result = fac1: Exit Sub
    Else
 
        If ex < 0 Then
            er = EXPO_ERR
            result = fac1: Exit Sub
        End If
 
        'clear fac3 mantissa
        For i = 0 To digits
            fac3(i) = 0
        Next
 
        den = digits
        While git(fac2.mantissa, den) = 0
            den = den - 1
        Wend
        num = digits
        While git(fac1.mantissa, num) = 0
            num = num - 1
        Wend
 
        If num < den Then
            Swap fac1, fac2
            'fac1=y
            'fac2=x
            Swap den, num
        End If
 
        For j = den To 0 Step -1
            carry = 0
            digit = git(fac2.mantissa, j)
            For i = num To 0 Step -1
                fac3(i) = fac3(i) + digit * git(fac1.mantissa, i)
            Next
            For i = num To 0 Step -1
                fac3(i) = fac3(i) + carry
                carry = fac3(i) \ 100000000
                fac3(i) = (fac3(i) Mod 100000000)
            Next
 
            For i = digits To 1 Step -1
                fac3(i) = fac3(i - 1)
            Next
            fac3(0) = carry
        Next
 
        For i = 0 To digits
            Call set(fac1.mantissa, i, fac3(i))
        Next
    End If
    'now determine exponent of result.
    'as you do...watch for overflow.
    ex = fac2.exponent - BIAS + fac1.exponent
    fac1.exponent = ex
    'determine the sign of the product
    fac1.sign = fac1.sign Xor fac2.sign
    Call NORM_FAC1(fac1, digits)
    result = fac1
End Sub
 
Sub fpmul_si (result As decfloat, x As decfloat, y As Integer64, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat fac1, fac2
    Dim As Long count, ex, er, i
    Dim As Integer64 carry, digit, prod, value
    fac1 = x
    digit = Abs(y)
    If digit > 99999999 Then
        Call si2fp(fac2, y, digits)
        Call fpmul(fac1, fac1, fac2, digits)
        result = fac1: Exit Sub
    End If
    'check exponents.  if either is zero,
    'the result is zero
    If fac1.exponent = 0 Or y = 0 Then 'result is zero...clear fac1.
        fac1.sign = 0
        fac1.exponent = 0
        For count = 0 To digits
            Call set(fac1.mantissa, count, 0)
        Next
        Call NORM_FAC1(fac1, digits)
        result = fac1: Exit Sub
    Else
        If digit = 1 Then
            If y < 0 Then
                fac1.sign = fac1.sign Xor &H8000
            End If
            result = fac1: Exit Sub
        End If
        'now determine exponent of result.
        'as you do...watch for overflow.
 
        If ex < 0 Then
            er = EXPO_ERR
            result = fac1: Exit Sub
        End If
 
        carry = 0
 
        For i = digits To 0 Step -1
            prod = digit * git(fac1.mantissa, i) + carry
            value = (prod Mod 100000000)
            Call set(fac1.mantissa, i, value)
            carry = prod \ 100000000
        Next
 
        If carry < 10 Then
            Call RSHIFT_1(fac1, digits)
            fac1.exponent = fac1.exponent + 1
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10000000)
        ElseIf carry < 100 Then
            Call RSHIFT_2(fac1, digits)
            fac1.exponent = fac1.exponent + 2
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 1000000)
        ElseIf carry < 1000 Then
            Call RSHIFT_3(fac1, digits)
            fac1.exponent = fac1.exponent + 3
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 100000)
        ElseIf carry < 10000 Then
            Call RSHIFT_4(fac1, digits)
            fac1.exponent = fac1.exponent + 4
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10000)
        ElseIf carry < 100000 Then
            Call RSHIFT_5(fac1, digits)
            fac1.exponent = fac1.exponent + 5
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 1000)
        ElseIf carry < 1000000 Then
            Call RSHIFT_6(fac1, digits)
            fac1.exponent = fac1.exponent + 6
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 100)
        ElseIf carry < 10000000 Then
            Call RSHIFT_7(fac1, digits)
            fac1.exponent = fac1.exponent + 7
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10)
        ElseIf carry < 100000000 Then
            Call RSHIFT_8(fac1, digits)
            fac1.exponent = fac1.exponent + 8
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry)
        End If
 
    End If
 
    Call NORM_FAC1(fac1, digits)
 
    If y < 0 Then
        fac1.sign = fac1.sign Xor &H8000
    End If
    result = fac1
End Sub
 
Sub fpmul_ll (result As decfloat, x As decfloat, y As Integer64, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat fac1, fac2
    Dim As Long count, ex, er, i
    Dim As Integer64 carry, digit, prod, value, n0, n1
    fac1 = x
    digit = Abs(y)
    If digit > 99999999 And digit < 10000000000000000 Then
        n0 = digit \ 100000000
        n1 = digit Mod 100000000
        Call fpmul_si(fac2, fac1, n1, digits)
        fac1.exponent = fac1.exponent + 8
        Call fpmul_si(fac1, fac1, n0, digits)
        Call fpadd(fac1, fac1, fac2, digits)
        If y < 0 Then fac1.sign = fac1.sign Xor &H8000
        result = fac1: Exit Sub
    End If
    If digit > 9999999999999999 Then
        Call str2fp(fac2, Str$(y))
        Call fpmul(fac1, fac1, fac2, digits)
        result = fac1: Exit Sub
    End If
    Call si2fp(fac2, y, digits)
 
    'check exponents.  if either is zero,
    'the result is zero
    If fac1.exponent = 0 Or y = 0 Then 'result is zero...clear fac1.
        fac1.sign = 0
        fac1.exponent = 0
        For count = 0 To digits
            Call set(fac1.mantissa, count, 0)
        Next
        Call NORM_FAC1(fac1, digits)
        result = fac1: Exit Sub
    Else
        If digit = 1 Then
            If y < 0 Then
                fac1.sign = fac1.sign Xor &H8000
            End If
            result = fac1: Exit Sub
        End If
        'now determine exponent of result.
        'as you do...watch for overflow.
 
        If ex < 0 Then
            er = EXPO_ERR
            result = fac1: Exit Sub
        End If
 
        carry = 0
 
        For i = digits To 0 Step -1
            prod = digit * git(fac1.mantissa, i) + carry
            value = (prod Mod 100000000)
            Call set(fac1.mantissa, i, value)
            carry = prod \ 100000000
        Next
 
        If carry < 10 Then
            Call RSHIFT_1(fac1, digits)
            fac1.exponent = fac1.exponent + 1
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10000000)
        ElseIf carry < 100 Then
            Call RSHIFT_2(fac1, digits)
            fac1.exponent = fac1.exponent + 2
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 1000000)
        ElseIf carry < 1000 Then
            Call RSHIFT_3(fac1, digits)
            fac1.exponent = fac1.exponent + 3
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 100000)
        ElseIf carry < 10000 Then
            Call RSHIFT_4(fac1, digits)
            fac1.exponent = fac1.exponent + 4
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10000)
        ElseIf carry < 100000 Then
            Call RSHIFT_5(fac1, digits)
            fac1.exponent = fac1.exponent + 5
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 1000)
        ElseIf carry < 1000000 Then
            Call RSHIFT_6(fac1, digits)
            fac1.exponent = fac1.exponent + 6
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 100)
        ElseIf carry < 10000000 Then
            Call RSHIFT_7(fac1, digits)
            fac1.exponent = fac1.exponent + 7
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry * 10)
        ElseIf carry < 100000000 Then
            Call RSHIFT_8(fac1, digits)
            fac1.exponent = fac1.exponent + 8
            Call set(fac1.mantissa, 0, git(fac1.mantissa, 0) + carry)
        End If
 
    End If
 
    Call NORM_FAC1(fac1, digits)
 
    If y < 0 Then
        fac1.sign = fac1.sign Xor &H8000
    End If
    result = fac1
End Sub
 
Sub fpinv (result As decfloat, m As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As Double x
    Dim As Long k, l, ex
    Dim As Long prec
    Dim As decfloat r, r2, one, n
    n = m
    l = Log((NUM_DIGITS + 8) * 0.0625) * 1.5
    Dim As String v, ts
    If n.exponent > 0 Then
        ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
    Else
        ex = 0
    End If
 
    If n.sign Then v = "-" Else v = " "
    ts = Str$(git(n.mantissa, 0))
    If Len(ts) < 8 Then
        ts = ts + String$(8 - Len(ts), "0")
    End If
    v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
    ts = Str$(git(n.mantissa, 1))
    If Len(ts) < 8 Then
        ts = String$(8 - Len(ts), "0") + ts
    End If
    v = v + ts
    x = Val(v)
    If x = 0 Then Print "Div 0": result = r: Exit Sub
    If x = 1 And ex = 0 Then
        Call str2fp(r, "1")
        result = r: Exit Sub
    ElseIf x = 1 Then
        x = 10
        ex = ex - 1
    End If
    ex = (-1) - ex
    x = 1 / x
    Call str2fp(r, Str$(x))
    r.exponent = ex + BIAS + 1
    Call set(one.mantissa, 0, 10000000)
    one.exponent = BIAS + 1
    r2 = r
    prec = 3
    For k = 1 To l
        prec = 2 * prec - 1
        Call fpmul(r2, n, r, prec)
        Call fpsub(r2, one, r2, prec)
        Call fpmul(r2, r, r2, prec)
        Call fpadd(r, r, r2, prec)
    Next
    result = r
End Sub
 
Function min& (a As Long, b As Long)
    If a < b Then min& = a Else min& = b
End Function
 
Function RealW# (w() As Double, j As Long)
    Dim wx As Double
    wx = ((w(j - 1) * 10000 + w(j)) * 10000 + w(j + 1)) * 10000
    If UBound(w) >= (j + 2) Then wx = wx + w(j + 2)
    RealW# = wx
End Function
 
Sub subtract (w() As Double, q As Long, d() As Double, ka As Long, kb As Long)
    Dim As Long j
    For j = ka To kb
        w(j) = w(j) - q * d(j - ka + 2)
    Next
End Sub
 
Sub normalize (w() As Double, ka As Long, q As Long)
    w(ka) = w(ka) + w(ka - 1) * 10000
    w(ka - 1) = q
End Sub
 
Sub finalnorm (w() As Double, kb As Long)
    Dim As Long carry, j
    For j = kb To 3 Step -1
        If w(j) < 0 Then
            carry = ((-w(j) - 1) \ 10000) + 1
        Else
            If w(j) >= 10000 Then
                carry = -(w(j) \ 10000)
            Else
                carry = 0
            End If
        End If
        w(j) = w(j) + carry * 10000
        w(j - 1) = w(j - 1) - carry
    Next
End Sub
 
Sub fpdiv (result As decfloat, x As decfloat, y As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat fac1, fac2
    Dim As Long i, er, is_power_of_ten
 
    fac1 = x
    fac2 = y
 
    If fac2.exponent = 0 Then ' if fac2 = 0, return
        ' a divide-by-zero error and
        ' bail out.
        For i = 0 To digits
            Call set(fac1.mantissa, i, 99999999)
        Next
        fac1.exponent = 99999 + BIAS + 1
        er = DIVZ_ERR
        result = fac1
        Exit Sub
    ElseIf fac1.exponent = 0 Then 'fact1=0, just return
        er = 0
        result = fac1
        Exit Sub
    Else
        'check to see if fac2 is a power of ten
        is_power_of_ten = 0
        If git(fac2.mantissa, 0) = 10000000 Then
            is_power_of_ten = 1
            For i = 1 To digits
                If git(fac2.mantissa, i) <> 0 Then
                    is_power_of_ten = 0
                    Exit For
                End If
            Next
        End If
        'if fac2 is a power of ten then all we need to do is to adjust the sign and exponent and we are finished
        If is_power_of_ten = 1 Then
            fac1.sign = fac1.sign Xor fac2.sign
            fac1.exponent = fac1.exponent - fac2.exponent + BIAS + 1
            result = fac1
            Exit Sub
        End If
 
        Dim As Double result(1 To 2 * digits + 3), n(1 To 2 * digits + 3), d(1 To 2 * digits + 3)
        Const b = 10000
        Dim As Long j, last, laststep, q, t
        Dim As Long stp
        Dim As Double xd, xn, rund
        Dim As Double w(1 To UBound(n) + 4)
 
        For j = 0 To digits
            n(2 * j + 2) = git(fac1.mantissa, j) \ 10000
            n(2 * j + 3) = git(fac1.mantissa, j) Mod 10000
            d(2 * j + 2) = git(fac2.mantissa, j) \ 10000
            d(2 * j + 3) = git(fac2.mantissa, j) Mod 10000
        Next
        n(1) = (fac1.exponent And &H7FFFFFFF) - BIAS - 1
        d(1) = (fac2.exponent And &H7FFFFFFF) - BIAS - 1
        For j = UBound(n) To UBound(w)
            w(j) = 0
        Next
        t = UBound(n) - 1
        w(1) = n(1) - d(1) + 1
        w(2) = 0
        For j = 2 To UBound(n)
            w(j + 1) = n(j)
        Next
        xd = (d(2) * b + d(3)) * b + d(4) + d(5) / b
        laststep = t + 2
        For stp = 1 To laststep
            xn = RealW(w(), (stp + 2))
            q = Int(xn / xd)
            last = min(stp + t + 1, UBound(w))
            Call subtract(w(), q, d(), (stp + 2), last)
            Call normalize(w(), (stp + 2), q)
        Next
        Call finalnorm(w(), (laststep + 1))
        If w(2) <> 0 Then laststep = laststep - 1
        rund = w(laststep + 1) / b
        If rund >= 0.5 Then w(laststep) = w(laststep) + 1
        If w(2) = 0 Then
            For j = 1 To t + 1
                result(j) = w(j + 1)
            Next
        Else
            For j = 1 To t + 1
                result(j) = w(j)
            Next
        End If
        If w(2) = 0 Then result(1) = w(1) - 1 Else result(1) = w(1)
 
        For j = 0 To digits
            Call set(fac1.mantissa, j, result(2 * j + 2) * 10000 + result(2 * j + 3))
        Next
        Call NORM_FAC1(fac1, digits)
        fac1.exponent = (result(1) + BIAS)
    End If
    fac1.sign = fac1.sign Xor fac2.sign
    result = fac1
End Sub
 
' sqrt(num)
Sub fpsqr (result As decfloat, num As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat r, r2, tmp, n, half
    Dim As Long ex, k, l, prec
    Dim As String ts, v
    Dim As Double x
    l = Log((NUM_DIGITS + 8) * 0.0625) * 1.5
    'l=estimated number of iterations needed
    'first estimate is accurate to about 16 digits
    'l is approximatly = to log2((NUM_DIGITS+9)/16)
    'NUM_DIGITS+9 because decfloat has an extra 9 guard digits
    n = num
    Call si2fp(tmp, 0, digits)
    If fpcmp(n, tmp, digits) = 0 Then
        Call si2fp(r, 0, digits)
        result = r
        Exit Sub
    End If
    Call si2fp(tmp, 1, digits)
    If fpcmp(n, tmp, digits) = 0 Then
        Call si2fp(r, 1, digits)
        result = r
        Exit Sub
    End If
    Call si2fp(tmp, 0, digits)
    If fpcmp(n, tmp, digits) < 0 Then
        Call si2fp(r, 0, digits)
        result = r
        Exit Sub
    End If
    '=====================================================================
    'hack to bypass the limitation of double exponent range
    'in case the number is larger than what a double can handle
    'for example, if the number is 2e500
    'we separate the exponent and mantissa in this case 2
    'if the exponent is odd then multiply the mantissa by 10
    'take the square root and assign it to decfloat
    'divide the exponent in half for square root
    'in this case 1.414213562373095e250
    If n.exponent > 0 Then
        ex = (n.exponent And &H7FFFFFFF) - BIAS - 1
    Else
        ex = 0
    End If
    ts = Trim$(Str$(git(n.mantissa, 0)))
    If Len(ts) < 8 Then
        ts = ts + String$(8 - Len(ts), "0")
    End If
    v = v + Left$(ts, 1) + "." + Mid$(ts, 2)
    ts = Trim$(Str$(git(n.mantissa, 1)))
    If Len(ts) < 8 Then
        ts = String$(8 - Len(ts), "0") + ts
    End If
    v = v + ts
    x = Val(v)
    If x = 0 Then Print "Div 0": result = r: Exit Sub
    If x = 1 And ex = 0 Then
        Call si2fp(r, 1, digits)
        result = r
        Exit Sub
    End If
    If Abs(ex) And 1 Then
        x = x * 10
        ex = ex - 1
    End If
    x = Sqr(x) 'approximation
    v = Trim$(Str$(x))
    k = InStr(v, ".")
    Call str2fp(r, v)
    r.exponent = ex \ 2 + BIAS + 1
    If Len(v) > 1 And k = 0 Then r.exponent = r.exponent + 1
    For k = 1 To digits
        Call set(half.mantissa, k, 0)
    Next
    Call set(half.mantissa, 0, 50000000)
    half.exponent = BIAS
    half.sign = 0
    '=====================================================================
    'Newton-Raphson method
    prec = 3
    For k = 1 To l + 1
        prec = 2 * prec - 1
        Call fpdiv(tmp, n, r, prec)
        Call fpadd(r2, r, tmp, prec)
        Call fpmul(r, r2, half, prec)
    Next
    result = r
End Sub
 
Sub fpdiv_si (result As decfloat, num As decfloat, den As Long, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat fac1
    Dim As Unsigned Integer64 carry, remder
    Dim As Integer64 i, divisor
    Dim As Integer64 quotient
    remder = 0
    divisor = Abs(den)
    fac1 = num
    If divisor = 0 Then
        Print "error: divisor = 0"
        result = fac1: Exit Sub
    End If
    If divisor > 99999999 Then
        Print "error: divisor too large"
        result = fac1: Exit Sub
    End If
 
    For i = 0 To digits
        quotient = git(fac1.mantissa, i) + remder * 100000000
        remder = quotient Mod divisor
        Call set(fac1.mantissa, i, quotient \ divisor)
    Next
    quotient = remder * 100000000
    quotient = quotient \ divisor
    carry = git(fac1.mantissa, 0)
 
    If carry = 0 Then
        Call LSHIFT_8(fac1, digits)
        fac1.exponent = fac1.exponent - 8
        Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient)
    ElseIf carry < 10 Then
        Call LSHIFT_7(fac1, digits)
        fac1.exponent = fac1.exponent - 7
        Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 10)
    ElseIf carry < 100 Then
        Call LSHIFT_6(fac1, digits)
        fac1.exponent = fac1.exponent - 6
        Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 100)
    ElseIf carry < 1000 Then
        Call LSHIFT_5(fac1, digits)
        fac1.exponent = fac1.exponent - 5
        Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 1000)
    ElseIf carry < 10000 Then
        Call LSHIFT_4(fac1, digits)
        fac1.exponent = fac1.exponent - 4
        Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 10000)
    ElseIf carry < 100000 Then
        Call LSHIFT_3(fac1, digits)
        fac1.exponent = fac1.exponent - 3
        Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 100000)
    ElseIf carry < 1000000 Then
        Call LSHIFT_2(fac1, digits)
        fac1.exponent = fac1.exponent - 2
        Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 1000000)
    ElseIf carry < 10000000 Then
        Call LSHIFT_1(fac1, digits)
        fac1.exponent = fac1.exponent - 1
        Call set(fac1.mantissa, digits, git(fac1.mantissa, digits) + quotient \ 10000000)
    End If
 
    'NORM_FAC1(fac1)
    If den < 0 Then
        fac1.sign = fac1.sign Xor &H8000
    End If
    result = fac1
End Sub
 
'fractional part of num
Sub fpfrac (result As decfloat, num As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat n
    Call fpfix(n, num)
    Call fpsub(n, num, n, digits)
    result = n
End Sub
 
'returns the positive of n
Sub fpabs (result As decfloat, n As decfloat)
    Dim As decfloat x
    x = n
    x.sign = 0
    result = x
End Sub
 
'changes the sign of n, if n is positive then n will be negative & vice versa
Sub fpneg (result As decfloat, n As decfloat)
    Dim As decfloat x
    x = n
    x.sign = x.sign Xor &H8000
    result = x
End Sub
 
'returns the negative of n regardless of the sign of n
Sub fpnegative (result As decfloat, n As decfloat)
    Dim As decfloat x
    x = n
    x.sign = &H8000
    result = x
End Sub
 
Sub fpfmod (quotient As decfloat, f_mod As decfloat, num As decfloat, denom As decfloat, digits As Long)
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat q, fm 'make copy in case the destination and source are the same
    Call fpdiv(fm, num, denom, digits): Call fpfix(q, fm)
    Call fpsub(fm, fm, q, digits): Call fpmul(fm, fm, denom, digits)
    quotient = q
    f_mod = fm
End Sub
 
Sub fpeps (result As decfloat, digits As Long)
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat ep
    Call si2fp(ep, 1, digits)
    ep.exponent = (-(NUM_DIGITS) + BIAS + 1)
    result = ep
End Sub
 
Sub fpipow (result As decfloat, x As decfloat, e As Integer64, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    'take x to an Long power
    Dim As decfloat y
    Dim As decfloat z
    Dim As Integer64 n, c
    Dim As Long i
    c = 0
    y = x
    n = Abs(e)
    z.sign = 0
    z.exponent = (BIAS + 1)
    Call set(z.mantissa, 0, 10000000)
    For i = 1 To NUM_DWORDS
        Call set(z.mantissa, i, 0)
    Next
    While n > 0
        While (n And 1) = 0
            n = n \ 2
            Call fpmul(y, y, y, digits)
            c = c + 1
        Wend
        n = n - 1
        Call fpmul(z, y, z, digits)
        c = c + 1
    Wend
    If e < 0 Then
        Call fpinv(z, z, digits)
    End If
    result = z
End Sub
 
Sub fpfactorial (result As decfloat, n As Long, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    Dim As Long i
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat f
    If n < 0 Then Print "inf": result = f: Exit Sub
    If n = 0 Or n = 1 Then
        Call si2fp(f, 1, digits)
        result = f: Exit Sub
    End If
    Call si2fp(f, 2, digits)
    If n = 2 Then result = f: Exit Sub
    For i = 3 To n
        Call fpmul_si(f, f, i, digits)
    Next
    result = f
End Sub
 
Sub fplogTaylor (result As decfloat, x As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    'taylor series
    '====================Log Guard==================
    Dim As decfloat g, zero
    Dim As Long i
    Call fpabs(g, x)
    zero.sign = 0
    zero.exponent = 0
    For i = 0 To NUM_DWORDS
        Call set(zero.mantissa, i, 0)
    Next
    If fpcmp(g, x, digits) <> 0 Then result = zero: Exit Sub
    If fpcmp(x, zero, digits) = 0 Then result = zero: Exit Sub
    '=============================================
    Dim As Long invflag
    Dim As decfloat XX, Term, Accum, x9, tmp, tmp2
    Dim As decfloat T, B, one, Q, two
 
    one.sign = 0
    one.exponent = (BIAS + 1)
    Call set(one.mantissa, 0, 10000000)
    two.sign = 0
    two.exponent = (BIAS + 1)
    Call set(two.mantissa, 0, 20000000)
    For i = 1 To NUM_DWORDS
        Call set(one.mantissa, i, 0)
        Call set(two.mantissa, i, 0)
    Next
    x9 = x
    If fpcmp(x, one, digits) < 0 Then
        invflag = 1
        Call fpdiv(x9, one, x9, digits)
    End If
    Call fpsub(T, x9, one, digits)
    Call fpadd(B, x9, one, digits)
    Call fpdiv(Accum, T, B, digits)
    Call fpdiv(Q, T, B, digits)
    tmp = Q
    Call fpmul(XX, Q, Q, digits)
    Dim As Long c
    c = 1
    Do
        c = c + 2
        tmp2 = tmp
        Call fpmul(Q, Q, XX, digits)
        Call fpdiv_si(Term, Q, c, digits)
        Call fpadd(Accum, tmp, Term, digits)
        Swap tmp, Accum
    Loop Until fpcmp(tmp, tmp2, digits) = 0
    Call fpmul_si(Accum, Accum, 2, digits)
    If invflag Then
        Call fpneg(Accum, Accum)
        result = Accum: Exit Sub
    End If
    result = Accum
End Sub
 
Sub fplog (result As decfloat, x As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    '====================Log Guard==================
    Dim As decfloat g, one, zero
    Dim As Long i, factor
    zero.sign = 0
    zero.exponent = 0
    For i = 0 To NUM_DWORDS
        Call set(zero.mantissa, i, 0)
    Next
    Call si2fp(one, 1, digits)
    Call fpabs(g, x)
    If fpcmp(g, x, digits) <> 0 Then result = zero: Exit Sub
    If fpcmp(x, zero, digits) = 0 Then result = zero: Exit Sub
    If fpcmp(x, one, digits) = 0 Then result = zero: Exit Sub
    '=============================================
    Dim As decfloat approx, ans, logx
    approx = x
    factor = 4096
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fpsqr(approx, approx, digits)
    Call fplogTaylor(logx, approx, digits)
    Call fpmul_si(ans, logx, factor, digits)
    result = ans
End Sub
 
Sub fpexp (result As decfloat, x As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    'taylor series
    Dim As decfloat fac, x9, temp, accum, p, term
    Dim As Long i, c
    Call si2fp(temp, 0, digits)
 
    fac.sign = 0
    fac.exponent = (BIAS + 1)
    Call set(fac.mantissa, 0, 10000000)
    For i = 1 To NUM_DWORDS
        Call set(fac.mantissa, i, 0)
    Next
    If fpcmp(x, temp, digits) = 0 Then result = fac: Exit Sub
    c = 1
    Call fpdiv_si(x9, x, 8192, digits) 'fpdiv_si(x, 67108864) '
    p = x9
    Call fpadd(accum, fac, x9, digits) '1 + x
 
    Do
        c = c + 1
        temp = accum
        Call fpdiv_si(fac, fac, c, digits)
        Call fpmul(p, p, x9, digits)
        Call fpmul(term, p, fac, digits)
        Call fpadd(accum, temp, term, digits)
    Loop Until fpcmp(accum, temp, digits) = 0
    For i = 1 To 13
        Call fpmul(accum, accum, accum, digits)
    Next
    result = accum
End Sub
 
Sub fpexp_cf (result As decfloat, x As decfloat, n As Long, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
 
    '                       x
    '    cf =   ------------------------------
    '                         x^2
    '           2 + -------------------------
    '                           x^2
    '              6 + ----------------------
    '                             x^2
    '                 10 + ------------------
    '                               x^2
    '                     14 + --------------
    '                                 x^2
    '                          18 + ---------
    '                                   x^2
    '                              22 + -----
    '                                   26 + ...
    '           (1 + cf)
    ' exp(x) = ----------
    '           (1 - cf)
 
    Dim i As Long
    Dim As decfloat s, x1, x2, tmp, tmp2
    Call fpdiv_si(x1, x, 8192, digits) ' 2^13
    Call fpmul(x2, x1, x1, digits)
    Call si2fp(s, 4 * n + 6, digits)
    For i = n To 0 Step -1
        Call si2fp(tmp, 4 * i + 2, digits)
        Call fpdiv(s, x2, s, digits)
        Call fpadd(s, s, tmp, digits)
    Next
    Call fpdiv(s, x1, s, digits)
    Call si2fp(tmp, 1, digits)
    tmp2 = tmp
    Call fpadd(tmp, tmp, s, digits)
    Call fpsub(tmp2, tmp2, s, digits)
    Call fpdiv(s, tmp, tmp2, digits)
    For i = 1 To 13 ' s^13
        Call fpmul(s, s, s, digits)
    Next
    result = s
End Sub
 
Sub fplog_r (result As decfloat, x As decfloat, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    '====================Log Guard==================
    Dim As decfloat x2, y, tmp
    Dim As Long i
    Dim As Long ex, terms, prec, n
    Dim As Double t, t2
 
    tmp.sign = 0
    tmp.exponent = (BIAS + 1)
    Call set(tmp.mantissa, 0, 10000000)
    For i = 1 To NUM_DWORDS
        Call set(tmp.mantissa, i, 0)
    Next
    If fpcmp(x, tmp, digits) = 0 Then result = x2: Exit Sub
    ex = (x.exponent And &H7FFFFFFF) - BIAS - 1
    t = git(x.mantissa, 0) + git(x.mantissa, 1) / 100000000
    t = t / 10000000
    t2 = Log(t) + Log(10) * ex
 
    n = Log((digits + 1) * 0.5) * 1.5
 
    Call str2fp(x2, Str$(t2))
    prec = 3
    terms = 2
    prec = 2 * prec - 1
    Call fpexp_cf(y, x2, terms, prec)
    Call fpsub(tmp, y, x, prec)
    Call fpdiv(tmp, tmp, y, prec)
    Call fpsub(x2, x2, tmp, prec)
    If digits > 4 And digits <= 9 Then
        prec = 2 * prec - 1
        terms = 2 * terms + 1
        Call fpexp_cf(y, x2, terms, prec)
        Call fpsub(tmp, y, x, prec)
        Call fpdiv(tmp, tmp, y, prec)
        Call fpsub(x2, x2, tmp, prec)
    ElseIf digits > 9 And digits < 18 Then
        For i = 1 To 3
            prec = 2 * prec - 1
            terms = 2 * terms + 1
            Call fpexp_cf(y, x2, terms, prec)
            Call fpsub(tmp, y, x, prec)
            Call fpdiv(tmp, tmp, y, prec)
            Call fpsub(x2, x2, tmp, prec)
        Next
    ElseIf digits > 17 And digits < 34 Then
        For i = 1 To 4
            prec = 2 * prec - 1
            terms = 2 * terms + 1
            Call fpexp_cf(y, x2, terms, prec)
            Call fpsub(tmp, y, x, prec)
            Call fpdiv(tmp, tmp, y, prec)
            Call fpsub(x2, x2, tmp, prec)
        Next
    ElseIf digits > 33 And digits < 65 Then
        For i = 1 To 5
            prec = 2 * prec - 1
            terms = 2 * terms + 1
            Call fpexp_cf(y, x2, terms, prec)
            Call fpsub(tmp, y, x, prec)
            Call fpdiv(tmp, tmp, y, prec)
            Call fpsub(x2, x2, tmp, prec)
        Next
    ElseIf digits > 64 Then
        For i = 1 To 6
            prec = 2 * prec - 1
            terms = 2 * terms + 1
            Call fpexp_cf(y, x2, terms, prec)
            Call fpsub(tmp, y, x, prec)
            Call fpdiv(tmp, tmp, y, prec)
            Call fpsub(x2, x2, tmp, prec)
        Next
    End If
    result = x2
End Sub
 
Sub fppow (result As decfloat, lhs As decfloat, rhs As decfloat)
    Dim As decfloat lhs2
    Call fplog(lhs2, lhs, NUM_DWORDS)
    Call fpmul(lhs2, lhs2, rhs, NUM_DWORDS)
    Call fpexp(lhs2, lhs2, NUM_DWORDS)
    result = lhs2
End Sub
 
Sub fpnroot (result As decfloat, x As decfloat, p As Long, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat ry, tmp, tmp2
    Dim As Double t, t2
    Dim As Long i, ex, l, prec
 
    l = Log((NUM_DIGITS + 8) * 0.0625) * 1.5
    ex = (x.exponent And &H7FFFFFFF) - BIAS - 1
    t = git(x.mantissa, 0) + git(x.mantissa, 1) / 100000000
    t = t / 10000000
    t2 = Log(t) / p + Log(10) * ex / p
    t2 = Exp(t2)
    Call str2fp(ry, Str$(t2))
    prec = 3
 
    Call fpipow(tmp, ry, p - 1, prec)
    Call fpdiv(tmp, x, tmp, prec)
    Call fpmul_si(tmp2, ry, p - 1, prec)
    Call fpadd(tmp2, tmp2, tmp, prec)
    Call fpdiv_si(ry, tmp2, p, prec)
    For i = 1 To l
        prec = 2 * prec - 1
        Call fpipow(tmp, ry, p - 1, prec)
        Call fpdiv(tmp, x, tmp, prec)
        Call fpmul_si(tmp2, ry, p - 1, prec)
        Call fpadd(tmp2, tmp2, tmp, prec)
        Call fpdiv_si(ry, tmp2, p, prec)
    Next
    result = ry
End Sub
 
Sub pi_brent_salamin (pi_bs As decfloat, digits_in As Unsigned Long)
    Dim As Unsigned Long digits
    digits = digits_in
    If digits > NUMBER_OF_DIGITS Then digits = NUMBER_OF_DIGITS
 
    Dim As Long limit2
    Dim As decfloat c0, c1, c2, c05
    Dim As decfloat a, b, sum
    Dim As decfloat ak, bk, ck
    Dim As decfloat ab, asq
    Dim As decfloat pow2, tmp
 
    limit2 = -digits + BIAS + 1
    Call si2fp(c0, 0, NUMBER_OF_DIGITS): ak = c0: bk = c0: ab = c0: asq = c0
    Call si2fp(c1, 1, NUMBER_OF_DIGITS): a = c1: ck = c1: pow2 = c1
    Call si2fp(c2, 2, NUMBER_OF_DIGITS): b = c2
    Call str2fp(c05, ".5"): sum = c05
    Call si2fp(pi_bs, 3, NUMBER_OF_DIGITS)
 
    Call fpsqr(b, b, NUMBER_OF_DIGITS)
    Call fpdiv(b, c1, b, NUMBER_OF_DIGITS)
    While fpcmp(ck, c0, NUMBER_OF_DIGITS) <> 0 And ck.exponent > limit2
        Call fpadd(ak, a, b, NUMBER_OF_DIGITS)
        Call fpmul(ak, c05, ak, NUMBER_OF_DIGITS)
        Call fpmul(ab, a, b, NUMBER_OF_DIGITS)
        Call fpsqr(bk, ab, NUMBER_OF_DIGITS)
        Call fpmul(asq, ak, ak, NUMBER_OF_DIGITS)
        Call fpsub(ck, asq, ab, NUMBER_OF_DIGITS)
        Call fpmul(pow2, pow2, c2, NUMBER_OF_DIGITS)
        Call fpmul(tmp, pow2, ck, NUMBER_OF_DIGITS)
        Call fpsub(sum, sum, tmp, NUMBER_OF_DIGITS)
        a = ak: b = bk
    Wend
    Call fpdiv(tmp, asq, sum, NUMBER_OF_DIGITS)
    Call fpmul(pi_bs, c2, tmp, NUMBER_OF_DIGITS)
End Sub
 
Sub fpsin (result As decfloat, x As decfloat, digits_in As Unsigned Long)
    Dim As decfloat XX, Term, Accum, p, temp2, fac, x_2
    Dim As decfloat pi2, circ, Ab
    Dim As Long precision
    precision = digits_in
 
    x_2 = x
    pi2 = pi_dec
    Call fpmul_si(circ, pi2, 2, precision)
    Call fpabs(Ab, x)
    If fpcmp(Ab, circ, precision) > 0 Then
        '======== CENTRALIZE ==============
        'floor/ceil to centralize
        Dim As decfloat tmp, tmp2
        If precision > 20 Then
            pi2 = pi_dec
        End If
        Call fpmul_si(pi2, pi2, 2, precision) 'got 2*pi
        Call fpdiv(tmp, x_2, pi2, precision)
        tmp2 = tmp
        Call fpfix(tmp, tmp) 'int part
        Call fpsub(tmp, tmp2, tmp, precision) 'frac part
        Call fpmul(tmp, tmp, pi2, precision)
        x_2 = tmp
    End If
 
    Dim As Long lm, limit2, i
    Dim As decfloat factor
    lm = precision
    limit2 = Int(-0.45344993886092585968 + 0.022333002852398072433 * lm + 5.0461814408333079844E-7 * lm * lm - 4.2338453039804235772E-11 * lm * lm * lm)
 
    Call si2fp(factor, 5, precision)
    Call fpipow(factor, factor, limit2, precision)
    Call fpdiv(x_2, x_2, factor, precision) 'x_=x_/5^limit2
    '==================================
    Dim As Long sign(3): sign(3) = 1
    Dim As Long c: c = 1
 
    Accum = x_2
    Call si2fp(fac, 1, precision)
    p = x_2
    Call fpmul(XX, x_2, x_2, precision)
    Do
        c = c + 2
        temp2 = Accum
        Call fpmul_si(fac, fac, c * (c - 1), precision)
        Call fpmul(p, p, XX, precision)
        Call fpdiv(Term, p, fac, precision)
        If sign(c And 3) Then
            Call fpsub(Accum, temp2, Term, precision)
        Else
            Call fpadd(Accum, temp2, Term, precision)
        End If
    Loop Until fpcmp(Accum, temp2, precision) = 0
    'multiply the result by 5^limit2
 
    For i = 1 To limit2
        Call fpmul(p, Accum, Accum, precision)
        Call fpmul(temp2, Accum, p, precision)
        '*** sin(5*x) = 5 * sin(x) - 20 * sin(x)^3 + 16 * sin(x)^5
        Call fpmul_si(Accum, Accum, 5, precision)
        Call fpmul_si(Term, temp2, 20, precision)
        Call fpmul_si(XX, temp2, 16, precision)
        Call fpmul(XX, XX, p, precision)
        Call fpsub(Accum, Accum, Term, precision)
        Call fpadd(Accum, Accum, XX, precision)
    Next i
    result = Accum
End Sub
 
Sub fpcos (result As decfloat, z As decfloat, digits_in As Unsigned Long)
    Dim As decfloat x_2, pi2
    Dim As Long precision
    precision = digits_in
    Call fpdiv_si(pi2, pi_dec, 2, precision)
    Call fpsub(x_2, pi2, z, precision)
 
    Call fpsin(result, x_2, precision)
End Sub
 
Sub fptan (result As decfloat, z As decfloat, digits_in As Unsigned Long)
    Dim As decfloat x_2, s, c
    Dim As Long precision
    precision = digits_in
    x_2 = z
    Call fpsin(s, x_2, precision)
    x_2 = z
    Call fpcos(c, x_2, precision)
    Call fpdiv(result, s, c, precision)
End Sub
 
Sub fpatn (result As decfloat, x As decfloat, digits_in As Unsigned Long)
    Dim As Long precision, z
    precision = digits_in
    Dim As Long sign(3): sign(3) = 1
    Dim As Unsigned Long c: c = 1
    Dim As decfloat XX, Term, Accum, strC, x_2, mt, mt2, p
    Dim As decfloat decnum, one, decnum2, factor
 
    decnum2 = x
    decnum2.sign = 0
    Call si2fp(one, 1, precision)
    If fpcmp(decnum2, one, precision) = 0 Then
        Call fpdiv_si(result, pi_dec, 4, precision)
        result.sign = x.sign
        Exit Sub
    End If
    decnum2.sign = x.sign
    Dim As Long limit2: limit2 = 16
    Call si2fp(factor, SHL(2, limit2 - 1), precision)
    For z = 1 To limit2
        Call fpmul(decnum, decnum2, decnum2, precision)
        Call fpadd(decnum, decnum, one, precision)
        Call fpsqr(decnum, decnum, precision)
        Call fpadd(decnum, decnum, one, precision)
        Call fpdiv(decnum, decnum2, decnum, precision)
        decnum2 = decnum
    Next z
 
    mt = decnum
    x_2 = decnum
    p = decnum
    Call fpmul(XX, x_2, x_2, precision)
    Do
        c = c + 2
        mt2 = mt
        Call si2fp(strC, c, precision)
        Call fpmul(p, p, XX, precision)
        Call fpdiv(Term, p, strC, precision)
        If sign(c And 3) Then
            Call fpsub(Accum, mt, Term, precision)
        Else
            Call fpadd(Accum, mt, Term, precision)
        End If
        Swap mt, Accum
    Loop Until fpcmp(mt, mt2, precision) = 0
    Call fpmul(result, factor, mt, precision)
End Sub
 
Sub fpasin (result As decfloat, x As decfloat, digits_in As Unsigned Long)
    Dim As Long precision
    precision = digits_in
    Dim As Double num
    ' ARCSIN = ATN(x / SQR(-x * x + 1))
    '============= ARCSIN GUARD =========
    num = fp2dbl(x)
    If num > 1 Then Exit Sub
    If num < -1 Then Exit Sub
    '========================
    Dim As decfloat one, T, B, term1, minusone
    Call si2fp(one, 1, precision)
    Call si2fp(minusone, -1, precision)
    T = x
    Call fpmul(B, x, x, precision) 'x*x
    'for 1 and -1
    If fpcmp(B, one, precision) = 0 Then
        Dim As decfloat two, atn1
        Call si2fp(two, 2, precision)
        Call fpatn(atn1, one, precision)
        If fpcmp(x, minusone, precision) = 0 Then
            Call fpmul(two, two, atn1, precision)
            Call fpmul(result, two, minusone, precision)
            Exit Sub
        Else
            Call fpmul(result, two, atn1, precision)
            Exit Sub
        End If
    End If
    Call fpsub(B, one, B, precision) '1-x*x
    Call fpsqr(B, B, precision) 'sqr(1-x*x)
    Call fpdiv(term1, T, B, precision)
    Call fpatn(result, term1, precision)
End Sub
 
Sub fpacos (result As decfloat, x As decfloat, digits_in As Unsigned Long)
    Dim As Long precision
    precision = digits_in
    Dim As decfloat one, minusone, two, atn1, tail, T, B, term1, atnterm1 ',_x,temp
    Dim As Double num
    'ARCCOS = ATN(-x / SQR(-x * x + 1)) + 2 * ATN(1)
    '============= ARCCOS GUARD =========
    num = fp2dbl(x)
    If num > 1 Then Exit Sub
    If num < -1 Then Exit Sub
    '========================
    Call si2fp(one, 1, precision)
    Call si2fp(minusone, -1, precision)
    Call si2fp(two, 2, precision)
    Call fpatn(atn1, one, precision)
    Call fpmul(tail, two, atn1, precision) '2*atn(1)
    Call fpmul(T, minusone, x, precision) '-x
    Call fpmul(B, x, x, precision) 'x*x
    If fpcmp(B, one, precision) = 0 Then
        'for 1 and -1
        If fpcmp(x, minusone, precision) = 0 Then
            Dim As decfloat four
            Call si2fp(four, 4, precision)
            Call fpmul(result, four, atn1, precision)
            Exit Sub
        Else
            Call si2fp(result, 0, precision)
            Exit Sub
        End If
    End If
    Call fpsub(B, one, B, precision) '1-x*x
    Call fpsqr(B, B, precision) 'sqr(1-x*x)
    Call fpdiv(term1, T, B, precision)
    Call fpatn(atnterm1, term1, precision)
    Call fpadd(result, atnterm1, tail, precision)
End Sub
 

test-decfloat.bas

Code: QB64: [Select]

'you must use $NoPrefix
$NoPrefix
$Console:Only
Dest Console
 
Option Explicit
 
Const NUMBER_OF_DIGITS = 80 ' you must have NUMBER_OF_DIGITS defined -- not to exceed 1024
'                             before including decfloat.bi
'$include: 'decfloat.bi'
 
'======================================================================
' your code follows
 
Dim As decfloat x, y, z
Dim As Long i, n
n = 5
Print "     x                y=sin(x)                                     asin(y)"
For i = 1 To n
    Call si2fp(x, i, NUMBER_OF_DIGITS)
    Call fpdiv_si(x, x, n, NUMBER_OF_DIGITS)
    Call fpsin(y, x, NUMBER_OF_DIGITS)
    Call fpasin(z, y, NUMBER_OF_DIGITS)
    Print fp2str(x, 10), fp2str(y, 30), fp2str(z, 30)
Next i
Print "     x                y=cos(x)                                     acos(y)"
For i = 1 To n
    Call si2fp(x, i, NUMBER_OF_DIGITS)
    Call fpdiv_si(x, x, n, NUMBER_OF_DIGITS)
    Call fpcos(y, x, NUMBER_OF_DIGITS)
    Call fpacos(z, y, NUMBER_OF_DIGITS)
    Print fp2str(x, 10), fp2str(y, 30), fp2str(z, 30)
Next i
Print "     x                y=tan(x)                                     atn(y)"
For i = 1 To n
    Call si2fp(x, i, NUMBER_OF_DIGITS)
    Call fpdiv_si(x, x, n, NUMBER_OF_DIGITS)
    Call fptan(y, x, NUMBER_OF_DIGITS)
    Call fpatn(z, y, NUMBER_OF_DIGITS)
    Print fp2str(x, 10), fp2str(y, 30), fp2str(z, 30)
Next i
 
' end of your code
' you must include decfloat.bm
'$include: 'decfloat.bm'
 

output

Code: [Select]

     x                y=sin(x)                                     asin(y)
 0.2000000000       0.198669330795061215459412627118        0.199999999999999999999999999999
 0.4000000000       0.389418342308650491666311756795        0.399999999999999999999999999999
 0.6000000000       0.564642473395035357200945445658        0.600000000000000000000000000000
 0.8000000000       0.717356090899522761627174610581        0.799999999999999999999999999999
 1.000000000        0.841470984807896506652502321630        1.00000000000000000000000000000
     x                y=cos(x)                                     acos(y)
 0.2000000000       0.980066577841241631124196516748        0.199999999999999999999999999999
 0.4000000000       0.921060994002885082798526732051        0.399999999999999999999999999999
 0.6000000000       0.825335614909678297240952498955        0.599999999999999999999999999999
 0.8000000000       0.696706709347165420920749981642        0.800000000000000000000000000000
 1.000000000        0.540302305868139717400936607442        0.999999999999999999999999999999
     x                y=tan(x)                                     atn(y)
 0.2000000000       0.202710035508672483321358271647        0.199999999999999999999999999999
 0.4000000000       0.422793218738161761981635427165        0.399999999999999999999999999999
 0.6000000000       0.684136808341692317070925417463        0.600000000000000000000000000000
 0.8000000000       1.02963855705036401274636117282         0.799999999999999999999999999999
 1.000000000        1.55740772465490223050697480745         1.00000000000000000000000000000

Press any key to continue

News:

Author Topic: decfloat 2.0 (Read 10363 times)

jack

decfloat 2.0

bplus

Re: decfloat 2.0

FellippeHeitor

Re: decfloat 2.0

jack

Re: decfloat 2.0

jack

Re: decfloat 2.0

bplus

Re: decfloat 2.0

jack

Re: decfloat 2.0

bplus

Re: decfloat 2.0

jack

Re: decfloat 2.0

SMcNeill

Re: decfloat 2.0

bplus

Re: decfloat 2.0

FellippeHeitor

Re: decfloat 2.0

jack

Re: decfloat 2.0

Cobalt

Re: decfloat 2.0

jack

Re: decfloat 2.0