Announcement

Collapse
No announcement yet.

arbitrary length number math

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • #21
    here my code, all seems to work except pi_brent_salamin
    Code:
    #COMPILE EXE
    #DIM ALL
    
    %NUMBER_OF_DIGITS = 35660
    %NUM_DIGITS = %NUMBER_OF_DIGITS
    %NUM_DWORDS = %NUM_DIGITS \ 8
    %BIAS = 1073741824 '2 ^ 30
    %MANTISSA_BYTES = (%NUM_DWORDS + 1) * 8 * 4
    
    TYPE decfloat
        sign AS LONG
        exponent AS DWORD
        mantissa(%NUM_DWORDS) AS DWORD
    END TYPE
    
    ' Error definitions
    
    %DIVZ_ERR = 1 'Divide by zero
    %EXPO_ERR = 2 'Exponent overflow error
    
    FUNCTION PBMAIN () AS LONG
    
    DIM m AS decfloat, x AS decfloat, y AS decfloat, z AS decfloat, v AS decfloat
    DIM s AS STRING
    DIM t AS DOUBLE
    DIM i AS QUAD
    
    str2fp x, "9."+STRING$(5999, "9")
    str2fp y, "9."+STRING$(5999, "9")
    t=TIMER
    'for i=1 to 1000
        'si2fp x, i, %NUM_DWORDS
        'fpsqr x, x, %NUM_DWORDS
        'fpipow x, x, 2, %NUM_DWORDS
        '? fp2str(x, %NUMBER_OF_DIGITS)
        'fpmul z, x, y, %NUMBER_OF_DIGITS
    'next
    fpfactorial z, 10000, %NUMBER_OF_DIGITS
    t=TIMER-t
    ? fp2str(z, %NUMBER_OF_DIGITS)
    PRINT "elapsed time ";t;" seconds"
    
    WAITKEY$
    
    END FUNCTION
    
    FUNCTION fp2str_exp$ (n AS decfloat, places AS LONG)
        DIM i AS LONG, ex AS LONG
        DIM v AS STRING, f AS STRING, ts AS STRING
        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$(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$(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 ip AS decfloat
        DIM ex AS LONG, ex2 AS LONG, j AS LONG, k AS LONG
    
        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
            ip.mantissa(ex2)=num.mantissa(ex2)
        WEND
        IF j = 1 THEN
            ip.mantissa(k)= 10000000 * (num.mantissa(k) \ 10000000)
        ELSEIF j = 2 THEN
            ip.mantissa(k)= 1000000 * (num.mantissa(k) \ 1000000)
        ELSEIF j = 3 THEN
            ip.mantissa(k)= 100000 * (num.mantissa(k) \ 100000)
        ELSEIF j = 4 THEN
            ip.mantissa(k)= 10000 * (num.mantissa(k) \ 10000)
        ELSEIF j = 5 THEN
            ip.mantissa(k)= 1000 * (num.mantissa(k) \ 1000)
        ELSEIF j = 6 THEN
            ip.mantissa(k)= 100 * (num.mantissa(k) \ 100)
        ELSEIF j = 7 THEN
            ip.mantissa(k)= 10 * (num.mantissa(k) \ 10)
        ELSEIF j = 8 THEN
            ip.mantissa(k)= num.mantissa(k)
        END IF
        ip.exponent = ex + %BIAS
        ip.sign = num.sign
        result = ip
    END SUB
    
    FUNCTION fpfix_is_odd& (num AS decfloat)
        DIM ex AS LONG, j AS LONG, k AS LONG
    
        ex = (num.exponent AND &H7FFFFFFF) - %BIAS
        IF ex < 1 THEN
            fpfix_is_odd = 0: EXIT FUNCTION
        END IF
        IF ex >= (%NUM_DIGITS) THEN
            PRINT "error in function fpfix_is_odd"
            fpfix_is_odd = 99999999: EXIT FUNCTION
        END IF
        k = ex \ 8
        j = ex MOD 8
    
        IF j = 1 THEN
            fpfix_is_odd = (num.mantissa(k) \ 10000000) AND 1: EXIT FUNCTION
        ELSEIF j = 2 THEN
            fpfix_is_odd = (num.mantissa(k) \ 1000000) AND 1: EXIT FUNCTION
        ELSEIF j = 3 THEN
            fpfix_is_odd = (num.mantissa(k) \ 100000) AND 1: EXIT FUNCTION
        ELSEIF j = 4 THEN
            fpfix_is_odd = (num.mantissa(k) \ 10000) AND 1: EXIT FUNCTION
        ELSEIF j = 5 THEN
            fpfix_is_odd = (num.mantissa(k) \ 1000) AND 1: EXIT FUNCTION
        ELSEIF j = 6 THEN
            fpfix_is_odd = (num.mantissa(k) \ 100) AND 1: EXIT FUNCTION
        ELSEIF j = 7 THEN
            fpfix_is_odd = (num.mantissa(k) \ 10) AND 1: EXIT FUNCTION
        ELSEIF j = 8 THEN
            fpfix_is_odd = num.mantissa(k) AND 1: EXIT FUNCTION
        END IF
        fpfix_is_odd = 0
    END FUNCTION
    
    FUNCTION fp2dbl# (n AS decfloat)
        DIM ex AS LONG
        DIM v AS STRING, f AS STRING, ts AS STRING
        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$(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$(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 i AS LONG, ex AS LONG
        DIM v AS STRING, ts AS STRING, s AS STRING
        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$(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$(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 ex AS LONG
        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 j AS LONG, s AS LONG, d AS LONG, e AS LONG, ep AS LONG
        DIM ex AS LONG, es AS LONG, i AS LONG, f AS LONG, fp AS LONG, fln AS LONG
        DIM c AS STRING, f1 AS STRING, f2 AS STRING, f3 AS STRING, ts AS STRING
        DIM ulng AS DWORD
        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)
            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 QUAD, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM fac1 AS decfloat
        DIM i AS LONG
        DIM n AS QUAD
        DIM ndw AS LONG
        n = ABS(m)
        ndw = n
        IF n > 9999999999999999 THEN
            str2fp fac1, STR$(m)
            result = fac1: EXIT SUB
        END IF
    
        FOR i = 1 TO digits
            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
                fac1.mantissa(0) = n * 10000000
                fac1.exponent = fac1.exponent + 1
            ELSEIF n < 100 THEN
                fac1.mantissa(0) = n * 1000000
                fac1.exponent = fac1.exponent + 2
            ELSEIF n < 1000 THEN
                fac1.mantissa(0) = n * 100000
                fac1.exponent = fac1.exponent + 3
            ELSEIF n < 10000 THEN
                fac1.mantissa(0) = n * 10000
                fac1.exponent = fac1.exponent + 4
            ELSEIF n < 100000 THEN
                fac1.mantissa(0) = n * 1000
                fac1.exponent = fac1.exponent + 5
            ELSEIF n < 1000000 THEN
                fac1.mantissa(0) = n * 100
                fac1.exponent = fac1.exponent + 6
            ELSEIF n < 10000000 THEN
                fac1.mantissa(0) = n * 10
                fac1.exponent = fac1.exponent + 7
            ELSEIF n < 100000000 THEN
                fac1.mantissa(0) = ndw
                fac1.exponent = fac1.exponent + 8
            END IF
        END IF
        IF n > 99999999 THEN
            fac1.exponent = fac1.exponent + 8
            IF n < 1000000000 THEN
                fac1.mantissa(0) = n \ 10
                fac1.mantissa(1) = (n MOD 10) * 10000000
                fac1.exponent = fac1.exponent + 1
            ELSEIF n < 100000000000 THEN
                fac1.mantissa(0) = n \ 100
                fac1.mantissa(1) = (n MOD 100) * 1000000
                fac1.exponent = fac1.exponent + 2
            ELSEIF n < 1000000000000 THEN
                fac1.mantissa(0) = n \ 1000
                fac1.mantissa(1) = (n MOD 1000) * 100000
                fac1.exponent = fac1.exponent + 3
            ELSEIF n < 10000000000000 THEN
                fac1.mantissa(0) = n \ 10000
                fac1.mantissa(1) = (n MOD 10000) * 10000
                fac1.exponent = fac1.exponent + 4
            ELSEIF n < 100000000000000 THEN
                fac1.mantissa(0) = n \ 100000
                fac1.mantissa(1) = (n MOD 100000) * 1000
                fac1.exponent = fac1.exponent + 5
            ELSEIF n < 1000000000000000 THEN
                fac1.mantissa(0) = n \ 1000000
                fac1.mantissa(1) = (n MOD 1000000) * 100
                fac1.exponent = fac1.exponent + 6
            ELSEIF n < 10000000000000000 THEN
                fac1.mantissa(0) = n \ 10000000
                fac1.mantissa(1) = (n MOD 10000000) * 10
                fac1.exponent = fac1.exponent + 7
            ELSEIF n < 100000000000000000 THEN
                fac1.mantissa(0) = n \ 100000000
                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 digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
    
        DIM v1 AS DWORD, v2 AS DWORD
        DIM i AS LONG
        FOR i = digits TO 1 STEP -1
            v1 = mantissa.mantissa(i) \ 10
            v2 = mantissa.mantissa(i - 1) MOD 10
            v2 = v2 * 10000000 + v1
            mantissa.mantissa(i) = v2
        NEXT
        mantissa.mantissa(0) = mantissa.mantissa(0) \ 10
    END SUB
    
    SUB LSHIFT_1 (mantissa AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
    
        DIM v1 AS DWORD, v2 AS DWORD
        DIM i AS LONG
        FOR i = 0 TO digits - 1
            v1 = mantissa.mantissa(i) MOD 10000000
            v2 = mantissa.mantissa(i + 1) \ 10000000
            mantissa.mantissa(i) = v1 * 10 + v2
            mantissa.mantissa(i + 1) = mantissa.mantissa(i + 1) MOD 10000000
        NEXT
        mantissa.mantissa(digits) = 10 * (mantissa.mantissa(digits) MOD 10000000)
    END SUB
    
    SUB RSHIFT_2 (mantissa AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
    
        DIM v1 AS DWORD, v2 AS DWORD
        DIM i AS LONG
        FOR i = digits TO 1 STEP -1
            v1 = mantissa.mantissa(i) \ 100
            v2 = mantissa.mantissa(i - 1) MOD 100
            v2 = v2 * 1000000 + v1
            mantissa.mantissa(i) = v2
        NEXT
        mantissa.mantissa(0) = mantissa.mantissa(0) \ 100
    END SUB
    
    SUB LSHIFT_2 (mantissa AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
    
        DIM v1 AS DWORD, v2 AS DWORD
        DIM i AS LONG
        FOR i = 0 TO digits - 1
            v1 = mantissa.mantissa(i) MOD 1000000
            v2 = mantissa.mantissa(i + 1) \ 1000000
            mantissa.mantissa(i) = v1 * 100 + v2
            mantissa.mantissa(i + 1) = mantissa.mantissa(i + 1) MOD 1000000
        NEXT
        mantissa.mantissa(digits) = 100 * (mantissa.mantissa(digits) MOD 1000000)
    END SUB
    
    SUB RSHIFT_3 (mantissa AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
    
        DIM v1 AS DWORD, v2 AS DWORD
        DIM i AS LONG
        FOR i = digits TO 1 STEP -1
            v1 = mantissa.mantissa(i) \ 1000
            v2 = mantissa.mantissa(i - 1) MOD 1000
            v2 = v2 * 100000 + v1
            mantissa.mantissa(i) = v2
        NEXT
        mantissa.mantissa(0) = mantissa.mantissa(0) \ 1000
    END SUB
    
    SUB LSHIFT_3 (mantissa AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
    
        DIM v1 AS DWORD, v2 AS DWORD
        DIM i AS LONG
        FOR i = 0 TO digits - 1
            v1 = mantissa.mantissa(i) MOD 100000
            v2 = mantissa.mantissa(i + 1) \ 100000
            mantissa.mantissa(i) = v1 * 1000 + v2
            mantissa.mantissa(i + 1) = mantissa.mantissa(i + 1) MOD 100000
        NEXT
        mantissa.mantissa(digits) = 1000 * (mantissa.mantissa(digits) MOD 100000)
    END SUB
    
    SUB RSHIFT_4 (mantissa AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
    
        DIM v1 AS DWORD, v2 AS DWORD
        DIM i AS LONG
        FOR i = digits TO 1 STEP -1
            v1 = mantissa.mantissa(i) \ 10000
            v2 = mantissa.mantissa(i - 1) MOD 10000
            v2 = v2 * 10000 + v1
            mantissa.mantissa(i) = v2
        NEXT
        mantissa.mantissa(0) = mantissa.mantissa(0) \ 10000
    END SUB
    
    SUB LSHIFT_4 (mantissa AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
    
        DIM v1 AS DWORD, v2 AS DWORD
        DIM i AS LONG
        FOR i = 0 TO digits - 1
            v1 = mantissa.mantissa(i) MOD 10000
            v2 = mantissa.mantissa(i + 1) \ 10000
            mantissa.mantissa(i) = v1 * 10000 + v2
            mantissa.mantissa(i + 1) = mantissa.mantissa(i + 1) MOD 10000
        NEXT
        mantissa.mantissa(digits) = 10000 * (mantissa.mantissa(digits) MOD 10000)
    END SUB
    
    SUB RSHIFT_5 (mantissa AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
    
        DIM v1 AS DWORD, v2 AS DWORD
        DIM i AS LONG
        FOR i = digits TO 1 STEP -1
            v1 = mantissa.mantissa(i) \ 100000
            v2 = mantissa.mantissa(i - 1) MOD 100000
            v2 = v2 * 1000 + v1
            mantissa.mantissa(i) = v2
        NEXT
        mantissa.mantissa(0) = mantissa.mantissa(0) \ 100000
    END SUB
    
    SUB LSHIFT_5 (mantissa AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
    
        DIM v1 AS DWORD, v2 AS DWORD
        DIM i AS LONG
        FOR i = 0 TO digits - 1
            v1 = mantissa.mantissa(i) MOD 1000
            v2 = mantissa.mantissa(i + 1) \ 1000
            mantissa.mantissa(i) = v1 * 100000 + v2
            mantissa.mantissa(i + 1) = mantissa.mantissa(i + 1) MOD 1000
        NEXT
        mantissa.mantissa(digits) = 100000 * (mantissa.mantissa(digits) MOD 1000)
    END SUB
    
    SUB RSHIFT_6 (mantissa AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
    
        DIM v1 AS DWORD, v2 AS DWORD
        DIM i AS LONG
        FOR i = digits TO 1 STEP -1
            v1 = mantissa.mantissa(i) \ 1000000
            v2 = mantissa.mantissa(i - 1) MOD 1000000
            v2 = v2 * 100 + v1
            mantissa.mantissa(i) = v2
        NEXT
        mantissa.mantissa(0) = mantissa.mantissa(0) \ 1000000
    END SUB
    
    SUB LSHIFT_6 (mantissa AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
    
        DIM v1 AS DWORD, v2 AS DWORD
        DIM i AS LONG
        FOR i = 0 TO digits - 1
            v1 = mantissa.mantissa(i) MOD 100
            v2 = mantissa.mantissa(i + 1) \ 100
            mantissa.mantissa(i) = v1 * 1000000 + v2
            mantissa.mantissa(i + 1) = mantissa.mantissa(i + 1) MOD 100
        NEXT
        mantissa.mantissa(digits) = 1000000 * (mantissa.mantissa(digits) MOD 100)
    END SUB
    
    SUB RSHIFT_7 (mantissa AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
    
        DIM v1 AS DWORD, v2 AS DWORD
        DIM i AS LONG
        FOR i = digits TO 1 STEP -1
            v1 = mantissa.mantissa(i) \ 10000000
            v2 = mantissa.mantissa(i - 1) MOD 10000000
            v2 = v2 * 10 + v1
            mantissa.mantissa(i) = v2
        NEXT
        mantissa.mantissa(0) = mantissa.mantissa(0) \ 10000000
    END SUB
    
    SUB LSHIFT_7 (mantissa AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
    
        DIM v1 AS DWORD, v2 AS DWORD
        DIM i AS LONG
        FOR i = 0 TO digits - 1
            v1 = mantissa.mantissa(i) MOD 10
            v2 = mantissa.mantissa(i + 1) \ 10
            mantissa.mantissa(i) = v1 * 10000000 + v2
            mantissa.mantissa(i + 1) = mantissa.mantissa(i + 1) MOD 10
        NEXT
        mantissa.mantissa(digits) = 10000000 * (mantissa.mantissa(digits) MOD 10)
    END SUB
    
    SUB RSHIFT_8 (mantissa AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM i AS LONG
        FOR i = digits TO 1 STEP -1
            mantissa.mantissa(i) = mantissa.mantissa(i - 1)
        NEXT
        mantissa.mantissa(0) = 0
    END SUB
    
    SUB LSHIFT_8 (mantissa AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM i AS LONG
        FOR i = 0 TO digits - 1
            mantissa.mantissa(i) = mantissa.mantissa(i + 1)
        NEXT
        mantissa.mantissa(digits) = 0
    END SUB
    
    FUNCTION fpcmp& (x AS decfloat, y AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM c AS LONG, i AS LONG
        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 = x.mantissa(i) - 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 sit meaning that x is negative then x < y
    
    END FUNCTION
    
    SUB NORM_FAC1 (fac1 AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        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, sit the exponent and sign equal to 0.
        DIM i AS LONG, er AS LONG, f AS LONG
        er = 0: f = 0
        FOR i = 0 TO digits
            IF 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 fac1.mantissa(0) > 99999999 THEN
            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 fac1.mantissa(0) = 0
                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 fac1.mantissa(0) < 10 THEN
                LSHIFT_7 fac1, digits
                fac1.exponent = fac1.exponent - 7
            ELSEIF fac1.mantissa(0) < 100 THEN
                LSHIFT_6 fac1, digits
                fac1.exponent = fac1.exponent - 6
            ELSEIF fac1.mantissa(0) < 1000 THEN
                LSHIFT_5 fac1, digits
                fac1.exponent = fac1.exponent - 5
            ELSEIF fac1.mantissa(0) < 10000 THEN
                LSHIFT_4 fac1, digits
                fac1.exponent = fac1.exponent - 4
            ELSEIF fac1.mantissa(0) < 100000 THEN
                LSHIFT_3 fac1, digits
                fac1.exponent = fac1.exponent - 3
            ELSEIF fac1.mantissa(0) < 1000000 THEN
                LSHIFT_2 fac1, digits
                fac1.exponent = fac1.exponent - 2
            ELSEIF fac1.mantissa(0) < 10000000 THEN
                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 digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM v AS LONG, c AS LONG, i AS LONG
        c = 0
        FOR i = digits TO 1 STEP -1
            v = fac2.mantissa(i) + fac1.mantissa(i) + c
            IF v > 99999999 THEN
                v = v - 100000000
                c = 1
            ELSE
                c = 0
            END IF
            fac1.mantissa(i) = v
        NEXT
        v = fac1.mantissa(0) + fac2.mantissa(0) + c
        fac1.mantissa(0) = v
    
        NORM_FAC1 fac1, digits
    END SUB
    
    SUB fpsub_aux (fac1 AS decfloat, fac2 AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM v AS LONG, c AS LONG, i AS LONG
        c = 0
        FOR i = digits TO 1 STEP -1
            v = fac1.mantissa(i) - fac2.mantissa(i) - c
            IF v < 0 THEN
                v = v + 100000000
                c = 1
            ELSE
                c = 0
            END IF
            fac1.mantissa(i) = v
        NEXT
        v = fac1.mantissa(0) - fac2.mantissa(0) - c
        fac1.mantissa(0) = v
    
        NORM_FAC1 fac1, digits
    END SUB
    
    SUB fpadd (result AS decfloat, x AS decfloat, y AS decfloat, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
    
        DIM fac1 AS decfloat, fac2 AS decfloat
        DIM i AS LONG, t AS LONG, c AS LONG, xsign AS LONG, ysign AS LONG
    
        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
                    RSHIFT_8 fac2, digits
                    t = t - 8
                    i = i - 1
                WEND
                IF t = 7 THEN
                    RSHIFT_7 fac2, digits
                ELSEIF t = 6 THEN
                    RSHIFT_6 fac2, digits
                ELSEIF t = 5 THEN
                    RSHIFT_5 fac2, digits
                ELSEIF t = 4 THEN
                    RSHIFT_4 fac2, digits
                ELSEIF t = 3 THEN
                    RSHIFT_3 fac2, digits
                ELSEIF t = 2 THEN
                    RSHIFT_2 fac2, digits
                ELSEIF t = 1 THEN
                    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
                fpadd_aux fac1, fac2, digits
            ELSE
                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 digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM fac1 AS decfloat, fac2 AS decfloat
        fac1 = x
        fac2 = y
        fac2.sign = fac2.sign XOR &H8000
        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 digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM fac1 AS decfloat, fac2 AS decfloat
        DIM i AS LONG, j AS LONG, ex AS LONG, er AS LONG, den AS LONG, num AS LONG
        DIM digit AS QUAD, carry AS QUAD
        DIM fac3(0 TO digits) AS QUAD
        DIM ndw AS LONG
        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
                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 fac2.mantissa(den) = 0
                den = den - 1
            WEND
            num = digits
            WHILE 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 = fac2.mantissa(j)
                FOR i = num TO 0 STEP -1
                    fac3(i) = fac3(i) + digit * 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
                ndw = fac3(i)
                fac1.mantissa(i) = ndw
            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
        NORM_FAC1 fac1, digits
        result = fac1
    END SUB
    
    SUB fpmul_si (result AS decfloat, x AS decfloat, y AS QUAD, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM fac1 AS decfloat, fac2 AS decfloat
        DIM count AS LONG, ex AS LONG, er AS LONG, i AS LONG
        DIM carry AS QUAD, digit AS QUAD, prod AS QUAD, value AS QUAD
        DIM ndw AS LONG
        fac1 = x
        digit = ABS(y)
        IF digit > 99999999 THEN
            si2fp fac2, y, digits
            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
                fac1.mantissa(count) = 0
            NEXT
            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 * fac1.mantissa(i) + carry
                value = (prod MOD 100000000)
                ndw = value
                fac1.mantissa(i) = ndw
                carry = prod \ 100000000
            NEXT
    
            IF carry < 10 THEN
                RSHIFT_1 fac1, digits
                fac1.exponent = fac1.exponent + 1
                fac1.mantissa(0) = fac1.mantissa(0) + carry * 10000000
            ELSEIF carry < 100 THEN
                RSHIFT_2 fac1, digits
                fac1.exponent = fac1.exponent + 2
                fac1.mantissa(0) = fac1.mantissa(0) + carry * 1000000
            ELSEIF carry < 1000 THEN
                RSHIFT_3 fac1, digits
                fac1.exponent = fac1.exponent + 3
                fac1.mantissa(0) = fac1.mantissa(0) + carry * 100000
            ELSEIF carry < 10000 THEN
                RSHIFT_4 fac1, digits
                fac1.exponent = fac1.exponent + 4
                fac1.mantissa(0) = fac1.mantissa(0) + carry * 10000
            ELSEIF carry < 100000 THEN
                RSHIFT_5 fac1, digits
                fac1.exponent = fac1.exponent + 5
                fac1.mantissa(0) = fac1.mantissa(0) + carry * 1000
            ELSEIF carry < 1000000 THEN
                RSHIFT_6 fac1, digits
                fac1.exponent = fac1.exponent + 6
                fac1.mantissa(0) = fac1.mantissa(0) + carry * 100
            ELSEIF carry < 10000000 THEN
                RSHIFT_7 fac1, digits
                fac1.exponent = fac1.exponent + 7
                fac1.mantissa(0) = fac1.mantissa(0) + carry * 10
            ELSEIF carry < 100000000 THEN
                RSHIFT_8 fac1, digits
                fac1.exponent = fac1.exponent + 8
                fac1.mantissa(0) = fac1.mantissa(0) + carry
            END IF
    
        END IF
    
        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 QUAD, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM fac1 AS decfloat, fac2 AS decfloat
        DIM count AS LONG, ex AS LONG, er AS LONG, i AS LONG
        DIM carry AS QUAD, digit AS QUAD, prod AS QUAD, value AS QUAD
        DIM n0 AS QUAD, n1 AS QUAD
        DIM ndw AS LONG
        fac1 = x
        digit = ABS(y)
        IF digit > 99999999 AND digit < 10000000000000000 THEN
            n0 = digit \ 100000000
            n1 = digit MOD 100000000
            fpmul_si fac2, fac1, n1, digits
            fac1.exponent = fac1.exponent + 8
            fpmul_si fac1, fac1, n0, digits
            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
            str2fp fac2, STR$(y)
            fpmul fac1, fac1, fac2, digits
            result = fac1: EXIT SUB
        END IF
        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
                fac1.mantissa(count) = 0
            NEXT
            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 * fac1.mantissa(i) + carry
                value = (prod MOD 100000000)
                ndw = value
                fac1.mantissa(i) = ndw
                carry = prod \ 100000000
            NEXT
    
            IF carry < 10 THEN
                RSHIFT_1 fac1, digits
                fac1.exponent = fac1.exponent + 1
                fac1.mantissa(0) = fac1.mantissa(0) + carry * 10000000
            ELSEIF carry < 100 THEN
                RSHIFT_2 fac1, digits
                fac1.exponent = fac1.exponent + 2
                fac1.mantissa(0) = fac1.mantissa(0) + carry * 1000000
            ELSEIF carry < 1000 THEN
                RSHIFT_3 fac1, digits
                fac1.exponent = fac1.exponent + 3
                fac1.mantissa(0) = fac1.mantissa(0) + carry * 100000
            ELSEIF carry < 10000 THEN
                RSHIFT_4 fac1, digits
                fac1.exponent = fac1.exponent + 4
                fac1.mantissa(0) = fac1.mantissa(0) + carry * 10000
            ELSEIF carry < 100000 THEN
                RSHIFT_5 fac1, digits
                fac1.exponent = fac1.exponent + 5
                fac1.mantissa(0) = fac1.mantissa(0) + carry * 1000
            ELSEIF carry < 1000000 THEN
                RSHIFT_6 fac1, digits
                fac1.exponent = fac1.exponent + 6
                fac1.mantissa(0) = fac1.mantissa(0) + carry * 100
            ELSEIF carry < 10000000 THEN
                RSHIFT_7 fac1, digits
                fac1.exponent = fac1.exponent + 7
                fac1.mantissa(0) = fac1.mantissa(0) + carry * 10
            ELSEIF carry < 100000000 THEN
                RSHIFT_8 fac1, digits
                fac1.exponent = fac1.exponent + 8
                fac1.mantissa(0) = fac1.mantissa(0) + carry
            END IF
    
        END IF
    
        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 digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM x AS DOUBLE
        DIM k AS LONG, l AS LONG, ex AS LONG
        DIM prec AS LONG
        DIM r AS decfloat, r2 AS decfloat, one AS decfloat, n AS decfloat
        n = m
        l = LOG((%NUM_DIGITS + 8) * 0.0625) * 1.5
        DIM v AS STRING, ts AS STRING
        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$(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$(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
            str2fp r, "1"
            result = r: EXIT SUB
        ELSEIF x = 1 THEN
            x = 10
            ex = ex - 1
        END IF
        ex = (-1) - ex
        x = 1 / x
        str2fp r, STR$(x)
        r.exponent = ex + %BIAS + 1
        one.mantissa(0) = 10000000
        one.exponent = %BIAS + 1
        r2 = r
        prec = 3
        FOR k = 1 TO l
            prec = 2 * prec - 1
            fpmul r2, n, r, prec
            fpsub r2, one, r2, prec
            fpmul r2, r, r2, prec
            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 j AS LONG
        FOR j = ka TO kb
            w(j) = w(j) - q * d(j - ka + 2)
        NEXT
    END SUB
    
    SUB normalise (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 carry AS LONG, j AS LONG
        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 digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM fac1 AS decfloat, fac2 AS decfloat
        DIM i AS LONG, er AS LONG, is_power_of_ten AS LONG
    
        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
                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 fac2.mantissa(0) = 10000000 THEN
                is_power_of_ten = 1
                FOR i = 1 TO digits
                    IF 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 result(1 TO 2 * digits + 3) AS DOUBLE
            DIM n(1 TO 2 * digits + 3) AS DOUBLE
            DIM d(1 TO 2 * digits + 3) AS DOUBLE
            DIM b AS LONG : b= 10000
            DIM j AS LONG, last AS LONG, laststep AS LONG, q AS LONG, t AS LONG
            DIM stp AS LONG
            DIM xd AS DOUBLE, xn AS DOUBLE, rund AS DOUBLE
            DIM w(1 TO UBOUND(n) + 4) AS DOUBLE
    
            FOR j = 0 TO digits
                n(2 * j + 2) = fac1.mantissa(j) \ 10000
                n(2 * j + 3) = fac1.mantissa(j) MOD 10000
                d(2 * j + 2) = fac2.mantissa(j) \ 10000
                d(2 * j + 3) = 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))
                subtract w(), q, d(), (stp + 2), last
                normalise w(), (stp + 2), q
            NEXT
            finalnorm w(), (laststep + 1)
            IF w(2) <> 0 THEN laststep = laststep - 1
            rund = w(laststep + 1) / b
            IF rund >= 0.5 THEN w(laststep) = w(laststep) + 1
            IF w(2) = 0 THEN
                FOR j = 1 TO t + 1
                    result(j) = w(j + 1)
                NEXT
            ELSE
                FOR j = 1 TO t + 1
                    result(j) = w(j)
                NEXT
            END IF
            IF w(2) = 0 THEN result(1) = w(1) - 1 ELSE result(1) = w(1)
    
            FOR j = 0 TO digits
                fac1.mantissa(j) = result(2 * j + 2) * 10000 + result(2 * j + 3)
            NEXT
            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 digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM r AS decfloat, r2 AS decfloat, tmp AS decfloat
        DIM n AS decfloat, half AS decfloat
        DIM ex AS LONG, k AS LONG, l AS LONG, prec AS LONG
        DIM ts AS STRING, v AS STRING
        DIM x AS DOUBLE
        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
        si2fp tmp, 0, digits
        IF fpcmp(n, tmp, digits) = 0 THEN
            si2fp r, 0, digits
            result = r
            EXIT SUB
        END IF
        si2fp tmp, 1, digits
        IF fpcmp(n, tmp, digits) = 0 THEN
            si2fp r, 1, digits
            result = r
            EXIT SUB
        END IF
        si2fp tmp, 0, digits
        IF fpcmp(n, tmp, digits) < 0 THEN
            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$(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$(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
            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, ".")
        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
            half.mantissa(k) = 0
        NEXT
        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
            fpdiv tmp, n, r, prec
            fpadd r2, r, tmp, prec
            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 digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM fac1 AS decfloat
        DIM carry AS QUAD, remder AS QUAD
        DIM i AS QUAD, divisor AS QUAD
        DIM quotient AS QUAD
        DIM ndw AS LONG
        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
            ndw = i
            quotient = fac1.mantissa(ndw) + remder * 100000000
            remder = quotient MOD divisor
            fac1.mantissa(ndw) = quotient \ divisor
        NEXT
        quotient = remder * 100000000
        quotient = quotient \ divisor
        carry = fac1.mantissa(0)
    
        IF carry = 0 THEN
            LSHIFT_8 fac1, digits
            fac1.exponent = fac1.exponent - 8
            fac1.mantissa(digits) = fac1.mantissa(digits) + quotient
        ELSEIF carry < 10 THEN
            LSHIFT_7 fac1, digits
            fac1.exponent = fac1.exponent - 7
            fac1.mantissa(digits) = fac1.mantissa(digits) + quotient \ 10
        ELSEIF carry < 100 THEN
            LSHIFT_6 fac1, digits
            fac1.exponent = fac1.exponent - 6
            fac1.mantissa(digits) = fac1.mantissa(digits) + quotient \ 100
        ELSEIF carry < 1000 THEN
            LSHIFT_5 fac1, digits
            fac1.exponent = fac1.exponent - 5
            fac1.mantissa(digits) = fac1.mantissa(digits) + quotient \ 1000
        ELSEIF carry < 10000 THEN
            LSHIFT_4 fac1, digits
            fac1.exponent = fac1.exponent - 4
            fac1.mantissa(digits) = fac1.mantissa(digits) + quotient \ 10000
        ELSEIF carry < 100000 THEN
            LSHIFT_3 fac1, digits
            fac1.exponent = fac1.exponent - 3
            fac1.mantissa(digits) = fac1.mantissa(digits) + quotient \ 100000
        ELSEIF carry < 1000000 THEN
            LSHIFT_2 fac1, digits
            fac1.exponent = fac1.exponent - 2
            fac1.mantissa(digits) = fac1.mantissa(digits) + quotient \ 1000000
        ELSEIF carry < 10000000 THEN
            LSHIFT_1 fac1, digits
            fac1.exponent = fac1.exponent - 1
            fac1.mantissa(digits) = 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 digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM n AS decfloat
        fpfix n, num
        fpsub n, num, n, digits
        result = n
    END SUB
    
    'returns the positive of n
    SUB fpabs (result AS decfloat, n AS decfloat)
        DIM x AS decfloat
        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 x AS decfloat
        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 x AS decfloat
        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 q AS decfloat, fm AS decfloat 'make copy in case the destination and source are the same
        fpdiv fm, num, denom, digits: fpfix q, fm
        fpsub fm, fm, q, digits: 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 ep AS decfloat
        si2fp ep, 1, digits
        ep.exponent = (-(%NUM_DIGITS) + %BIAS + 1)
        result = ep
    END SUB
    
    SUB fpipow (result AS decfloat, x AS decfloat, e AS QUAD, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        'take x to an Long power
        DIM y AS decfloat
        DIM z AS decfloat
        DIM n AS QUAD, c AS QUAD
        DIM i AS LONG
        c = 0
        y = x
        n = ABS(e)
        z.sign = 0
        z.exponent = (%BIAS + 1)
        z.mantissa(0) = 10000000
        FOR i = 1 TO %NUM_DWORDS
            z.mantissa(i) = 0
        NEXT
        WHILE n > 0
            WHILE (n AND 1) = 0
                n = n \ 2
                fpmul y, y, y, digits
                c = c + 1
            WEND
            n = n - 1
            fpmul z, y, z, digits
            c = c + 1
        WEND
        IF e < 0 THEN
            fpinv z, z, digits
        END IF
        result = z
    END SUB
    
    SUB fpfactorial (result AS decfloat, n AS LONG, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        DIM i AS QUAD
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM f AS decfloat
        IF n < 0 THEN PRINT "inf": result = f: EXIT SUB
        IF n = 0 OR n = 1 THEN
            si2fp f, 1, digits
            result = f: EXIT SUB
        END IF
        si2fp f, 2, digits
        IF n = 2 THEN result = f: EXIT SUB
        FOR i = 3 TO n
            fpmul_si f, f, i, digits
        NEXT
        result = f
    END SUB
    
    SUB fpnroot (result AS decfloat, x AS decfloat, p AS LONG, digits_in AS LONG)
        DIM digits AS LONG
        digits = digits_in
        IF digits > %NUM_DWORDS THEN digits = %NUM_DWORDS
        DIM ry AS decfloat, tmp AS decfloat, tmp2 AS decfloat
        DIM t AS DOUBLE, t2 AS DOUBLE
        DIM i AS LONG, ex AS LONG, l AS LONG, prec AS LONG
    
        l = LOG((%NUM_DIGITS + 8) * 0.0625) * 1.5
        ex = (x.exponent AND &H7FFFFFFF) - %BIAS - 1
        t = x.mantissa(0) + x.mantissa(1) / 100000000
        t = t / 10000000
        t2 = LOG(t) / p + LOG(10) * ex / p
        t2 = EXP(t2)
        str2fp ry, STR$(t2)
        prec = 3
    
        fpipow tmp, ry, p - 1, prec
        fpdiv tmp, x, tmp, prec
        fpmul_si tmp2, ry, p - 1, prec
        fpadd tmp2, tmp2, tmp, prec
        fpdiv_si ry, tmp2, p, prec
        FOR i = 1 TO l
            prec = 2 * prec - 1
            fpipow tmp, ry, p - 1, prec
            fpdiv tmp, x, tmp, prec
            fpmul_si tmp2, ry, p - 1, prec
            fpadd tmp2, tmp2, tmp, prec
            fpdiv_si ry, tmp2, p, prec
        NEXT
        result = ry
    END SUB
    
    SUB pi_brent_salamin (pi_bs AS decfloat, digits_in AS DWORD)
        DIM digits AS DWORD
        digits = digits_in
        IF digits > %NUMBER_OF_DIGITS THEN digits = %NUMBER_OF_DIGITS
        DIM i AS LONG
        DIM limit2 AS LONG
        DIM c0 AS decfloat, c1 AS decfloat, c2 AS decfloat, c05 AS decfloat
        DIM a AS decfloat, b AS decfloat, sum AS decfloat
        DIM ak AS decfloat, bk AS decfloat, ck AS decfloat
        DIM ab AS decfloat, asq AS decfloat
        DIM pow2 AS decfloat, tmp AS decfloat
    
        limit2 = -digits + %BIAS + 1
        si2fp c0, 0, %NUMBER_OF_DIGITS: ak = c0: bk = c0: ab = c0: asq = c0
        si2fp c1, 1, %NUMBER_OF_DIGITS: a = c1: ck = c1: pow2 = c1
        si2fp c2, 2, %NUMBER_OF_DIGITS: b = c2
        str2fp c05, ".5": sum = c05
        si2fp pi_bs, 3, %NUMBER_OF_DIGITS
    
        fpsqr b, b, %NUMBER_OF_DIGITS
        fpdiv b, c1, b, %NUMBER_OF_DIGITS
        WHILE fpcmp(ck, c0, %NUMBER_OF_DIGITS) <> 0 AND ck.exponent > limit2
            fpadd ak, a, b, %NUMBER_OF_DIGITS
            fpmul ak, c05, ak, %NUMBER_OF_DIGITS
            fpmul ab, a, b, %NUMBER_OF_DIGITS
            fpsqr bk, ab, %NUMBER_OF_DIGITS
            fpmul asq, ak, ak, %NUMBER_OF_DIGITS
            fpsub ck, asq, ab, %NUMBER_OF_DIGITS
            fpmul pow2, pow2, c2, %NUMBER_OF_DIGITS
            fpmul tmp, pow2, ck, %NUMBER_OF_DIGITS
            fpsub sum, sum, tmp, %NUMBER_OF_DIGITS
            a = ak: b = bk
        WEND
        fpdiv tmp, asq, sum, %NUMBER_OF_DIGITS
        fpmul pi_bs, c2, tmp, %NUMBER_OF_DIGITS
    END SUB
    as the salmon fish is compelled to go back to it's birthplace to spawn, so comes a day when man is compelled to go back to it's source.. GOD

    Comment


    • #22
      in the sub pi_brent_salamin all goes well until the second round in the while loop, the subtraction is wrong but the same numbers subtracted in pbmain are ok, what's going on?
      Code:
      #COMPILE EXE
      #DIM ALL
      
      %NUMBER_OF_DIGITS = 100
      %NUM_DIGITS = %NUMBER_OF_DIGITS
      %NUM_DWORDS = %NUM_DIGITS \ 8
      %BIAS = 1073741824 '2 ^ 30
      %MANTISSA_BYTES = (%NUM_DWORDS + 1) * 8 * 4
      
      TYPE decfloat
          sign AS LONG
          exponent AS DWORD
          mantissa(%NUM_DWORDS) AS DWORD
      END TYPE
      
      ' Error definitions
      
      %DIVZ_ERR = 1 'Divide by zero
      %EXPO_ERR = 2 'Exponent overflow error
      
      FUNCTION PBMAIN () AS LONG
      
      DIM x AS decfloat, y AS decfloat, z AS decfloat
      DIM s AS STRING
      pi_brent_salamin z, 100
      str2fp x, "0.717790036133724090747785501167681602903443553780351068051309250456926861978485274721876331803523"
      str2fp y, "0.717749986377537538194937730756726426343260154294346608661363599167329044533680409874857338234180"
      fpsub z, x, y, %NUMBER_OF_DIGITS
      PRINT fp2str(z, %NUMBER_OF_DIGITS)
      WAITKEY$
      
      END FUNCTION
      
      SUB fpcopy(result AS decfloat, num AS decfloat)
          DIM i AS LONG
          result.sign=num.sign
          result.exponent=num.exponent
          FOR i=0 TO %NUM_DWORDS
              result.mantissa(i)=num.mantissa(i)
          NEXT
      END SUB
      and pi_brent_salamin
      Code:
      SUB pi_brent_salamin (pi_bs AS decfloat, digits_in AS DWORD)
          DIM digits AS DWORD
          digits = digits_in
          IF digits > %NUMBER_OF_DIGITS THEN digits = %NUMBER_OF_DIGITS
          DIM i AS LONG
          DIM limit2 AS LONG
          DIM c0 AS decfloat, c1 AS decfloat, c2 AS decfloat, c05 AS decfloat
          DIM a AS decfloat, b AS decfloat, sum AS decfloat
          DIM ak AS decfloat, bk AS decfloat, ck AS decfloat
          DIM ab AS decfloat, asq AS decfloat
          DIM pow2 AS decfloat, tmp AS decfloat
      
          limit2 = -digits + %BIAS + 1
          si2fp c0, 0, %NUMBER_OF_DIGITS: ak = c0: bk = c0: ab = c0: asq = c0
          si2fp c1, 1, %NUMBER_OF_DIGITS: a = c1: ck = c1: pow2 = c1
          si2fp c2, 2, %NUMBER_OF_DIGITS: b = c2
          str2fp c05, ".5": sum = c05
          si2fp pi_bs, 3, %NUMBER_OF_DIGITS
      
          fpsqr b, b, %NUMBER_OF_DIGITS
          fpdiv b, c1, b, %NUMBER_OF_DIGITS
          'WHILE fpcmp(ck, c0, %NUMBER_OF_DIGITS) <> 0 AND ck.exponent > limit2
          FOR i=1 TO 2
              fpadd ak, a, b, %NUMBER_OF_DIGITS
              fpmul ak, c05, ak, %NUMBER_OF_DIGITS
              fpmul ab, a, b, %NUMBER_OF_DIGITS
              fpsqr bk, ab, %NUMBER_OF_DIGITS
              fpmul asq, ak, ak, %NUMBER_OF_DIGITS
              fpsub ck, asq, ab, %NUMBER_OF_DIGITS
              PRINT fp2str(asq, %NUMBER_OF_DIGITS)
              PRINT fp2str(ab, %NUMBER_OF_DIGITS)
              PRINT fp2str(ck, %NUMBER_OF_DIGITS)
              ? "-------------------------------------------------------------------"
              fpmul pow2, pow2, c2, %NUMBER_OF_DIGITS
              fpmul tmp, pow2, ck, %NUMBER_OF_DIGITS
              fpsub sum, sum, tmp, %NUMBER_OF_DIGITS
              '? fp2str(sum, %NUMBER_OF_DIGITS)
              fpcopy a, ak: fpcopy b, bk
          NEXT
          'WEND
          fpdiv tmp, asq, sum, %NUMBER_OF_DIGITS
          fpmul pi_bs, c2, tmp, %NUMBER_OF_DIGITS
      END SUB
      Last edited by Johan Klassen; 15 May 2022, 06:57 AM.
      as the salmon fish is compelled to go back to it's birthplace to spawn, so comes a day when man is compelled to go back to it's source.. GOD

      Comment


      • #23
        Johann this is great i really appreciate the chance to look at your code. I was able to do some testing. Your code is about 8 times faster than mine, which makes sense as you are taking the data 8 digits at a time vs my one digit at a time. For example on my machine, to do 1000 multiplications of 5000 digit numbers takes 95 seconds, your code on my machine for the exact same calculation took 11 seconds, about 8x faster. The factor of 8 held up for different numbers of decimals.

        Comment


        • #24
          just to remind you that something is not right, the pi_brent_salamin sub doesn't work, and I give up trying to figure out what the problem might be, it works without a problem in other basic's
          as the salmon fish is compelled to go back to it's birthplace to spawn, so comes a day when man is compelled to go back to it's source.. GOD

          Comment

          Working...
          X