A quick google search reveals that in C, it's possible to do this algorithm with the GMP library, which FB does have a .bi include file for. Should be trivial to convert the code in the example to FB (1:1 transliteration) contained in this article:
http://beej.us/blog/data/pi-chudnovsky-gmp/
GMP supports arbitrary precision numbers, and is optimized to be as fast as possible, including likely using the tricks mentioned by jj2007.
EDIT: Here's a port of the C program in that web site:
Code: Select all
#include "gmp.bi"
'how many to display if the user doesn't specify:
#define DEFAULT_DIGITS 60
'how many decimal digits the algorithm generates per iteration:
#define DIGITS_PER_ITERATION 14.1816474627254776555
'/**
' * Compute pi to the specified number of decimal digits using the
' * Chudnovsky Algorithm.
' *
' * http://en.wikipedia.org/wiki/Pi#Rapidly_convergent_series
' *
' * NOTE: this function returns a malloc()'d string!
' *
' * @param digits number of decimal digits to compute
' *
' * @return a malloc'd string result (with no decimal marker)
' */
type mpf_t as __mpf_struct
type mpz_t as __mpz_struct
function chudnovsky(digits as culong) as zstring ptr
dim as mpf_t con, AA, BB, FF, sum
dim as mpz_t a, b, c, d, e
dim as zstring ptr calc_output
dim as mp_exp_t exp_
dim as double bits_per_digit
dim as culong k, threek
dim as culong iterations
dim as culong precision_bits
iterations = digits / DIGITS_PER_ITERATION + 1
'roughly compute how many bits of precision we need for
'this many digit:
bits_per_digit = 3.32192809488736234789 'log2(10)
precision_bits = (digits * bits_per_digit) + 1
mpf_set_default_prec(precision_bits)
'allocate GMP variables
mpf_inits(@con, @AA, @BB, @FF, @sum, NULL)
mpz_inits(@a, @b, @c, @d, @e, NULL)
mpf_set_ui(@sum, 0) 'sum already zero at this point, so just FYI
'first the constant sqrt part
mpf_sqrt_ui(@con, 10005)
mpf_mul_ui(@con, @con, 426880)
' now the fun bit
for k = 0 to iterations-1
threek = 3*k
mpz_fac_ui(@a, 6*k) '(6k)!
mpz_set_ui(@b, 545140134) '13591409 + 545140134k
mpz_mul_ui(@b, @b, k)
mpz_add_ui(@b, @b, 13591409)
mpz_fac_ui(@c, threek) '(3k)!
mpz_fac_ui(@d, k) '(k!)^3
mpz_pow_ui(@d, @d, 3)
mpz_ui_pow_ui(@e, 640320, threek) ' -640320^(3k)
if (threek and 1) = 1 then mpz_neg(@e, @e)
'numerator (in A)
mpz_mul(@a, @a, @b)
mpf_set_z(@AA, @a)
'denominator (in B)
mpz_mul(@c, @c, @d)
mpz_mul(@c, @c, @e)
mpf_set_z(@BB, @c)
'result
mpf_div(@FF, @AA, @BB)
'add on to sum
mpf_add(@sum, @sum, @FF)
next
'final calculations (solve for pi)
mpf_ui_div(@sum, 1, @sum) ' invert result
mpf_mul(@sum, @sum, @con) ' multiply by constant sqrt part
'get result base-10 in a string:
calc_output = mpf_get_str(NULL, @exp_, 10, digits, @sum) 'calls malloc()
' free GMP variables
mpf_clears(@con, @AA, @BB, @FF, @sum, NULL)
mpz_clears(@a, @b, @c, @d, @e, NULL)
return calc_output
end function
dim as zstring ptr result
if len(command(1)) then
result = chudnovsky(val(command(1)))
else
result = chudnovsky(DEFAULT_DIGITS)
end if
'add a decimal after the first digit
print mid(result[0],1,1);".";result[1]
Deallocate(result)
The types are a little odd in that code, such as culong. That's just because I was doing a 1:1 translation, and also because the gmp calls were using unsigned ints in C, which is culong in FB.