Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
KISS MY ASCII GOOD PI!
#11
out of curiosity I calculated Pi to 1,000,000 digits using my math routines, it took 16.5 minutes, about 330 times longer than your routine Steve Big Grin
Reply
#12
(08-22-2022, 11:44 PM)Jack Wrote: out of curiosity I calculated Pi to 1,000,000 digits using my math routines, it took 16.5 minutes, about 330 times longer than your routine Steve Big Grin

You just need a bigger seed to quickly get you started on your way faster.  Wink
Reply
#13
(08-22-2022, 10:56 PM)SMcNeill Wrote: Here's my little attempt at finding Pi.  This relies on some very impressive math formulas and such, but it can calculate up to a million digits of PI for us in about 1/10th of a second!  (Just expect QB64 to spend quite a bit longer to print the results to screen for you than that -- our print routine isn't very quick at all!)


[Image: image.png]

The program here is a little too large to fit and play nicely inside a code box, so I'll attached it below for ease of download and reference.  

I hope you appreciate all the time and effort I put into this for you @Pete!  I doubt you'll ever find another method any more accurate or faster than this one!!  Tongue

Yeah, that code deserves a pi in the face! Lemon comes to mind. Rolleyes 

After doing some reading, I found Ramanujan's method to be interesting, but to program it I'd have to go back in time, 50 years, and re-learn how to do factorials in linear equations. Of course, calculating it to just 6 decimal places is a snap, as at n=zero, all the factorials drop out! (exclamation point after out is a pun at factorials.) Iterations past that require two factorial computations per iteration.

Other interesting things I came across online was the other guys pi formula I experimented with in my string math routine works, but get this, it apparently takes 5 million iterations to do what Ramanujan's method does in just one; and apparently pi to six places is good enough to measure the globe within 1-meter of accuracy. Bad news, @SMcNeill . That 3-feet is coming out of your farm!!! Consider it cropped.

Pete Big Grin
Reply
#14
Ramanujan's formula gives about 8 digits per iteration whereas the Gauss–Legendre algorithm doubles the accurate digits per iteration, the first iteration gives 2 digits then it doubles on subsequent iterations, for high precision calculation it's much faster
Reply
#15
just for you Pete, here's the Ramanujan formula simplified Cool
Code: (Select All)
Function Ramanujan# ()
    Dim As Double sum, f, f4, f4k, c1, c2, c3, c4, c34k
    Dim As Long k, k4

    c1 = 1103
    c2 = 26390
    c3 = 396
    f = 1
    f4k = 1
    sum = 1103
    c34k = 1
    k4 = 0
    c4 = c3 * c3: c4 = c4 * c4
    For k = 1 To 2
        f = f * k
        f4 = f * f: f4 = f4 * f4
        c34k = c34k * c4
        f4k = f4k * (k4 + 1) * (k4 + 2) * (k4 + 3) * (k4 + 4): k4 = k4 + 4
        sum = sum + (f4k * (c1 + c2 * k)) / (f4 * c34k)
    Next
    Ramanujan = 1 / (2 * Sqr(2) / 9801 * sum)
End Function
to calculate Pi to 1 million digits using Ramanujan would take 125000 iterations, it would take 20 iterations with the Gauss–Legendre algorithm
Reply
#16
(08-24-2022, 02:20 AM)Jack Wrote: just for you Pete, here's the Ramanujan formula simplified Cool
Code: (Select All)
Function Ramanujan# ()
    Dim As Double sum, f, f4, f4k, c1, c2, c3, c4, c34k
    Dim As Long k, k4

    c1 = 1103
    c2 = 26390
    c3 = 396
    f = 1
    f4k = 1
    sum = 1103
    c34k = 1
    k4 = 0
    c4 = c3 * c3: c4 = c4 * c4
    For k = 1 To 2
        f = f * k
        f4 = f * f: f4 = f4 * f4
        c34k = c34k * c4
        f4k = f4k * (k4 + 1) * (k4 + 2) * (k4 + 3) * (k4 + 4): k4 = k4 + 4
        sum = sum + (f4k * (c1 + c2 * k)) / (f4 * c34k)
    Next
    Ramanujan = 1 / (2 * Sqr(2) / 9801 * sum)
End Function
to calculate Pi to 1 million digits using Ramanujan would take 125000 iterations, it would take 20 iterations with the Gauss–Legendre algorithm

@Jack

Thanks a ton for posting this. I'm not sure how long it would take me to re-acquaint myself the factoring when applied to linear equations.

Something I found interesting when converting it to string math was the precision lacking in numeric computer math balanced out the results between the square root calculation and the other larger digit division operations. For fun, I included a way to make the string square root the same digits as the numeric one. If you un-remark that '===========================> line, you will see what I mean. The numeric and string pi numbers will no longer match to the non-greyed out decimal places.

Out of curiosity, is there a way this routine can be used to calculate pi to the next 8 digits and so on? s when in the formula, n = 0, n = 1, n = 2, etc. Increasing the k loops just doesn't seem to produce that effect tot he correct output. For instance, Ram's 22-digit 2 iteration value as per an online calculator is reported as: 3.141592653589793238462

Code: (Select All)
WIDTH 160, 42
_SCREENMOVE 0, 0
DIM SHARED Ramjan$, limit&&, beta
limit&& = 16
'beta = -1
LOCATE 1, 1: PRINT "Jack's Numeric Results:    ";
IF beta THEN LOCATE 5: PRINT
PRINT Ramanujan#
PRINT
LOCATE 3, 1: PRINT "Unrounded String Math Results: "; MID$(Ramjan$, 1, 17);: COLOR 8, 0: PRINT MID$(Ramjan$, 18)
COLOR 7, 0
END

FUNCTION Ramanujan# ()
    DIM AS DOUBLE sum, f, f4, f4k, c1, c2, c3, c34k
    DIM AS LONG k, k4

    sqrt$ = "": CALL square_root("8", sqrt$, limit&&)
    'sqrt$ = MID$(sqrt$, 1, 8) '=============================>

    stringmatha$ = sqrt$
    stringmathb$ = "9801"
    operator$ = "/"
    CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
    numerator$ = runningtotal$

    c1$ = "1103"
    c2$ = "26390"
    c3$ = "396"
    f$ = "1"
    f4k$ = "1"
    sum$ = "1103"
    c34k$ = "1"
    k4$ = "0"

    '----------------------
    c1 = 1103
    c2 = 26390
    c3 = 396
    f = 1
    f4k = 1
    sum = 1103
    c34k = 1
    k4 = 0
    '----------------------

    stringmatha$ = c3$
    FOR i = 1 TO 3
        stringmathb$ = c3$
        operator$ = "*"
        CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
        stringmatha$ = runningtotal$
    NEXT
    c3$ = runningtotal$: c3 = c3 * c3 * c3 * c3
    IF beta THEN PRINT "c3^4 = "; runningtotal$, c3

    stringmatha$ = c3$
    stringmathb$ = "9801"
    operator$ = "/"
    CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
    numerator$ = runningtotal$
    IF beta THEN PRINT "c3^4 / 9801 = "; runningtotal$, c3 / 9801

    FOR k = 1 TO 2
        stringmatha$ = f$
        stringmathb$ = LTRIM$(STR$(k))
        operator$ = "*"
        CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
        f$ = runningtotal$: f = f * (k)

        stringmatha$ = f$
        stringmathb$ = f$
        operator$ = "*"
        CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
        f4$ = runningtotal$
        stringmathb$ = runningtotal$
        CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
        f4$ = runningtotal$: f4 = f * f: f4 = f4 * f4

        stringmatha$ = c34k$
        stringmathb$ = c3$
        operator$ = "*"
        CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
        c34k$ = runningtotal$: c34k = c34k * c3

        stringmatha$ = f$
        stringmathb$ = LTRIM$(STR$(k))
        operator$ = "*"
        CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
        f$ = runningtotal$: f = f * k

        '-------------------------
        REDIM k4$(4): k4$(0) = k4$
        FOR i = 1 TO 4
            stringmatha$ = k4$(i - 1)
            stringmathb$ = "1"
            operator$ = "+"
            CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
            k4$(i) = runningtotal$
        NEXT

        FOR i = 1 TO 4
            stringmatha$ = f4k$
            stringmathb$ = k4$(i)
            operator$ = "*"
            CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
            f4k$ = runningtotal$
        NEXT
        f4k = f4k * (k4 + 1) * (k4 + 2) * (k4 + 3) * (k4 + 4)

        ' Increase k4$ variable.
        stringmatha$ = k4$
        stringmathb$ = "4"
        operator$ = "+"
        CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
        k4$ = runningtotal$: k4 = k4 + 4

        ' Calculate sum.
        stringmatha$ = c2$
        stringmathb$ = LTRIM$(STR$(k))
        operator$ = "*"
        CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
        stringmatha$ = runningtotal$
        stringmathb$ = c1$
        IF beta THEN PRINT: PRINT "String variables c2$ k$: "; runningtotal$, c1$, "   Numeric variables:"; (c2 * k), c1
        operator$ = "+"
        CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
        stringmatha$ = runningtotal$
        stringmathb$ = f4k$
        operator$ = "*"
        CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
        term1$ = runningtotal$

        stringmatha$ = f4$
        stringmathb$ = c34k$
        operator$ = "*"
        CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
        term2$ = runningtotal$

        stringmatha$ = term1$
        stringmathb$ = term2$
        operator$ = "/"
        CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)

        IF beta THEN
            PRINT
            COLOR 4
            PRINT "term1$ / term2$: "; term1$; " / "; term2$
            PRINT "term1 / term2:  "; ((c2 * k) + c1) * f4k; "/"; c34k * f4
            PRINT "String division = "; runningtotal$, "Numeric division ="; ((c2 * k) + c1) / c34k * f4
            COLOR 7, 0
        END IF

        stringmatha$ = runningtotal$
        stringmathb$ = sum$
        operator$ = "+"
        CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
        sum$ = runningtotal$

        sum = sum + (f4k * (c1 + c2 * (k))) / (f4 * c34k)

        IF beta THEN
            PRINT
            COLOR 2, 0: PRINT "String SQR(8) = "; sqrt$, "   Numeric 2 * SQR(2) ="; 2 * SQR(2)
            PRINT "String sum$ = "; sum$, "    Numeric sum ="; sum
            stringmatha$ = "9801"
            stringmathb$ = sum$
            operator$ = "*"
            CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
            PRINT "String 9801 * sum$ = "; runningtotal$, "Numeric 9801 * sum ="; 9801 * sum
            COLOR 7, 0
        END IF
    NEXT

    stringmatha$ = sqrt$
    stringmathb$ = "9801"
    operator$ = "/"
    CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
    stringmatha$ = runningtotal$
    IF beta THEN COLOR 6, 0: PRINT: PRINT "String sqr 8 / 9801: "; runningtotal$, "   Numeric 2 * SQR(2) / 9801 ="; (2 * SQR(2)) / 9801
    stringmathb$ = sum$
    operator$ = "*"
    CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
    IF beta THEN PRINT "String denominator:"; runningtotal$, "    Numeric denominator: "; (2 * SQR(2) / 9801 * sum)
    stringmatha$ = "1"
    stringmathb$ = runningtotal$
    operator$ = "/"
    CALL string_math(stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
    IF beta THEN PRINT "String pi = "; runningtotal$, "    Numeric pi ="; 1 / (2 * SQR(2) / 9801 * sum): COLOR 7, 0
    Ramjan$ = runningtotal$
    IF beta THEN LOCATE 1, 28
    Ramanujan = 1 / (SQR(8) / 9801 * sum)
END FUNCTION

SUB square_root (x$, sqrt$, limit&&)
    oldy$ = ""

    IF INSTR(x$, ".") THEN
        decx$ = MID$(x$, 1, INSTR(x$, ".") - 1)
        x$ = MID$(x$, 1, INSTR(x$, ".") - 1) + MID$(x$, INSTR(x$, ".") + 1)
        IF LEN(x$) = 1 THEN x$ = x$ + "0"
    ELSE
        decx$ = x$
    END IF

    j&& = LEN(decx$)

    ' VAL() okay, one character eval.
    IF VAL(RIGHT$(LTRIM$(STR$(j&&)), 1)) / 2 = VAL(RIGHT$(LTRIM$(STR$(j&&)), 1)) \ 2 THEN
        i&& = 1 ' Even number length.
    ELSE
        i&& = 0 ' Odd number length.
    END IF

    DO
        stringmatha$ = z$: stringmathb$ = k$

        string_math z$, "-", k$, runningtotal$, terminating_decimal%, limit&&
        z$ = runningtotal$ + (MID$(x$, i&&, 2))
        IF LEFT$(z$, 1) = "0" THEN z$ = MID$(z$, 2) ' Remove leading zeros

        oldy$ = ""
        FOR j&& = 1 TO 10
            IF i&& > 1 THEN
                string_math sqrt$, "*", "2", y$, terminating_decimal%, limit&&
                y$ = y$ + LTRIM$(STR$(j&&))
            ELSE
                y$ = LTRIM$(STR$(j&&))
            END IF

            string_math y$, "*", LTRIM$(STR$(j&&)), runningtotal$, terminating_decimal%, limit&&

            string_compare runningtotal$, z$, gl%
            IF gl% > -1 THEN
                IF gl% = 0 THEN
                    h% = 0: oldy$ = y$ ' Perfect square division.
                ELSE
                    h% = 1
                END IF
                string_math oldy$, "*", LTRIM$(STR$(j&& - h%)), runningtotal$, terminating_decimal%, limit&&
                IF STRING$(LEN(z$), "0") = z$ AND runningtotal$ = "0" AND i&& >= LEN(decx$) THEN EXIT DO

                IF dp&& = 0 THEN ' Limited to && size unless converted to string.
                    IF i&& >= LEN(decx$) THEN
                        dp&& = INT(LEN(decx$) / 2 + .5)
                        IF dp&& = 0 THEN dp&& = -1
                    END IF
                END IF

                IF betatest% THEN PRINT "Sqrt "; sqrt$; " * 2 = ";: COLOR 2, 0: PRINT LTRIM$(STR$(VAL(sqrt$) * 2));: COLOR 7, 0: PRINT LTRIM$(STR$(j&& - h%)); " * "; LTRIM$(STR$(j&& - h%)); " ="; VAL(oldy$) * (j&& - h%)
                sqrt$ = sqrt$ + LTRIM$(STR$(j&& - h%))

                string_math oldy$, "*", LTRIM$(STR$(j&& - h%)), runningtotal$, terminating_decimal%, limit&&
                k$ = runningtotal$

                IF betatest% THEN PRINT "Remainder "; z$; " minus "; k$; " = ";
                EXIT FOR
            END IF
            oldy$ = y$
        NEXT

        IF betatest% THEN
            string_math stringmatha$, "-", stringmathb$, runningtotal$, terminating_decimal%, limit&&
            PRINT runningtotal$; " sqrt = "; sqrt$
        END IF

        i&& = i&& + 2
        IF LEN(z$) >= limit&& THEN EXIT DO
        x$ = x$ + "00"
    LOOP

    IF dp&& THEN
        sqrt$ = MID$(sqrt$, 0, dp&& + 1) + "." + MID$(sqrt$, dp&& + 1)
    END IF

END SUB

SUB string_math (stringmatha$, operator$, stringmathb$, runningtotal$, terminating_decimal%, limit&&)
    DIM AS _INTEGER64 a, c, aa, cc, s, ss

    SELECT CASE operator$
        CASE "+", "-"
            GOSUB string_add_subtract_new
        CASE "*"
            GOSUB string_multiply_new
        CASE "/"
            GOSUB string_divide
        CASE ELSE
            PRINT "Error, no operator selected. operator$ = "; operator$: END
    END SELECT
    EXIT SUB

    string_divide:
    terminating_decimal% = 0: divsign% = 0: divremainder& = 0: divremainder$ = "": divplace& = 0: divplace2& = 0: quotient$ = "": divcarry& = 0
    divbuffer& = LEN(stringmathb$) - LEN(stringmatha$)
    IF divbuffer& < 0 THEN divbuffer& = 0
    d2dividend$ = stringmatha$
    d1divisor$ = stringmathb$
    IF LEFT$(d1divisor$, 1) = "0" AND LEN(d1divisor$) = 1 THEN PRINT "Division by zero not allowed.": divsign% = 0: EXIT SUB
    IF LEFT$(d1divisor$, 1) = "-" THEN divsign% = -1: d1divisor$ = MID$(d1divisor$, 2)
    IF LEFT$(d2dividend$, 1) = "-" THEN
        IF divsign% THEN
            divsign% = 0
        ELSE
            divsign% = -1
        END IF
        d2dividend$ = MID$(d2dividend$, 2)
    END IF
    IF INSTR(d1divisor$, ".") <> 0 THEN
        DO UNTIL RIGHT$(d1divisor$, 1) <> "0"
            d1divisor$ = MID$(d1divisor$, 1, LEN(d1divisor$) - 1) ' Strip off trailing zeros
        LOOP
        divplace& = LEN(d1divisor$) - INSTR(d1divisor$, ".")
        d1divisor$ = MID$(d1divisor$, 1, INSTR(d1divisor$, ".") - 1) + MID$(d1divisor$, INSTR(d1divisor$, ".") + 1) ' Strip off decimal point.
        DO UNTIL LEFT$(d1divisor$, 1) <> "0"
            d1divisor$ = MID$(d1divisor$, 2) ' Strip off leading zeros for divisors smaller than .1
        LOOP
    END IF

    IF INSTR(d2dividend$, ".") <> 0 THEN
        d2dividend$ = d2dividend$ + STRING$(divplace& - LEN(d2dividend$) - INSTR(d2dividend$, "."), "0") ' Add any zeros based on the length of dividend at decimal - length of divisor at decimal. If less than zero, nothing added.
        divplace2& = INSTR(d2dividend$, ".")
        DO UNTIL RIGHT$(d2dividend$, 1) <> "0"
            d2dividend$ = MID$(d2dividend$, 1, LEN(d2dividend$) - 1) ' Strip off trailing zeros
        LOOP
        d2dividend$ = MID$(d2dividend$, 1, INSTR(d2dividend$, ".") - 1) + MID$(d2dividend$, INSTR(d2dividend$, ".") + 1) ' Strip off decimal point.
    ELSE
        d2dividend$ = d2dividend$ + STRING$(divplace&, "0") ' Add any zeros based on the length of dividend at decimal - length of divisor at decimal. If less than zero, nothing added.
        divplace& = 0
    END IF
    DO
        DO
            divremainder& = divremainder& + 1: divremainder$ = divremainder$ + MID$(d2dividend$, divremainder&, 1)
            IF MID$(d2dividend$, divremainder&, 1) = "" THEN
                IF divremainder$ = STRING$(LEN(divremainder$), "0") AND LEN(quotient$) > LEN(d2dividend$) THEN
                    divflag% = -1
                    terminating_decimal% = -1
                    EXIT DO
                END IF
                divcarry& = divcarry& + 1
                IF divcarry& = 1 THEN divplace3& = divremainder& - 1
                IF divcarry& > limit&& + 1 + divbuffer& THEN
                    divflag% = -2: EXIT DO
                END IF
                divremainder$ = divremainder$ + "0" ' No more digits to bring down.
            END IF
            IF LEN(divremainder$) > LEN(d1divisor$) OR LEN(divremainder$) = LEN(d1divisor$) AND divremainder$ >= d1divisor$ THEN EXIT DO
            quotient$ = quotient$ + "0"
        LOOP
        IF divflag% THEN divflag% = 0: EXIT DO

        FOR div_i% = 9 TO 1 STEP -1
            stringmatha$ = LTRIM$(STR$(div_i%)): stringmathb$ = d1divisor$
            GOSUB string_multiply_new ' Gets runningtotal$
            tempcutd$ = divremainder$ ' divremainder$ can be 00 or other leading zero values.
            DO
                IF LEN(tempcutd$) = 1 THEN EXIT DO
                IF LEFT$(tempcutd$, 1) = "0" THEN
                    tempcutd$ = MID$(tempcutd$, 2)
                ELSE
                    EXIT DO
                END IF
            LOOP
            IF LEN(tempcutd$) > LEN(runningtotal$) OR LEN(tempcutd$) = LEN(runningtotal$) AND runningtotal$ <= tempcutd$ THEN EXIT FOR
        NEXT
        quotient$ = quotient$ + LTRIM$(STR$(div_i%))
        stringmatha$ = LTRIM$(STR$(div_i%)): stringmathb$ = d1divisor$
        GOSUB string_multiply_new ' Gets runningtotal$
        stringmatha$ = divremainder$: stringmathb$ = runningtotal$
        operator$ = "-": GOSUB string_add_subtract_new
        divremainder$ = runningtotal$
    LOOP

    IF divplace& = 0 AND divplace2& = 0 THEN divplace& = divplace3&
    IF divplace2& THEN divplace& = divplace& + divplace2& - 1
    IF quotient$ = "" THEN divplace& = 0 ' dividend is zero.
    IF divplace& OR divplace2& THEN
        quotient$ = MID$(quotient$, 1, divplace&) + "." + MID$(quotient$, divplace& + 1)
        DO UNTIL RIGHT$(quotient$, 1) <> "0"
            quotient$ = MID$(quotient$, 1, LEN(quotient$) - 1) ' Strip off trailing zeros
        LOOP
        IF RIGHT$(quotient$, 1) = "." THEN quotient$ = MID$(quotient$, 1, LEN(quotient$) - 1) ' Strip off abandoned decimal.
    END IF
    DO UNTIL LEFT$(quotient$, 1) <> "0"
        quotient$ = MID$(quotient$, 2) ' Strip off leading zeros
    LOOP
    IF quotient$ = "" THEN quotient$ = "0": divsign% = 0
    stringmathb$ = quotient$: quotient$ = ""

    IF stringmathb$ = "overflow" THEN divsign% = 0: EXIT SUB

    runningtotal$ = stringmathb$: stringmathb$ = ""
    IF divsign% THEN runningtotal$ = "-" + runningtotal$

    IF stringmathround$ <> "" THEN runningtotal$ = runningtotal$ + stringmathround$
    RETURN

    string_add_subtract_new:
    a1$ = stringmatha$: b1$ = stringmathb$
    s = 18: i&& = 0: c = 0

    a$ = stringmatha$: b$ = stringmathb$: op$ = operator$

    IF op$ = "-" THEN
        IF LEFT$(b$, 1) = "-" THEN b$ = MID$(b$, 2) ELSE b$ = "-" + b$
    END IF

    IF INSTR(a$, ".") <> 0 OR INSTR(b$, ".") <> 0 THEN
        decimal% = -1
        IF INSTR(a$, ".") <> 0 THEN
            dec_a&& = LEN(MID$(a$, INSTR(a$, ".") + 1))
            a$ = MID$(a$, 1, INSTR(a$, ".") - 1) + MID$(a$, INSTR(a$, ".") + 1)
        END IF
        IF INSTR(b$, ".") <> 0 THEN
            dec_b&& = LEN(MID$(b$, INSTR(b$, ".") + 1))
            b$ = MID$(b$, 1, INSTR(b$, ".") - 1) + MID$(b$, INSTR(b$, ".") + 1)
        END IF
        ' Line up decimal places by inserting trailing zeros.
        IF dec_b&& > dec_a&& THEN
            j&& = dec_b&&
            a$ = a$ + STRING$(dec_b&& - dec_a&&, "0")
        ELSE
            j&& = dec_a&&
            b$ = b$ + STRING$(dec_a&& - dec_b&&, "0")
        END IF
    END IF

    IF LEFT$(a$, 1) = "-" OR LEFT$(b$, 1) = "-" THEN
        IF LEFT$(a$, 1) = "-" AND LEFT$(b$, 1) = "-" THEN
            sign$ = "--": a$ = MID$(a$, 2): b$ = MID$(b$, 2)
        ELSE
            IF LEFT$(a$, 1) = "-" THEN a$ = MID$(a$, 2): sign_a$ = "-"
            IF LEFT$(b$, 1) = "-" THEN b$ = MID$(b$, 2): sign_b$ = "-"

            IF LEFT$(a1$, 1) = "-" THEN a1_x$ = MID$(a1$, 2) ELSE a1_x$ = a1$
            IF LEFT$(b1$, 1) = "-" THEN b1_x$ = MID$(b1$, 2) ELSE b1_x$ = b1$

            string_compare a1_x$, b1_x$, gl%

            IF gl% < 0 THEN
                IF LEN(sign_b$) THEN sign$ = "-": SWAP a$, b$
            ELSE
                IF LEN(sign_a$) THEN sign$ = "-": SWAP sign_a$, sign_b$
            END IF
        END IF
    END IF

    z$ = ""

    DO
        i&& = i&& + s
        x1$ = MID$(a$, LEN(a$) - i&& + 1, s)
        x2$ = MID$(b$, LEN(b$) - i&& + 1, s)
        zeros% = LEN(x1$): IF LEN(x2$) > zeros% THEN zeros% = LEN(x2$)
        a = VAL(sign_a$ + x1$) + VAL(sign_b$ + x2$) + c
        IF x1$ + x2$ = "" AND c = 0 THEN EXIT DO
        c = 0
        IF a > VAL(STRING$(s, "9")) THEN a = a - 10 ^ s: c = 1
        IF a < 0 THEN a = a + 10 ^ s: c = -1
        tmp$ = LTRIM$(STR$(a))
        z$ = STRING$(zeros% - LEN(tmp$), "0") + tmp$ + z$
    LOOP

    IF decimal% THEN
        z$ = MID$(z$, 1, LEN(z$) - j&&) + "." + MID$(z$, LEN(z$) - j&& + 1)
    END IF

    ' Remove any leading zeros.
    DO
        IF LEFT$(z$, 1) = "0" THEN z$ = MID$(z$, 2) ELSE EXIT DO
    LOOP

    IF z$ = "" OR z$ = "0" THEN z$ = "0" ELSE z$ = LEFT$(sign$, 1) + z$

    runningtotal$ = z$

    sign$ = "": sign_a$ = "": sign_b$ = "": i&& = 0: j&& = 0: decimal% = 0: c = 0
    RETURN

    string_multiply_new:
    z$ = "": sign$ = "": mult&& = 0: h&& = 0: i&& = 0: j&& = 0: c = 0: decimal% = 0
    zz$ = "": ii&& = 0: jj&& = 0
    s = 8: ss = 18

    a$ = stringmatha$: b$ = stringmathb$

    IF INSTR(a$, "-") <> 0 OR INSTR(b$, "-") <> 0 THEN
        IF INSTR(a$, "-") <> 0 AND INSTR(b$, "-") <> 0 THEN
            a$ = MID$(a$, 2): b$ = MID$(b$, 2)
        ELSE
            IF INSTR(a$, "-") <> 0 THEN a$ = MID$(a$, 2) ELSE b$ = MID$(b$, 2)
            sign$ = "-"
        END IF
    END IF

    IF INSTR(a$, ".") <> 0 OR INSTR(b$, ".") <> 0 THEN
        decimal% = -1
        IF INSTR(a$, ".") <> 0 THEN
            dec_a&& = LEN(MID$(a$, INSTR(a$, ".") + 1))
            a$ = MID$(a$, 1, INSTR(a$, ".") - 1) + MID$(a$, INSTR(a$, ".") + 1)
        END IF
        IF INSTR(b$, ".") <> 0 THEN
            dec_b&& = LEN(MID$(b$, INSTR(b$, ".") + 1))
            b$ = MID$(b$, 1, INSTR(b$, ".") - 1) + MID$(b$, INSTR(b$, ".") + 1)
        END IF
    END IF

    IF LEN(a$) < LEN(b$) THEN SWAP a$, b$ ' Needed so x1$ is always the largest for leading zero replacements.

    DO
        h&& = h&& + s: i&& = 0
        x2$ = MID$(b$, LEN(b$) - h&& + 1, s)
        WHILE -1
            i&& = i&& + s
            x1$ = MID$(a$, LEN(a$) - i&& + 1, s)
            a = VAL(sign_a$ + x1$) * VAL(sign_b$ + x2$) + c
            c = 0
            tmp$ = LTRIM$(STR$(a))
            IF LEN(tmp$) > s THEN c = VAL(MID$(tmp$, 1, LEN(tmp$) - s)): tmp$ = MID$(tmp$, LEN(tmp$) - s + 1)
            z$ = STRING$(LEN(x1$) - LEN(tmp$), "0") + tmp$ + z$
            IF i&& >= LEN(a$) AND c = 0 THEN EXIT WHILE
        WEND

        jj&& = jj&& + 1

        IF jj&& > 1 THEN
            ii&& = 0: cc = 0
            aa$ = holdaa$
            bb$ = z$ + STRING$((jj&& - 1) * s, "0")

            DO
                ii&& = ii&& + ss
                xx1$ = MID$(aa$, LEN(aa$) - ii&& + 1, ss)
                xx2$ = MID$(bb$, LEN(bb$) - ii&& + 1, ss)
                aa = VAL(xx1$) + VAL(xx2$) + cc
                IF xx1$ + xx2$ = "" AND cc = 0 THEN EXIT DO ' Prevents leading zeros.
                cc = 0
                IF aa > VAL(STRING$(ss, "9")) THEN aa = aa - 10 ^ ss: cc = 1
                tmp$ = LTRIM$(STR$(aa))
                zz$ = STRING$(LEN(xx1$) - LEN(tmp$), "0") + tmp$ + zz$
            LOOP

            DO WHILE LEFT$(zz$, 1) = "0"
                IF LEFT$(zz$, 1) = "0" THEN zz$ = MID$(zz$, 2)
            LOOP
            IF zz$ = "" THEN zz$ = "0"

            holdaa$ = zz$
        ELSE
            holdaa$ = z$ + STRING$(jj&& - 1, "0")
        END IF

        z$ = "": zz$ = ""

    LOOP UNTIL h&& >= LEN(b$)

    z$ = holdaa$

    IF decimal% THEN
        DO UNTIL LEN(z$) >= dec_a&& + dec_b&&
            z$ = "0" + z$
        LOOP

        z$ = MID$(z$, 0, LEN(z$) - (dec_a&& + dec_b&& - 1)) + "." + MID$(z$, LEN(z$) - (dec_a&& + dec_b&&) + 1)

        DO UNTIL RIGHT$(z$, 1) <> "0" AND RIGHT$(z$, 1) <> "."
            z$ = MID$(z$, 1, LEN(z$) - 1)
        LOOP
    END IF

    IF z$ = "" OR z$ = "0" THEN z$ = "0" ELSE z$ = sign$ + z$

    decimal% = 0: sign$ = ""

    runningtotal$ = z$
    RETURN
END SUB

SUB string_compare (compa$, compb$, gl%)
    DO
        ' Remove trailing zeros after a decimal point.
        IF INSTR(compa$, ".") THEN
            DO UNTIL RIGHT$(compa$, 1) <> "0" AND RIGHT$(compa$, 1) <> "." AND RIGHT$(compa$, 1) <> "-"
                compa$ = MID$(compa$, 1, LEN(compa$) - 1)
            LOOP
        END IF
        IF INSTR(compb$, ".") THEN
            DO UNTIL RIGHT$(compb$, 1) <> "0" AND RIGHT$(compb$, 1) <> "." AND RIGHT$(compb$, 1) <> "-"
                compb$ = MID$(compb$, 1, LEN(compb$) - 1)
            LOOP
        END IF

        IF MID$(compa$, 1, 2) = "-0" OR compa$ = "" OR compa$ = "-" THEN compa$ = "0"
        IF MID$(compb$, 1, 2) = "-0" OR compb$ = "" OR compb$ = "-" THEN compb$ = "0"

        ' A - and +
        IF LEFT$(compa$, 1) = "-" THEN j% = -1
        IF LEFT$(compb$, 1) = "-" THEN k% = -1
        IF k% = 0 AND j% THEN gl% = -1: EXIT DO
        IF j% = 0 AND k% THEN gl% = 1: EXIT DO

        ' A decimal and non-decimal.
        j% = INSTR(compa$, ".")
        k% = INSTR(compb$, ".")

        IF j% = 0 AND k% THEN
            IF compa$ = "0" THEN gl% = -1 ELSE gl% = 1
            EXIT DO
        END IF
        IF k% = 0 AND j% THEN
            IF compb$ = "0" THEN gl% = 1 ELSE gl% = -1
            EXIT DO
        END IF

        ' Both decimals.
        IF j% THEN
            IF compa$ > compb$ THEN
                gl% = 1
            ELSEIF compa$ = compb$ THEN gl% = 0
            ELSEIF compa$ < compb$ THEN gl% = -1
            END IF
            EXIT DO
        END IF

        ' Both positive or both negative whole numbers.
        SELECT CASE LEN(compa$)
            CASE IS < LEN(compb$)
                gl% = -1
            CASE IS = LEN(compb$)
                IF compa$ = compb$ THEN
                    gl% = 0
                ELSEIF compa$ > compb$ THEN gl% = 1
                ELSEIF compa$ < compb$ THEN gl% = -1
                END IF
            CASE IS > LEN(compb$)
                gl% = 1
        END SELECT
        EXIT DO
    LOOP
END SUB

Pete
Reply
#17
Pete
I converted the Ramanujan Pi routine to use Treebeard's string-math, maybe you can adapt it to use your string-math routines
Code: (Select All)
Dim As String n, m, c, sum, f, f4, f4k, c1, c2, c3, c34k, t1, t2, t3
Dim As Long k, k4
Dim t As Double
t = Timer
digits% = 100
bInit
c1 = "1103"
c2 = "26390"
c3 = "396"
f = "1"
f4k = "1"
sum = "1103"
c34k = "1"
k4 = 0
t1 = c3
t2 = c3
bMul t1, t2, c3 'result in c3
t1 = c3
t2 = c3
bMul t1, t2, c3 'result in c3
For k = 1 To digits% / 7.984
    t1 = f: bMul Str$(k), t1, f 'result in f
    t1 = f: t2 = f
    bMul t1, t2, f4 'result in f4
    t1 = f4: t2 = f4
    bMul t1, t2, f4 'result in f4
    t1 = c34k
    bMul c3, t1, c34k 'result in c34k
    t1 = Str$(k4 + 1)
    t2 = f4k
    bMul t1, t2, f4k 'result in f4k
    t1 = Str$(k4 + 2)
    t2 = f4k
    bMul t1, t2, f4k 'result in f4k
    t1 = Str$(k4 + 3)
    t2 = f4k
    bMul t1, t2, f4k 'result in f4k
    t1 = Str$(k4 + 4)
    t2 = f4k
    bMul t1, t2, f4k 'result in f4k
    k4 = k4 + 4
    t1 = Str$(k)
    bMul t1, c2, t2 'result in t2
    bAdd c1, t2, t1 'result in t1
    bMul f4k, t1, t2 'result in t2
    bMul f4, c34k, t1 'result in t1
    bDiv t2, t1, t3 'result in t3
    t1 = sum
    bAdd t1, t3, sum 'result in sum
Next
bSqr "8", t1 'result in t1
t2 = "9801"
bDiv t1, t2, t3 'result in t3
bMul t3, sum, t2 'result in t2
bDiv t1, t2, t3 'result in t3
bDiv "1", t2, t1 'result in t1
Print t1
t = Timer - t
Reply
#18
Your endless string formulas are too complicated and cumbersome. There I get headdache.  Angry

If one implement the formulas it is short and concise.  Wink  I have now implemented a guide for C++ (running) in C. It's Chudnovsky's formula for pi, but it won't work in QB64. Does anyone see where the error is? - I suspect it has to do with "pi += ..." - maybe.

Chudnovskys Pi formula in C
Code: (Select All)
//PI nach Cudnowsky - https://stackoverflow.com/questions/12028313/c-chudnovsky-formula-for-pi
//24. Aug. 2022, in C 25. Aug. 2022

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

long double fac(double num);

int main(void)
{
  long double pi = 0.0;
 
    for (double k = 0.0; k < 10.0; k++)
    {
    pi += (pow(-1.0, k) * fac(6.0 * k) * (13591409.0 + (545140134.0 * k)))
       / (fac(3.0 * k) * pow(fac(k), 3.0) * pow(640320.0, 3.0 * k + 3.0/2.0));        
  }
  pi *= 12.0;
    __mingw_printf("%.12Lf\n\n", 1.0 / pi);
 
  return 0;
}

long double fac(double num)
{
    double result = 1.0;
  for (double i = 2.0; i < num; i++)
        { result *= i; }
  return result;
}

That has to be! Unfortunately there is no explanation why.
Code: (Select All)
__mingw_printf("%.12Lf\n\n", 1.0 / pi);

Chudnovskys Pi formula in QB64
Code: (Select All)
'Pi nach Chudnovsky berechnen - 25. Aug. 2022

Option _Explicit

Declare Function fac (num as Double) as Double

Dim As Double pi, k

pi = 0.0: k = 0.0
While k < 10.0
  k = k + 1.0
  pi = pi + (-1.0 ^ k) * fac(6.0 * k) * (13591409.0 + (545140134.0 * k)) / fac(3.0 * k) * fac(k) ^ 3.0 * 640320.0 ^ 3.0 * k + (3.0 / 2.0)
  'pi = pi + pi
Wend
pi = pi * 12.0

Print
Print 1.0 / pi

End

Function fac (num As Double)

  Dim As Double resultat, i
  resultat = 1.0: i = 2.0

  While i < num
    i = i + i
    resultat = resultat * i
  Wend
  fac = resultat

End Function
Reply
#19
^
|
|
Why not just "printf" instead of "__mingw_printf"? Because you're including "stdio.h" then "printf()" should be available. I tried to compile your program as it was, "gcc" asked me if I meant "builtin_printf()" or alike. After I forgot what option for linker so it links to math library...
Reply
#20
(08-25-2022, 11:43 AM)Jack Wrote: Pete
I converted the Ramanujan Pi routine to use Treebeard's string-math, maybe you can adapt it to use your string-math routines
Code: (Select All)
Dim As String n, m, c, sum, f, f4, f4k, c1, c2, c3, c34k, t1, t2, t3
Dim As Long k, k4
Dim t As Double
t = Timer
digits% = 100
bInit
c1 = "1103"
c2 = "26390"
c3 = "396"
f = "1"
f4k = "1"
sum = "1103"
c34k = "1"
k4 = 0
t1 = c3
t2 = c3
bMul t1, t2, c3 'result in c3
t1 = c3
t2 = c3
bMul t1, t2, c3 'result in c3
For k = 1 To digits% / 7.984
    t1 = f: bMul Str$(k), t1, f 'result in f
    t1 = f: t2 = f
    bMul t1, t2, f4 'result in f4
    t1 = f4: t2 = f4
    bMul t1, t2, f4 'result in f4
    t1 = c34k
    bMul c3, t1, c34k 'result in c34k
    t1 = Str$(k4 + 1)
    t2 = f4k
    bMul t1, t2, f4k 'result in f4k
    t1 = Str$(k4 + 2)
    t2 = f4k
    bMul t1, t2, f4k 'result in f4k
    t1 = Str$(k4 + 3)
    t2 = f4k
    bMul t1, t2, f4k 'result in f4k
    t1 = Str$(k4 + 4)
    t2 = f4k
    bMul t1, t2, f4k 'result in f4k
    k4 = k4 + 4
    t1 = Str$(k)
    bMul t1, c2, t2 'result in t2
    bAdd c1, t2, t1 'result in t1
    bMul f4k, t1, t2 'result in t2
    bMul f4, c34k, t1 'result in t1
    bDiv t2, t1, t3 'result in t3
    t1 = sum
    bAdd t1, t3, sum 'result in sum
Next
bSqr "8", t1 'result in t1
t2 = "9801"
bDiv t1, t2, t3 'result in t3
bMul t3, sum, t2 'result in t2
bDiv t1, t2, t3 'result in t3
bDiv "1", t2, t1 'result in t1
Print t1
t = Timer - t

Nice. Tested it to 50-digits. All but the last digit matches, but that's just rounding on the online comparison. The final digit is always suspect, so the results match through digit 49, perfectly.

Code: (Select All)
WIDTH 160, 42
_SCREENMOVE 0, 0
' Treebeard's String Math +-*/
CONST neg$ = "-"
CONST negative = -1
CONST positive = 1
CONST asc0 = 48
CONST dp$ = "."
CONST zero$ = "0"
CONST one$ = "1"
CONST two$ = "2"
CONST three$ = "3"
CONST four$ = "4"
CONST five$ = "5"
CONST False = 0
CONST True = -1
CONST basechr = "@"
CONST basesep$ = ","
CONST maxlongdig = 8
CONST emem = 32
CONST memget = 0
CONST memput = 1
CONST defaultdigits = 30
CONST maxmem = 35
CONST maxstack = 10
CONST minconst = 30
CONST maxconst = 35
CONST pimem = 30
CONST pi2mem = 31
CONST phimem = 33
CONST ln10mem = 34
CONST ln2mem = 35
CONST memclr = 2

'useful shared stuff, initialize these in bInit()
DIM SHARED errormsg$, abortmsg$, Error$, bmem$(maxmem), out$
DIM SHARED zmem$(maxstack), cname$(maxconst)
DIM SHARED bncpath$, prmcntfile$
DIM SHARED digits%, zstack%

'Prime count table data
DIM maxprmcnt%
DIM prmcnt&
'--------------------------------------------
DIM AS STRING n, m, c, sum, f, f4, f4k, c1, c2, c3, c34k, t1, t2, t3
DIM AS LONG k, k4
DIM t AS DOUBLE

t = TIMER
digits% = 50
bInit
c1 = "1103"
c2 = "26390"
c3 = "396"
f = "1"
f4k = "1"
sum = "1103"
c34k = "1"
k4 = 0
t1 = c3
t2 = c3
bMul t1, t2, c3 'result in c3
t1 = c3
t2 = c3
bMul t1, t2, c3 'result in c3

FOR k = 1 TO digits% / 7.984
    t1 = f: bMul STR$(k), t1, f 'result in f
    t1 = f: t2 = f
    bMul t1, t2, f4 'result in f4
    t1 = f4: t2 = f4
    bMul t1, t2, f4 'result in f4
    t1 = c34k
    bMul c3, t1, c34k 'result in c34k
    t1 = STR$(k4 + 1)
    t2 = f4k
    bMul t1, t2, f4k 'result in f4k
    t1 = STR$(k4 + 2)
    t2 = f4k
    bMul t1, t2, f4k 'result in f4k
    t1 = STR$(k4 + 3)
    t2 = f4k
    bMul t1, t2, f4k 'result in f4k
    t1 = STR$(k4 + 4)
    t2 = f4k
    bMul t1, t2, f4k 'result in f4k
    k4 = k4 + 4
    t1 = STR$(k)
    bMul t1, c2, t2 'result in t2
    bAdd c1, t2, t1 'result in t1
    bMul f4k, t1, t2 'result in t2
    bMul f4, c34k, t1 'result in t1
    bDiv t2, t1, t3 'result in t3
    t1 = sum
    bAdd t1, t3, sum 'result in sum
    oldt1$ = t1: oldt3$ = t3: oldsum$ = sum
    CALL pi(t1, sum, k)
    t1 = oldt1$: t3 = oldt3$: sum = oldsum$
NEXT
COLOR 14, 0: PRINT "Ramanujan pi = 3.1415926535897932384626433832795028841971693993751"
COLOR 7, 0
t = TIMER - t
END

SUB pi (t1$, sum$, k)
    bSqr "8", t1$ 'result in t1
    t2$ = "9801"
    bDiv t1$, t2$, t3$ 'result in t3
    bMul t3$, sum$, t2$ 'result in t2
    bDiv t1$, t2$, t3$ 'result in t3
    bDiv "1", t2$, t1$ 'result in t1
    PRINT "loop #"; LTRIM$(STR$(k));: LOCATE , 10: PRINT " pi = "; t1$
END SUB

SUB bAdd (s1$, s2$, out$)
    DIM last1%, dp1%, sign1%, last2%, dp2%, sign2%
    DIM last%, d1%, d2%, dpt%, carry%
    DIM i%, n%

    'strip the numbers
    bStripDp s1$, last1%, dp1%, sign1%
    bStripDp s2$, last2%, dp2%, sign2%

    'treat different signs as subtraction and exit
    IF sign1% = negative AND sign2% = positive THEN
        bSub s2$, s1$, out$
        bNeg s1$
        EXIT SUB
    ELSEIF sign1% = positive AND sign2% = negative THEN
        bSub s1$, s2$, out$
        bNeg s2$
        EXIT SUB
    END IF

    'align the decimal points and digit pointers
    last% = bMaxInt%(last1% - dp1%, last2% - dp2%)
    d1% = last% + dp1%
    d2% = last% + dp2%
    dpt% = bMaxInt%(dp1%, dp2%)
    last% = dpt% + last%
    out$ = SPACE$(last%)
    carry% = 0

    'do the addition right to left
    FOR i% = last% TO 1 STEP -1
        IF i% <> dpt% THEN
            n% = carry%
            IF d1% > 0 THEN n% = n% + VAL(MID$(s1$, d1%, 1))
            IF d2% > 0 THEN n% = n% + VAL(MID$(s2$, d2%, 1))
            carry% = n% \ 10
            MID$(out$, i%, 1) = CHR$(asc0 + (n% MOD 10))
        ELSE
            MID$(out$, i%, 1) = dp$
        END IF
        d1% = d1% - 1
        d2% = d2% - 1
    NEXT i%
    IF carry% THEN out$ = one$ + out$

    'clean up
    IF sign1% = negative THEN s1$ = neg$ + s1$: s2$ = neg$ + s2$: out$ = neg$ + out$
    bClean s1$
    bClean s2$
    bClean out$
END SUB

SUB bSub (s1$, s2$, out$)
    DIM last1%, dp1%, sign1%
    DIM last2%, dp2%, sign2%
    DIM last%, d1%, d2%, dpt%, borrow%, swapflag%
    DIM i%, n%

    'strip the numbers
    bStripDp s1$, last1%, dp1%, sign1%
    bStripDp s2$, last2%, dp2%, sign2%

    'treat different signs as addition
    IF sign1% = negative AND sign2% = positive THEN
        bNeg s1$
        bNeg s2$
        bAdd s1$, s2$, out$
        bNeg s2$
        EXIT SUB
    ELSEIF sign1% = positive AND sign2% = negative THEN
        bAdd s1$, s2$, out$
        bNeg s2$
        EXIT SUB
    END IF

    'align the decimal points and digit pointers
    last% = bMaxInt%(last1% - dp1%, last2% - dp2%)
    d1% = last% + dp1%
    d2% = last% + dp2%
    dpt% = bMaxInt%(dp1%, dp2%)
    last% = dpt% + last%
    out$ = SPACE$(last%)
    borrow% = 0

    'always subtract smaller from bigger to avoid complements
    IF bIsMore%(s2$, s1$) THEN
        bSwapString s1$, s2$
        bSwapInt d2%, d1%
        swapflag% = True
    END IF

    'do the subtraction right to left
    FOR i% = last% TO 1 STEP -1
        IF i% <> dpt% THEN
            IF d1% > 0 THEN n% = VAL(MID$(s1$, d1%, 1)) ELSE n% = 0
            IF d2% > 0 THEN n% = n% - VAL(MID$(s2$, d2%, 1))
            n% = n% - borrow%
            IF n% >= 0 THEN borrow% = 0 ELSE borrow% = 1: n% = n% + 10
            MID$(out$, i%, 1) = CHR$(asc0 + n%)
        ELSE
            MID$(out$, i%, 1) = dp$
        END IF
        d1% = d1% - 1
        d2% = d2% - 1
    NEXT i%

    'clean up
    IF sign1% = negative THEN s1$ = neg$ + s1$: s2$ = neg$ + s2$
    IF swapflag% THEN
        bSwapString s1$, s2$
        sign1% = -sign1%
    END IF
    IF sign1% = negative THEN out$ = neg$ + out$
    bClean s1$
    bClean s2$
    bClean out$

END SUB

SUB bSqr (s$, out$)
    DIM dvd$, div$, dig$, newdiv$, t$, z$
    DIM slog%, ssign%, slen%, spt%, olddigits%, n%, m%

    IF bIsNeg%(s$) THEN out$ = errormsg$: EXIT SUB

    'strip to whole number + group digits by 2 left or right of decimal
    bLogGet s$, slog%, ssign%, True
    slen% = LEN(s$)
    IF slog% MOD 2 THEN spt% = 2 ELSE spt% = 1

    'Force at least enough digits to show integer of root
    olddigits% = digits%
    n% = 1 + slog% \ 2
    IF digits% < n% THEN digits% = n%

    'figure first digit and setup loop
    n% = VAL(LEFT$(s$ + "0", spt%))
    m% = INT(SQR(n%))
    out$ = LTRIM$(STR$(m%))
    dvd$ = LTRIM$(STR$(n% - m% * m%))
    spt% = spt% + 1

    DO
        'all done?
        IF (spt% > slen% AND bIsZero%(dvd$)) OR LEN(out$) >= digits% THEN EXIT DO

        'append next 2 digits (or 0s) to dividend
        dvd$ = dvd$ + LEFT$(MID$(s$, spt%, 2) + "00", 2)
        spt% = spt% + 2

        'divisor=twice the root * 10
        z$ = out$
        bAdd out$, z$, div$
        bShift div$, 1

        'estimate divisor, and adjust if too big.  Unit is next digit of root.
        bDivInt dvd$, div$, dig$
        DO
            bAdd div$, dig$, newdiv$
            bMul newdiv$, dig$, t$
            IF NOT bIsMore%(t$, dvd$) THEN EXIT DO
            bInc dig$, -1
        LOOP
        out$ = out$ + dig$

        'form new divisor
        z$ = dvd$
        bSub z$, t$, dvd$

    LOOP

    'clean up
    bLogPut s$, slog%, ssign%
    IF slog% < 0 THEN slog% = slog% - 1
    bLogPut out$, slog% \ 2, ssign%
    digits% = olddigits%

END SUB

SUB bMul (s1$, s2$, out$)
    DIM t$
    DIM slog1%, sign1%, slog2%, sign2%, outdp%, outsign%, outlog%, swapflag%

    'strip multiplier
    t$ = s2$
    bLogGet t$, slog2%, sign2%, True

    'times 0
    IF t$ = zero$ THEN
        out$ = zero$

        'do powers of 10 with shifts
    ELSEIF t$ = one$ THEN
        out$ = s1$
        sign1% = bSign%(out$)
        IF sign1% = negative THEN bAbs out$
        bShift out$, slog2%
        IF sign1% <> sign2% THEN bNeg out$

        'the hard way
    ELSE
        'strip all
        s2$ = t$: t$ = ""
        bLogGet s1$, slog1%, sign1%, True

        'figure decimal point and sign of answer
        outdp% = bLogDp%(s1$, slog1%) + bLogDp%(s2$, slog2%)
        IF sign1% <> sign2% THEN outsign% = negative ELSE outsign% = positive

        'always multiply by the shorter number
        IF LEN(s2$) > LEN(s1$) THEN bSwapString s1$, s2$: swapflag% = True

        'do it
        IF LEN(s2$) <= maxlongdig THEN bMulLong s1$, s2$, out$ ELSE bMulChar s1$, s2$, out$

        'clean up
        outlog% = bLogDp%(out$, outdp%)
        bLogPut out$, outlog%, outsign%
        IF swapflag% THEN bSwapString s1$, s2$
        bLogPut s1$, slog1%, sign1%
        bLogPut s2$, slog2%, sign2%

    END IF

END SUB

SUB bMulChar (s1$, s2$, out$)
    DIM last1%, last2%, last%
    DIM i%, j%, k%, sj%, ej%
    DIM product&

    last1% = LEN(s1$)
    last2% = LEN(s2$)
    last% = last1% + last2%
    out$ = SPACE$(last%)
    product& = 0
    FOR i% = 0 TO last% - 1
        k% = last1% - i%
        sj% = 1 - k%: IF sj% < 0 THEN sj% = 0
        ej% = last1% - k%: IF ej% > last2% - 1 THEN ej% = last2% - 1
        FOR j% = sj% TO ej%
            product& = product& + VAL(MID$(s1$, k% + j%, 1)) * VAL(MID$(s2$, last2% - j%, 1))
        NEXT j%
        MID$(out$, last% - i%, 1) = CHR$(asc0 + CINT(product& MOD 10&))
        product& = product& \ 10&
    NEXT i%
    IF product& THEN out$ = LTRIM$(STR$(product&)) + out$
END SUB

SUB bMulLong (s1$, s2$, out$)
    DIM last1%, i%
    DIM s2val&, product&

    last1% = LEN(s1$)
    s2val& = VAL(s2$)
    out$ = SPACE$(last1%)
    FOR i% = last1% TO 1 STEP -1
        product& = product& + VAL(MID$(s1$, i%, 1)) * s2val&
        MID$(out$, i%, 1) = CHR$(asc0 + CINT(product& MOD 10&))
        product& = product& \ 10&
    NEXT i%
    IF product& THEN out$ = LTRIM$(STR$(product&)) + out$
END SUB

SUB bDivLong (s1$, s2$, quotient$, remainder$)
    DIM rmdr&, dividend&, divisor&
    DIM dig%, i%

    quotient$ = ""
    rmdr& = 0
    divisor& = VAL(s2$)

    FOR i% = 1 TO digits%
        dividend& = rmdr& * 10& + VAL(MID$(s1$, i%, 1))
        dig% = dividend& \ divisor&
        quotient$ = quotient$ + CHR$(asc0 + dig%)
        rmdr& = dividend& - dig% * divisor&
    NEXT i%

    IF LEN(quotient$) = 0 THEN quotient$ = zero$
    remainder$ = LTRIM$(STR$(rmdr&))

END SUB

SUB bDiv (s1$, s2$, out$)
    DIM t$
    DIM slog1%, sign1%, slog2%, sign2%
    DIM outlog%, outsign%, olddigits%

    'strip divisor
    t$ = s2$
    bLogGet t$, slog2%, sign2%, True

    'divide by zero?
    IF t$ = zero$ THEN
        out$ = Error$

        'do powers of 10 with shifts
    ELSEIF t$ = one$ THEN
        out$ = s1$
        sign1% = bSign%(out$)
        IF sign1% = negative THEN bAbs out$
        bShift out$, -slog2%
        IF sign1% <> sign2% THEN bNeg out$

        'the hard way
    ELSE
        'strip all
        s2$ = t$: t$ = ""
        bLogGet s1$, slog1%, sign1%, True

        'figure decimal point and sign of answer
        outlog% = slog1% + bLogDp%(s2$, slog2%)
        IF sign1% <> sign2% THEN outsign% = negative ELSE outsign% = positive

        'bump digits past leading zeros and always show whole quotient
        olddigits% = digits%
        digits% = digits% + LEN(s2$)
        IF digits% < outlog% + 1 THEN digits% = outlog% + 1

        'do it, ignore remainder
        IF LEN(s2$) <= maxlongdig THEN bDivLong s1$, s2$, out$, t$ ELSE bDivChar s1$, s2$, out$, t$

        'clean up
        bLogPut out$, outlog%, outsign%
        bLogPut s1$, slog1%, sign1%
        bLogPut s2$, slog2%, sign2%
        digits% = olddigits%
    END IF

END SUB

SUB bDivChar (s1$, s2$, quotient$, remainder$)
    DIM last1%, last2%, ldvd%, lrem%, dig%, borrow%
    DIM i%, j%, n%
    DIM dvd$

    last1% = LEN(s1$) 'length of the dividend
    last2% = LEN(s2$) 'length of the divisor
    quotient$ = ""
    remainder$ = ""

    FOR i% = 1 TO digits%
        'get next digit of dividend or zero$ if past end
        IF i% <= last1% THEN
            dvd$ = remainder$ + MID$(s1$, i%, 1)
        ELSE
            dvd$ = remainder$ + zero$
        END IF

        'if dividend < divisor then digit%=0 else have to calculate it.
        'do fast compare using string operations. see bComp%()
        bStripZero dvd$
        ldvd% = LEN(dvd$)
        IF (ldvd% < last2%) OR ((ldvd% = last2%) AND (dvd$ < s2$)) THEN
            'divisor is bigger, so digit is 0, easy!
            dig% = 0
            remainder$ = dvd$

        ELSE
            'dividend is bigger, but no more than 9 times bigger.
            'subtract divisor until we get remainder less than divisor.
            'time hog, average is 5 tries through j% loop.  There's a better way.
            FOR dig% = 1 TO 9
                remainder$ = ""
                borrow% = 0
                FOR j% = 0 TO ldvd% - 1
                    n% = last2% - j%
                    IF n% < 1 THEN n% = 0 ELSE n% = VAL(MID$(s2$, n%, 1))
                    n% = VAL(MID$(dvd$, ldvd% - j%, 1)) - n% - borrow%
                    IF n% >= 0 THEN borrow% = 0 ELSE borrow% = 1: n% = n% + 10
                    remainder$ = CHR$(asc0 + n%) + remainder$
                NEXT j%

                'if remainder < divisor then exit
                bStripZero remainder$
                lrem% = LEN(remainder$)
                IF (lrem% < last2%) OR ((lrem% = last2%) AND (remainder$ < s2$)) THEN EXIT FOR

                dvd$ = remainder$
                ldvd% = LEN(dvd$)
            NEXT dig%

        END IF
        quotient$ = quotient$ + CHR$(asc0 + dig%)
    NEXT i%

END SUB

SUB bLogGet (s$, slog%, sign%, zeroflag%)
    DIM dpt%, n%

    IF LEFT$(s$, 1) = neg$ THEN s$ = MID$(s$, 2): sign% = negative ELSE sign% = positive
    bStripZero s$
    dpt% = INSTR(s$, dp$)
    SELECT CASE dpt%
        CASE 0
            slog% = LEN(s$) - 1
        CASE 1
            n% = dpt% + 1
            DO WHILE MID$(s$, n%, 1) = zero$
                n% = n% + 1
            LOOP
            s$ = MID$(s$, n%)
            slog% = dpt% - n%
        CASE ELSE
            s$ = LEFT$(s$, dpt% - 1) + MID$(s$, dpt% + 1)
            slog% = dpt% - 2
    END SELECT

    'remove trailing 0's if zeroflag%
    IF zeroflag% THEN bStripTail s$

END SUB


SUB bLogPut (s$, slog%, sign%)
    DIM last%

    last% = LEN(s$)
    IF LEN(s$) = 0 OR s$ = zero$ THEN
        s$ = zero$
    ELSEIF slog% < 0 THEN
        s$ = dp$ + STRING$(-slog% - 1, zero$) + s$
    ELSEIF slog% > last% - 1 THEN
        s$ = s$ + STRING$(slog% - last% + 1, zero$) + dp$
    ELSE
        s$ = LEFT$(s$, slog% + 1) + dp$ + MID$(s$, slog% + 2)
    END IF
    bClean s$
    IF sign% = negative THEN s$ = neg$ + s$
END SUB

SUB bInt (s$)
    DIM n%

    n% = INSTR(s$, dp$)
    IF n% THEN
        IF n% = 1 THEN s$ = zero$ ELSE s$ = LEFT$(s$, n% - 1)
        IF s$ = neg$ OR LEFT$(s$, 2) = "-." THEN s$ = zero$
    END IF

END SUB

SUB bInit ()
    DIM i%

    'a few defaults
    'digits% = defaultdigits
    errormsg$ = "error"
    abortmsg$ = "abort"

    'clear memory
    zstack% = 0
    FOR i% = 0 TO maxmem
        bmem$(i%) = zero$
    NEXT i%
    FOR i% = 1 TO maxstack
        zmem$(i%) = zero$
    NEXT i%

    'useful constants
    cname$(pimem) = "pi": bmem$(pimem) = "3.14159265358979323846264338327"
    cname$(pi2mem) = "2pi": bmem$(pi2mem) = "6.28318530717958647692528676654"
    cname$(emem) = "e": bmem$(emem) = "2.71828182845904523536028747135"
    cname$(phimem) = "phi": bmem$(phimem) = "1.61803398874989484820458683436"
    cname$(ln10mem) = "ln(10)": bmem$(ln10mem) = "2.30258509299404568401799145468"
    cname$(ln2mem) = "ln(2)": bmem$(ln2mem) = ".693147180559945309417232121458"

    bncpath$ = "" 'path for files (or current dir if null)
    prmcntfile$ = "BNPRMCNT.DAT" 'prime count table
    '    LoadPrimeTable

END SUB

SUB bStr (s$, out$)
    DIM t$
    DIM n%, i%

    n% = INSTR(s$, ".")
    IF n% THEN t$ = MID$(s$, n% + 1) ELSE t$ = RIGHT$(s$, 1)
    out$ = ""
    FOR i% = 1 TO VAL(s$)
        out$ = t$ + out$
    NEXT i%
    IF LEN(out$) = 0 THEN out$ = zero$

END SUB

'Trim leading spaces, add decimal points, eliminate signs.
'Returns last%=length of string, dpt%=decimal place, sign%=-1 or 1.
'Called only by bAdd() and bSub() which needs a final decimal point.
'
SUB bStripDp (s$, last%, dpt%, sign%)
    IF LEFT$(s$, 1) = neg$ THEN s$ = MID$(s$, 2): sign% = negative ELSE sign% = positive
    bStripZero s$
    IF INSTR(s$, dp$) = 0 THEN s$ = s$ + dp$
    IF s$ = dp$ THEN s$ = "0."
    dpt% = INSTR(s$, dp$)
    last% = LEN(s$)
END SUB

'Strip trailing 0s to "." (but leave something)
'
SUB bStripTail (s$)
    DIM n%

    n% = LEN(s$)
    DO WHILE MID$(s$, n%, 1) = zero$
        n% = n% - 1
        IF n% <= 1 THEN EXIT DO
    LOOP
    IF n% THEN IF MID$(s$, n%, 1) = dp$ THEN n% = n% - 1
    s$ = LEFT$(s$, n%)
    IF LEN(s$) = 0 THEN s$ = zero$
END SUB

'Strip leading 0s and final "." (but leave something)
'
SUB bStripZero (s$)
    DIM n%

    n% = 1
    DO WHILE MID$(s$, n%, 1) = zero$
        n% = n% + 1
    LOOP
    IF n% > 1 THEN s$ = MID$(s$, n%)
    IF RIGHT$(s$, 1) = dp$ THEN s$ = LEFT$(s$, LEN(s$) - 1)
    IF LEN(s$) = 0 THEN s$ = zero$
END SUB

SUB bNeg (s$)
    IF LEFT$(s$, 1) = neg$ THEN s$ = MID$(s$, 2) ELSE s$ = neg$ + s$
END SUB

SUB bClean (s$)
    DIM sign%

    IF LEFT$(s$, 1) = neg$ THEN s$ = MID$(s$, 2): sign% = True
    bStripZero s$
    IF INSTR(s$, dp$) THEN bStripTail s$
    IF sign% AND s$ <> zero$ THEN s$ = neg$ + s$

END SUB


SUB bSwapInt (s1%, s2%)
    DIM t%

    t% = s1%
    s1% = s2%
    s2% = t%
END SUB

SUB bSwapString (s1$, s2$)
    DIM t$

    t$ = s1$
    s1$ = s2$
    s2$ = t$
END SUB

SUB bShift (s$, n%)
    DIM slog%, sign%

    bLogGet s$, slog%, sign%, False
    bLogPut s$, slog% + n%, sign%
END SUB

SUB bDivInt (s1$, s2$, out$)
    DIM t$

    bDivIntMod s1$, s2$, out$, t$
END SUB

SUB bDivIntMod (s1$, s2$, quotient$, remainder$)
    DIM slog1%, sign1%, slog2%, sign2%
    DIM olddigits%, outlog%, outsign%

    olddigits% = digits%

    'strip the numbers, set flag false to NOT trim zeros, slower but needed
    bLogGet s2$, slog2%, sign2%, False
    IF s2$ = zero$ THEN quotient$ = Error$: remainder$ = Error$: EXIT SUB
    bLogGet s1$, slog1%, sign1%, False

    'figure decimal point and sign of answer
    outlog% = slog1% + bLogDp%(s2$, slog2%)
    IF sign1% <> sign2% THEN outsign% = negative ELSE outsign% = positive

    'a trick: figure the decimal and only find that many digits
    digits% = outlog% + 1

    'send the work out
    IF LEN(s2$) <= maxlongdig THEN bDivLong s1$, s2$, quotient$, remainder$ ELSE bDivChar s1$, s2$, quotient$, remainder$

    'clean up
    bLogPut s1$, slog1%, sign1%
    bLogPut s2$, slog2%, sign2%
    bClean quotient$
    bClean remainder$
    IF sign1% <> sign2% THEN bNeg quotient$
    digits% = olddigits%

END SUB

SUB bInc (s$, num%)
    DIM dig%, n%, borrow%

    IF num% = 0 THEN EXIT SUB
    dig% = INSTR(s$, dp$)
    IF dig% THEN dig% = dig% - 1 ELSE dig% = LEN(s$)
    n% = num%
    IF n% > 0 THEN 'increment (n>0)
        DO WHILE n%
            IF dig% < 1 THEN
                s$ = LTRIM$(STR$(n%)) + s$
                n% = 0
            ELSE
                n% = n% + VAL(MID$(s$, dig%, 1))
                MID$(s$, dig%, 1) = CHR$(asc0 + (n% MOD 10))
                n% = n% \ 10
                dig% = dig% - 1
            END IF
        LOOP
    ELSE 'decrement (n<0)
        n% = -n%
        DO WHILE n%
            IF dig% < 1 THEN s$ = zero$: EXIT DO
            borrow% = 0
            n% = VAL(MID$(s$, dig%, 1)) - n%
            DO WHILE n% < 0
                n% = n% + 10: borrow% = borrow% + 1
            LOOP
            MID$(s$, dig%, 1) = CHR$(asc0 + n%)
            n% = borrow%
            dig% = dig% - 1
        LOOP
    END IF
    bStripZero s$
END SUB

SUB bAbs (s$)
    IF LEFT$(s$, 1) = neg$ THEN s$ = MID$(s$, 2)
END SUB

SUB bMod (s1$, s2$, out$)
    DIM t$

    bDivIntMod s1$, s2$, t$, out$
END SUB

SUB bModPower (s1$, s2$, s3$, out$)
    'Use variation of "Russian Peasant Method" to figure m=(c^d) mod n.
    'Byte, Jan 83, p.206.
    'test value: (71611947 ^ 63196467) mod 94815109 = 776582

    'm=1
    'do
    '  if d is odd then m=(m*c) mod n
    '  c=(c*c) mod n
    '  d=int(d/2)
    'loop while d>0
    'm is the answer

    DIM c$, d$, z$, w$
    STATIC n$ 'remember modulus for next call

    'positive numbers only, modulus must be >1!  Find mod inverse if s2=-1.
    out$ = errormsg$
    IF LEN(s3$) THEN n$ = s3$
    IF bIsNeg%(s1$) OR bIsNeg%(n$) THEN EXIT SUB
    IF bIsNeg%(s2$) THEN
        IF bIsEqual%(s2$, "-1") THEN bModInv s1$, n$, out$
        EXIT SUB
    END IF

    c$ = s1$
    d$ = s2$
    out$ = one$

    DO
        IF bIsOdd%(d$) THEN
            z$ = out$
            bMul z$, c$, out$
            z$ = out$
            bMod z$, n$, out$
        END IF
        z$ = c$
        w$ = c$
        bMul z$, w$, c$
        z$ = c$
        bMod z$, n$, c$
        z$ = d$
        bDivInt z$, two$, d$
    LOOP UNTIL bIsZero%(d$)

END SUB

SUB bModInv (s1$, s2$, out$)
    DIM g0$, g1$, g2$, v0$, v1$, v2$, y$, t$, z$

    IF NOT bIsRelPrime%(s1$, s2$) THEN out$ = zero$: EXIT SUB

    g0$ = s2$: g1$ = s1$
    v0$ = zero$: v1$ = one$

    DO UNTIL bIsZero%(g1$)
        bDivInt g0$, g1$, y$
        bMul y$, g1$, t$
        bSub g0$, t$, g2$
        bMul y$, v1$, t$
        bSub v0$, t$, v2$
        g0$ = g1$: g1$ = g2$
        v0$ = v1$: v1$ = v2$
    LOOP

    out$ = v0$
    IF bIsNeg%(out$) THEN
        z$ = out$
        bAdd z$, s2$, out$
    END IF
END SUB

SUB bGCD (s1$, s2$, out$)
    DIM div$, dvd$, t$

    'work with copies
    div$ = s1$
    dvd$ = s2$
    IF bIsMore%(div$, dvd$) THEN bSwapString div$, dvd$

    DO UNTIL bIsZero%(div$)
        bMod dvd$, div$, t$
        dvd$ = div$
        div$ = t$
    LOOP
    out$ = dvd$

END SUB


SUB bSqrInt (s$, out$)
    DIM t$
    DIM olddigits%

    IF bIsNeg%(s$) THEN out$ = errormsg$: EXIT SUB
    t$ = s$
    bInt t$

    'a trick: let bSqr() figure the decimal and only find that many digits
    olddigits% = digits%
    digits% = 0
    bSqr t$, out$
    digits% = olddigits%

END SUB

FUNCTION bIsBase% (s$)
    bIsBase% = INSTR(UCASE$(s$), basechr$)
END FUNCTION

'return true if s1 divides s2
'
FUNCTION bIsDiv% (s1$, s2$)
    DIM t$

    bMod s2$, s1$, t$
    bIsDiv% = (t$ = zero$)
END FUNCTION

'return true if s1 = s2
'
FUNCTION bIsEqual% (s1$, s2$)
    bIsEqual% = (s1$ = s2$)
END FUNCTION

'return true if s$ is even, no decimals!
'
FUNCTION bIsEven% (s$)
    bIsEven% = (VAL(RIGHT$(s$, 1)) MOD 2 = 0)
END FUNCTION

'return true if s in an integer (no decimal point).
'
FUNCTION bIsInteger% (s$)
    bIsInteger% = (INSTR(s$, dp$) = 0)
END FUNCTION

'return true if s1 < s2
'
FUNCTION bIsLess% (s1$, s2$)
    bIsLess% = (bComp%(s1$, s2$) = -1)
END FUNCTION

FUNCTION bComp% (s1$, s2$)
    DIM s1flag%, s2flag%, sign1%, sign2%
    DIM dp1%, dp2%, arg%

    'kludge to fix 0<.1
    IF LEFT$(s1$, 1) = dp$ THEN s1$ = zero$ + s1$: s1flag% = True
    IF LEFT$(s2$, 1) = dp$ THEN s2$ = zero$ + s2$: s2flag% = True

    sign1% = (LEFT$(s1$, 1) = neg$)
    sign2% = (LEFT$(s2$, 1) = neg$)
    dp1% = INSTR(s1$, dp$): IF dp1% = 0 THEN dp1% = LEN(s1$) + 1
    dp2% = INSTR(s2$, dp$): IF dp2% = 0 THEN dp2% = LEN(s2$) + 1

    IF sign1% <> sign2% THEN
        IF sign1% THEN arg% = -1 ELSE arg% = 1
    ELSEIF s1$ = s2$ THEN
        arg% = 0
    ELSEIF (dp1% < dp2%) OR ((dp1% = dp2%) AND (s1$ < s2$)) THEN
        arg% = -1
    ELSE
        arg% = 1
    END IF

    IF sign1% AND sign2% THEN arg% = -arg%
    IF s1flag% THEN s1$ = MID$(s1$, 2)
    IF s2flag% THEN s2$ = MID$(s2$, 2)
    bComp% = arg%

END FUNCTION

'return true if s1 > s2
'
FUNCTION bIsMore% (s1$, s2$)
    bIsMore% = (bComp%(s1$, s2$) = 1)
END FUNCTION

'return true if s is negative
'
FUNCTION bIsNeg% (s$)
    bIsNeg% = (LEFT$(s$, 1) = neg$)
END FUNCTION

FUNCTION bIsNotZero% (s$)
    DIM flag%, i%

    flag% = False
    FOR i% = 1 TO LEN(s$)
        IF INSTR("0-. ", MID$(s$, i%, 1)) = False THEN flag% = True: EXIT FOR
    NEXT i%
    bIsNotZero% = flag%
END FUNCTION

'return true if odd
'
FUNCTION bIsOdd% (s$)
    bIsOdd% = (VAL(RIGHT$(s$, 1)) MOD 2 <> 0)
END FUNCTION

'return true if s is prime
'
FUNCTION bIsPrime% (s$)
    bIsPrime% = (bPrmDiv$(s$, False) = s$)
END FUNCTION

's is pseudoprime to base b if (b,s)=1 and b^(s-1)=1 (mod s).  Integers only!
'
FUNCTION bIsPseudoPrime% (s$, bas$)
    DIM t$, smin$
    DIM flag%

    flag% = False
    IF bIsRelPrime%(s$, bas$) THEN
        smin$ = s$: bInc smin$, -1
        bModPower bas$, smin$, s$, t$
        flag% = (t$ = one$)
    END IF
    bIsPseudoPrime% = flag%
END FUNCTION

'return true if s1 and s2 are relatively prime, ie share no factor
'
FUNCTION bIsRelPrime% (s1$, s2$)
    DIM gcd$

    bGCD s1$, s2$, gcd$
    bIsRelPrime% = bIsEqual%(gcd$, one$)
END FUNCTION

'Return true if s$ is zero$ or null, s$ needn't be clean.
'
FUNCTION bIsZero% (s$)
    DIM flag%, i%

    flag% = True
    FOR i% = 1 TO LEN(s$)
        IF INSTR("0-. ", MID$(s$, i%, 1)) = False THEN flag% = False: EXIT FOR
    NEXT i%
    bIsZero% = flag%
END FUNCTION

FUNCTION bSign% (s$)
    IF bIsNeg%(s$) THEN bSign% = negative ELSE bSign% = positive
END FUNCTION

FUNCTION bLogDp% (s$, logdp%)
    bLogDp% = LEN(s$) - 1 - logdp%
END FUNCTION

FUNCTION bPrmDiv$ (s$, dspflag%)
    DIM num$, sfac$, maxfac$, t$
    DIM lfac&, lnum&, lmaxfac&, ldfac&
    DIM i%, cnt%, flag%, dfac%

    num$ = s$
    bInt num$
    bAbs num$
    IF LEN(num$) <= maxlongdig THEN GOSUB bpdLong ELSE GOSUB bpdChar
    EXIT FUNCTION

    bpdChar:
    'try some classic divisibility tests for small factors.
    'Cf Gardner, Unexpected Hanging, p.160.

    'by 2?
    '  If dspflag% Then
    '  frmBncFactor.lblTryNum.Caption = two$
    '  frmBncFactor.lblTryNum.Refresh
    'End If
    IF VAL(RIGHT$(num$, 1)) MOD 2 = 0 THEN bPrmDiv$ = two$: RETURN

    'by 3?
    'IF dspflag% THEN LOCATE , dspflag%: PRINT three$;
    '  If dspflag% Then
    '  frmBncFactor.lblTryNum.Caption = three$
    '  frmBncFactor.lblTryNum.Refresh
    'End If

    lfac& = 0
    FOR i% = 1 TO LEN(num$)
        lfac& = lfac& + ASC(MID$(num$, i%, 1)) - asc0
    NEXT i%
    IF lfac& MOD 3 = 0 THEN bPrmDiv$ = three$: RETURN

    'by 5?
    'IF dspcol% THEN LOCATE , dspcol%: PRINT five$;
    '  If dspflag% Then
    '  frmBncFactor.lblTryNum.Caption = five$
    '  frmBncFactor.lblTryNum.Refresh
    'End If

    IF VAL(RIGHT$(num$, 1)) MOD 5 = 0 THEN bPrmDiv$ = five$: RETURN

    'by 7, 11, or 13?
    'IF dspcol% THEN LOCATE , dspcol%: PRINT "7+";
    '  If dspflag% Then
    '  frmBncFactor.lblTryNum.Caption = "7+"
    '  frmBncFactor.lblTryNum.Refresh
    'End If

    lfac& = 0
    i% = LEN(num$) + 1
    cnt% = 3
    flag% = True
    DO
        i% = i% - 3: IF i% < 1 THEN cnt% = i% + 2: i% = 1
        IF flag% THEN
            lfac& = lfac& + VAL(MID$(num$, i%, cnt%))
        ELSE
            lfac& = lfac& - VAL(MID$(num$, i%, cnt%))
        END IF
        flag% = NOT flag%
    LOOP WHILE i% > 1
    IF lfac& MOD 7 = 0 THEN bPrmDiv$ = "7": RETURN
    IF lfac& MOD 11 = 0 THEN bPrmDiv$ = "11": RETURN
    IF lfac& MOD 13 = 0 THEN bPrmDiv$ = "13": RETURN

    'main loop, increment factor by 2 or 4.
    sfac$ = "17"
    dfac% = 2
    bSqrInt num$, maxfac$

    DO
        'IF dspcol% THEN LOCATE , dspcol%: PRINT sfac$;
        '    If dspflag% Then
        '  frmBncFactor.lblTryNum.Caption = sfac$
        '  frmBncFactor.lblTryNum.Refresh
        'End If

        bMod num$, sfac$, t$
        IF bIsZero%(t$) THEN EXIT DO
        bInc sfac$, dfac%
        dfac% = 6 - dfac%
        IF bIsMore%(sfac$, maxfac$) THEN sfac$ = num$: EXIT DO
        'If INKEY$ = esc$ Then sfac$ = zero$: Exit Do
    LOOP
    bPrmDiv$ = sfac$
    RETURN

    bpdLong:
    lnum& = VAL(num$)
    IF lnum& <= 1 THEN
        lfac& = 1&
    ELSEIF lnum& MOD 2& = 0& THEN
        lfac& = 2&
    ELSEIF lnum& MOD 3& = 0& THEN
        lfac& = 3&
    ELSE
        lmaxfac& = INT(SQR(lnum&))
        lfac& = 5&
        ldfac& = 2&
        DO
            'IF dspcol% THEN LOCATE , dspcol%: PRINT lfac&;
            '      If dspflag% Then
            '  frmBncFactor.lblTryNum.Caption = LTrim$(Str$(lfac&))
            '  frmBncFactor.lblTryNum.Refresh
            'End If

            IF lnum& MOD lfac& = 0& THEN EXIT DO
            lfac& = lfac& + ldfac&
            ldfac& = 6& - ldfac&
            IF lfac& > lmaxfac& THEN lfac& = lnum&: EXIT DO
        LOOP
    END IF
    bPrmDiv$ = LTRIM$(STR$(lfac&))
    RETURN

END FUNCTION

FUNCTION bMaxInt% (n1%, n2%)
    IF n1% >= n2% THEN bMaxInt% = n1% ELSE bMaxInt% = n2%
END FUNCTION

@Jack

+2 to you, and if anyone wants Treebeard's string math just for +-*/ and square root, this is the stripped down routine, which provides that.

Pete
Reply




Users browsing this thread: 1 Guest(s)