calculate pi to different bases

[jack]

converted from a post in the freebasic forum
example output

Code: [Select]

pi in base 10 = 3.1415926535897932384626433832795028841971693993751
pi in base 16 = 3.243F6A8885A308D313198A2E03707344A4093822299F31D01
pi in base 2 = 11.0010010000111111011010101000100010000101101001000


Code: [Select]

Declare Function find_pi$(Byval n As Integer, Byval nbase As Integer)

PRINT "waiting"
PRINT find_pi$(50, 10)
PRINT find_pi$(50, 16)
PRINT find_pi$(50, 2)

FUNCTION convert$ (n AS _INTEGER64, nbase AS INTEGER)
    IF nbase < 2 OR nbase > 16 THEN
        convert$ = "base must be between 2 and 16"
        EXIT FUNCTION
    END IF
    DIM sn AS STRING
    DIM r AS _INTEGER64
    DIM h(15) AS STRING * 1
    h(0) = "0": h(1) = "1": h(2) = "2": h(3) = "3": h(4) = "4": h(5) = "5": h(6) = "6": h(7) = "7": h(8) = "8"
    h(9) = "9": h(10) = "A": h(11) = "B": h(12) = "C": h(13) = "D": h(14) = "E": h(15) = "F"
    IF n = 0 THEN
        convert$ = "0"
        EXIT FUNCTION
    ELSEIF n < 0 THEN
        convert$ = "positive number only"
        EXIT FUNCTION
    END IF
    sn = ""
    WHILE n >= nbase
        r = n MOD nbase
        n = n \ nbase
        sn = h(r) + sn
    WEND
    sn = h(n) + sn
    convert$ = sn
END FUNCTION

SUB array_mul (a() AS _INTEGER64, f AS INTEGER, n AS INTEGER, nbase AS INTEGER)
    DIM c AS _INTEGER64, p AS _INTEGER64, h AS _INTEGER64
    c = 0
    FOR p = n TO 0 STEP -1
        h = a(p) * f + c
        c = h \ nbase
        a(p) = h MOD nbase
    NEXT p
END SUB

SUB array_div (a() AS _INTEGER64, f AS INTEGER, n AS INTEGER, nbase AS INTEGER)
    DIM b AS _INTEGER64, p AS _INTEGER64, h AS _INTEGER64
    b = 0
    FOR p = 0 TO n
        h = a(p) + b * nbase
        a(p) = h \ f
        b = h MOD f
    NEXT p
END SUB

SUB array_add (a() AS _INTEGER64, b() AS _INTEGER64, n AS INTEGER, nbase AS INTEGER)
    DIM c AS _INTEGER64, p AS _INTEGER64, h AS _INTEGER64
    c = 0
    FOR p = n TO 0 STEP -1
        h = a(p) + b(p) + c
        c = h \ nbase
        a(p) = h MOD nbase
    NEXT p
END SUB

SUB array_sub (a() AS _INTEGER64, b() AS _INTEGER64, n AS INTEGER, nbase AS INTEGER)
    DIM c AS _INTEGER64, p AS _INTEGER64, h AS _INTEGER64
    c = 0
    FOR p = n TO 0 STEP -1
        h = a(p) - b(p) + nbase
        c = h \ nbase
        a(p) = h MOD nbase
        IF (c = 0) THEN a(p - 1) = a(p - 1) - 1
    NEXT p
END SUB

FUNCTION is_zero% (a() AS _INTEGER64, n AS INTEGER)
    DIM p AS INTEGER

    FOR p = 0 TO n
        IF a(p) <> 0 THEN
            is_zero% = 0
            EXIT FUNCTION
        END IF
    NEXT p
    is_zero% = 1
END FUNCTION

SUB arctan (t() AS _INTEGER64, s() AS _INTEGER64, div AS INTEGER, n AS INTEGER, nbase AS INTEGER)
    DIM w AS _INTEGER64, i AS INTEGER
    s(0) = 1
    i = 1
    w = div
    array_div s(), w, n, nbase
    array_add t(), s(), n, nbase
    DO
        array_mul s(), i, n, nbase
        w = div * div
        array_div s(), w, n, nbase
        i = i + 2
        w = i
        array_div s(), w, n, nbase
        array_sub t(), s(), n, nbase
        array_mul s(), i, n, nbase
        w = div * div
        array_div s(), w, n, nbase
        i = i + 2
        w = i
        array_div s(), w, n, nbase
        array_add t(), s(), n, nbase
    LOOP UNTIL is_zero%(s(), n) = 1
END SUB

FUNCTION find_pi$ (n AS INTEGER, nbase AS INTEGER)
    IF nbase < 2 OR nbase > 16 THEN
        find_pi$ = "base must be between 2 and 16"
        EXIT FUNCTION
    END IF
    IF n < 2 THEN n = 2
    'DIM tm AS DOUBLE: tm = TIMER
    DIM digit AS STRING, uit AS STRING
    uit = "pi in base " + STR$(nbase) + " = " + convert(3, nbase) + "."
    DIM s(n + 2) AS _INTEGER64
    DIM t(n + 2) AS _INTEGER64
    DIM i AS INTEGER, f AS INTEGER
    f = 2
    arctan t(), s(), f, n, nbase
    f = 3
    arctan t(), s(), f, n, nbase
    array_mul t(), 4, n, nbase
    'tm = TIMER - tm
    FOR i = 1 TO n - 1
        digit = convert(t(i), nbase)
        uit = uit + digit
    NEXT i
    'PRINT "elapsed time = "; tm; " seconds"
    find_pi$ = uit
END FUNCTION


[bplus]

Oh man base 2, so cool!

They are doing this at JB forum (also inspired by FB) and are a bit stuck in base 10.

So how do you check if results are correct, say in 100th digit?

[jack]

hi bplus
there's no check, but from trial runs and comparison with, for example pari/gp, the result is accurate, the last digit may be rounded-up
can you provide a link to the discussion?

[bplus]

Well I did manage to check a base 10 result for some algo by comparing digits to a PI posted over Internet, turned out that method was correct only to like the first 75% of digits, so if you wanted 100 accurate digits you had to go past 4/3 past 100 ~ 133 places.

Can you translate a base n result to a base 10 accurately? because then you can check the base 10 result.

[jack]

I am not sure I understand you correctly, the calculations are done in the number base selected, there's no conversion other than converting the digits to the appropriate base, for example, in base 16, 10 is mapped to A and so on
for base 10 and digits 100, I get "pi in base 10 = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117070"
with an error of -2.0178519134867177149 E-99 because the result is rounded up
pari/gp gives Pi = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

[bplus]

I am not sure I understand you correctly, the calculations are done in the number base selected,

Yes, and then I am saying you could check a base 2 calculation if you could accurately convert a base 2 expansion to a base 10 expansion, then you could just see if the base 10 conversion was accurate to a reference source. True though you couldn't verify digit by digit, you could verify the whole thing were matched x amount of digits base 10 which would then imply the digits up to some other point were correct.

I don't know, maybe an algo to convert a base x expansion to a base y expansion might be more complicated than the original pi calculation? though it might be more useful! ;-)) ie computers are converting binary to base 10 all the time.

[jack]

I think that to verify the result you would need to calculate pi using two different methods and then compare the results, the program I posted calculates pi in the given base, no conversion is necessary except to display the result.
let me explain, you could instead of using the convert function, if the base is between 4 and 10 then simply use str$ to convert the digits to string, if the base is 2 then the first two digits are 11, if the base is 3, 10 thereafter just use str$
if the base is greater than 10 you could use hex$ instead of str$

[bplus]

To verify, one could also compare a result to a well respected, well researched, peer reviewed reference, which is what I am talking about. And to do that I might have to convert a base n result to a base 10 result... and if that can be done accurately then I can have verification.

Of course, if they use same method or algo to find pi to x amount of digits.... not a great verification but at least I have respectable company from whom we have agreement. :)

Update:

Our posts have crisscrossed, you have edited and started picking up on what I mean, conversion, better word than translate I suppose.

[jack]

one reason for using _integer64 was to have the option to increase the number of digits per array element, but I forgot to declare nbase as _integer64
you can have upto 7 digits per element, so you divide n by 7 and add 1, the problem is that you are not alway going to get n digits, so we add 2 instead
the changes are in the function find_pi$
sample output

Code: [Select]

pi in base 10 = 3.14159265358979323846264338327950288419716939937510582097
pi in base 16 = 3.243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98
pi in base 2 = 11.00100100001111110110101010001000100001011010001100001000


Code: [Select]

Declare Function find_pi$(Byval n As Integer, Byval nbase As _Integer64)

PRINT "waiting"
PRINT find_pi$(50, 10)
PRINT find_pi$(50, 16)
PRINT find_pi$(50, 2)

FUNCTION convert$ (n AS _INTEGER64, nbase AS _INTEGER64)
    IF nbase < 2 OR nbase > 16 THEN
        convert$ = "base must be between 2 and 16"
        EXIT FUNCTION
    END IF
    DIM sn AS STRING
    DIM r AS _INTEGER64
    DIM h(15) AS STRING * 1
    h(0) = "0": h(1) = "1": h(2) = "2": h(3) = "3": h(4) = "4": h(5) = "5": h(6) = "6": h(7) = "7": h(8) = "8"
    h(9) = "9": h(10) = "A": h(11) = "B": h(12) = "C": h(13) = "D": h(14) = "E": h(15) = "F"
    IF n = 0 THEN
        convert$ = "0"
        EXIT FUNCTION
    ELSEIF n < 0 THEN
        convert$ = "positive number only"
        EXIT FUNCTION
    END IF
    sn = ""
    WHILE n >= nbase
        r = n MOD nbase
        n = n \ nbase
        sn = h(r) + sn
    WEND
    sn = h(n) + sn
    convert$ = sn
END FUNCTION

SUB array_mul (a() AS _INTEGER64, f AS INTEGER, n AS INTEGER, nbase AS _INTEGER64)
    DIM c AS _INTEGER64, p AS _INTEGER64, h AS _INTEGER64
    c = 0
    FOR p = n TO 0 STEP -1
        h = a(p) * f + c
        c = h \ nbase
        a(p) = h MOD nbase
    NEXT p
END SUB

SUB array_div (a() AS _INTEGER64, f AS INTEGER, n AS INTEGER, nbase AS _INTEGER64)
    DIM b AS _INTEGER64, p AS _INTEGER64, h AS _INTEGER64
    b = 0
    FOR p = 0 TO n
        h = a(p) + b * nbase
        a(p) = h \ f
        b = h MOD f
    NEXT p
END SUB

SUB array_add (a() AS _INTEGER64, b() AS _INTEGER64, n AS INTEGER, nbase AS _INTEGER64)
    DIM c AS _INTEGER64, p AS _INTEGER64, h AS _INTEGER64
    c = 0
    FOR p = n TO 0 STEP -1
        h = a(p) + b(p) + c
        c = h \ nbase
        a(p) = h MOD nbase
    NEXT p
END SUB

SUB array_sub (a() AS _INTEGER64, b() AS _INTEGER64, n AS INTEGER, nbase AS _INTEGER64)
    DIM c AS _INTEGER64, p AS _INTEGER64, h AS _INTEGER64
    c = 0
    FOR p = n TO 0 STEP -1
        h = a(p) - b(p) + nbase
        c = h \ nbase
        a(p) = h MOD nbase
        IF (c = 0) THEN a(p - 1) = a(p - 1) - 1
    NEXT p
END SUB

FUNCTION is_zero% (a() AS _INTEGER64, n AS INTEGER)
    DIM p AS INTEGER

    FOR p = 0 TO n
        IF a(p) <> 0 THEN
            is_zero% = 0
            EXIT FUNCTION
        END IF
    NEXT p
    is_zero% = 1
END FUNCTION

SUB arctan (t() AS _INTEGER64, s() AS _INTEGER64, div AS INTEGER, n AS INTEGER, nbase AS _INTEGER64)
    DIM w AS _INTEGER64, i AS INTEGER
    s(0) = 1
    i = 1
    w = div
    array_div s(), w, n, nbase
    array_add t(), s(), n, nbase
    DO
        array_mul s(), i, n, nbase
        w = div * div
        array_div s(), w, n, nbase
        i = i + 2
        w = i
        array_div s(), w, n, nbase
        array_sub t(), s(), n, nbase
        array_mul s(), i, n, nbase
        w = div * div
        array_div s(), w, n, nbase
        i = i + 2
        w = i
        array_div s(), w, n, nbase
        array_add t(), s(), n, nbase
    LOOP UNTIL is_zero%(s(), n) = 1
END SUB

FUNCTION find_pi$ (n AS INTEGER, nbase AS _INTEGER64)
    IF nbase < 2 OR nbase > 16 THEN
        find_pi$ = "base must be between 2 and 16"
        EXIT FUNCTION
    END IF
    IF n < 2 THEN n = 2
    n = n / 7 + 2
    'DIM tm AS DOUBLE: tm = TIMER
    DIM digit AS STRING, uit AS STRING
    uit = "pi in base " + STR$(nbase) + " = " + convert(3, nbase) + "."
    DIM s(n + 2) AS _INTEGER64
    DIM t(n + 2) AS _INTEGER64
    DIM i AS INTEGER, f AS INTEGER
    f = 2
    arctan t(), s(), f, n, nbase * nbase * nbase * nbase * nbase * nbase * nbase
    f = 3
    arctan t(), s(), f, n, nbase * nbase * nbase * nbase * nbase * nbase * nbase
    array_mul t(), 4, n, nbase * nbase * nbase * nbase * nbase * nbase * nbase
    'tm = TIMER - tm
    FOR i = 1 TO n - 1
        digit = convert(t(i), nbase)
        IF LEN(digit) < 7 THEN digit = STRING$(7 - LEN(digit), "0") + digit
        uit = uit + digit
    NEXT i
    'PRINT "elapsed time = "; tm; " seconds"
    find_pi$ = uit
END FUNCTION
in case you want to check the accuracy https://pi2e.ch/blog/2017/03/10/pi-digits-download/#download