' DOUBLE Precision Math — pi, e, sqrt convergence
#COMPILE EXE
#DIM ALL

FUNCTION ComputePiLeibniz(terms AS LONG) AS DOUBLE
    LOCAL i AS LONG
    LOCAL pi AS DOUBLE
    LOCAL denom AS DOUBLE
    pi = 0.0
    FOR i = 0 TO terms - 1
        denom = 2.0 * i + 1.0
        IF i MOD 2 = 0 THEN
            pi = pi + 1.0 / denom
        ELSE
            pi = pi - 1.0 / denom
        END IF
    NEXT i
    FUNCTION = 4.0 * pi
END FUNCTION

FUNCTION ComputePiNilakantha(iterations AS LONG) AS DOUBLE
    LOCAL i AS LONG
    LOCAL pi AS DOUBLE
    LOCAL n AS DOUBLE
    pi = 3.0
    FOR i = 1 TO iterations
        n = 2.0 * i
        IF i MOD 2 = 1 THEN
            pi = pi + 4.0 / (n * (n + 1.0) * (n + 2.0))
        ELSE
            pi = pi - 4.0 / (n * (n + 1.0) * (n + 2.0))
        END IF
    NEXT i
    FUNCTION = pi
END FUNCTION

FUNCTION ComputeE(terms AS LONG) AS DOUBLE
    LOCAL i AS LONG
    LOCAL e AS DOUBLE
    LOCAL fact AS DOUBLE
    e = 1.0
    fact = 1.0
    FOR i = 1 TO terms
        fact = fact * i
        e = e + 1.0 / fact
    NEXT i
    FUNCTION = e
END FUNCTION

FUNCTION SqrtNewton(x AS DOUBLE, iterations AS LONG) AS DOUBLE
    LOCAL i AS LONG
    LOCAL guess AS DOUBLE
    IF x <= 0 THEN FUNCTION = 0
    guess = x / 2.0
    FOR i = 1 TO iterations
        guess = (guess + x / guess) / 2.0
    NEXT i
    FUNCTION = guess
END FUNCTION

SUB ConvergenceDemo()
    LOCAL i AS LONG
    LOCAL piLeib AS DOUBLE
    LOCAL piNila AS DOUBLE
    LOCAL eVal AS DOUBLE
    LOCAL sqrtVal AS DOUBLE
    PRINT "=== CONVERGENCE DEMONSTRATION ==="
    PRINT
    PRINT "Pi by Leibniz series:"
    PRINT "Terms    Value               Error"
    PRINT STRING$(50, 61)
    FOR i = 1 TO 10
        piLeib = ComputePiLeibniz(i * 100)
        PRINT FORMAT$(i * 100, "#####") & "   " & FORMAT$(piLeib, "##.##############") & "   " & FORMAT$(piLeib - 3.14159265358979, "##.##############")
    NEXT i
    PRINT "Reference: 3.14159265358979"
    PRINT
    PRINT "Pi by Nilakantha series:"
    PRINT "Terms    Value               Error"
    PRINT STRING$(50, 61)
    FOR i = 1 TO 10
        piNila = ComputePiNilakantha(i * 50)
        PRINT FORMAT$(i * 50, "#####") & "   " & FORMAT$(piNila, "##.##############") & "   " & FORMAT$(piNila - 3.14159265358979, "##.##############")
    NEXT i
    PRINT
    PRINT "Euler's number e:"
    PRINT "Terms    Value               Error"
    PRINT STRING$(50, 61)
    FOR i = 1 TO 10
        eVal = ComputeE(i * 5)
        PRINT FORMAT$(i * 5, "#####") & "   " & FORMAT$(eVal, "##.##############") & "   " & FORMAT$(eVal - 2.71828182845905, "##.##############")
    NEXT i
    PRINT
    PRINT "Newton sqrt convergence for x=12345.6789:"
    PRINT "Iter     Value               Error"
    PRINT STRING$(50, 61)
    FOR i = 1 TO 8
        sqrtVal = SqrtNewton(12345.6789, i)
        PRINT FORMAT$(i, "####") & "   " & FORMAT$(sqrtVal, "###.##############") & "   " & FORMAT$(sqrtVal - SQR(12345.6789), "##.##############")
    NEXT i
    PRINT
END SUB

SUB ExtendedPrecisionDemo()
    LOCAL a AS DOUBLE
    LOCAL b AS DOUBLE
    LOCAL c AS DOUBLE
    a = 1.234567890123456789
    b = 9.876543210987654321
    c = a * b
    PRINT "=== DOUBLE PRECISION ==="
    PRINT "a = " & FORMAT$(a, "#.################")
    PRINT "b = " & FORMAT$(b, "#.################")
    PRINT "a * b = " & FORMAT$(c, "###.################")
    PRINT "a / b = " & FORMAT$(a / b, "#.################")
    PRINT "SIN(a) = " & FORMAT$(SIN(a), "#.################")
    PRINT "COS(b) = " & FORMAT$(COS(b), "#.################")
    PRINT "SIZEOF(EXT) = " & STR$(SIZEOF(a))
    PRINT
END SUB

SUB SeriesSums()
    LOCAL n AS LONG
    LOCAL harmonic AS DOUBLE
    LOCAL altHarmonic AS DOUBLE
    LOCAL geom AS DOUBLE
    LOCAL i AS LONG
    harmonic = 0.0
    altHarmonic = 0.0
    geom = 0.0
    PRINT "=== SERIES SUMS ==="
    n = 10000
    FOR i = 1 TO n
        harmonic = harmonic + 1.0 / i
        IF i MOD 2 = 1 THEN
            altHarmonic = altHarmonic + 1.0 / i
        ELSE
            altHarmonic = altHarmonic - 1.0 / i
        END IF
    NEXT i
    LOCAL pow2 AS DOUBLE
    pow2 = 1.0
    FOR i = 0 TO 20
        geom = geom + 1.0 / pow2
        pow2 = pow2 * 2.0
    NEXT i
    PRINT "Harmonic sum (n=" & STR$(n) & "): " & FORMAT$(harmonic, "##.########")
    PRINT "ln(" & STR$(n) & ") + gamma ~ " & FORMAT$(LOG(n) + 0.5772156649, "##.########")
    PRINT "Alternating harmonic: " & FORMAT$(altHarmonic, "##.########")
    PRINT "ln(2) = " & FORMAT$(LOG(2.0), "##.########")
    PRINT "Geometric series (1/2^n): " & FORMAT$(geom, "##.########")
    PRINT "Expected: 2.0"
    PRINT
END SUB

FUNCTION PBMAIN() AS LONG
    PRINT "=================================="
    PRINT " DOUBLE Precision Math Library"
    PRINT "=================================="
    PRINT
    CALL ExtendedPrecisionDemo()
    CALL ConvergenceDemo()
    CALL SeriesSums()
    PRINT "DOUBLE precision demos complete."
    FUNCTION = 0
END FUNCTION
