Getting past the 15 decimal place limit of a double in FreeBASIC?

[blackburnfjames]

Currently I am attempting to use the Chudnovsky Algorithm within FreeBASIC. I am using the slightly more dated formula: https://wikimedia.org/api/rest_v1/media ... ccb1519887 yet i am encountering the limitation of a 15 d.p. double. The formula I am using in FreeBASIC for a single generation:

Code: Select all

 1/((((factorial(6*k)) * (13591409 + (545140134*k)))/((factorial(3*k))*((factorial(k))^3)*((-640320)^(3*k))))*(12/(640320^(3/2))))

Yet when i place this into a variable. For example my full code is:

Code: Select all

function factorial( byval i as integer ) as integer
    var answer = 1
    while( i > 1 )
        answer *= i
        i -= 1
    wend
    return answer
end function
dim pi as double
dim k as integer

k = 0
pi = 1/((((factorial(6*k)) * (13591409 + (545140134*k)))/((factorial(3*k))*((factorial(k))^3)*((-640320)^(3*k))))*(12/(640320^(3/2)))) 

When I run this i get:

Code: Select all

3.141592653589734

This isn't exactly enough digits when playing with Pi. Is there anyway to bypass this limitation and increase the digit amount?
Sorry if this is in the wrong area I'm new to this forum please redirect me to the correct area if so. Also if this already has an answer elsewhere please redirect me to where.
Thanks in advance,
James

[jj2007]

There is no simple way to overcome the limits of a double. REAL8 precision is 15 significant digits. You can add 3-4 digits by using the FPU, though: FB has built-in assembler support. You could use the FPU for the non-integer parts of your formula, but it requires some knowledge of assembler. Simple example:

Code: Select all

dim tmp as integer 
k=5
tmp=factorial(6*k)
asm
  int 3
  fild dword ptr [tmp]	; ST0 = 1409286144.0000000000 at REAL10 precision
end asm

[caseih]

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.

[jj2007]

FreeBasic\bin\win32\ld.exe: cannot find -lgmp

[Josep Roca]

> FreeBasic\bin\win32\ld.exe: cannot find -lgmp

When trying to use third party libraries, specially the open source ones, you have to surf the web to find them.

See: https://gmplib.org/

Be aware that the GNU licenses set firm restrictions on the use with non-free programs.

[dodicat]

srvaldez supplied a gmp 6 a while back.
http://www.mediafire.com/file/6hj8ry1jj ... -6.1.1.zip

(The freebasic .bi file is OK to use with it)

Also an older 32 bit gmp
https://sourceforge.net/projects/fbc/fi ... Libraries/

[caseih]

Sure. All third-party libraries, be they open source or proprietary, have licenses that must be abided by and establish firm restrictions. Some are royalty-free, some are not. Open source is not special in this regard. GMP is available under the LGPL license, which does allow for use in proprietary software, provided it's dynamically linked. So no reason to fear the GNU! :) Also the GNU licenses are unique compared to proprietary licenses in that you need not agree with the license to simply use the GPL or LGPL code. The license only comes into effect when you distribute your derivative work to others. So for messing around and trying things, in this case the license does not even matter.

I can't think of any application of the chudonovsky algortithm in a commercial program (hard-coding the value of pi would be the most optimal), so I expect that he is exploring it for fun, so the license of GMP would be nearly irrelevant. But always good to be aware of license terms.

[caseih]

FreeBasic\bin\win32\ld.exe: cannot find -lgmp

I always try to get binaries of third-party open-source libraries from the MingW compiler project project if possible. Here's a link to the mingw page on libgmp:
https://sourceforge.net/projects/mingw/ ... gmp-6.1.2/
The dll itself is in a tarball at https://sourceforge.net/projects/mingw/ ... z/download

A good archiver program like 7-zip can open tar.xz files and let you extract the dll.

EDIT: sigh, Windows doesn't make anything easy, particularly when it comes to versioning dlls. The extracted dll will have to be renamed to libgmp.dll (or just gmp.dll). Might be simpler to just use the Windows Subsystem for Linux where everything can be pulled in with a simple package manager install command! But I digress.

[dodicat]

Much easier to use a static GMP lib.
Here is a pi I dug up.

Code: Select all

 Function pie(n As Long) As String
    #define AddaDigit(d)  pi[inc]=d:inc+=1
    Dim As Long  _len = 20*n\4,inc,k,x,m,nines,digit,q
    Dim As String pi=String(n,0)
    Redim As Long a(_len)
    For j As Long  = 1 To _len:a(j-1) = 2:Next
        For j As Long = 1 To n
            For i As Long = _len To 1 Step -1
                x   = 10*a(i-1) + q*i
                m=(i Shl 1 - 1)
                a(i-1) = x Mod m
                q= x\m
            Next i
            a(0) = q Mod 10
            q =q \ 10
            If q = 9 Then
                nines += 1
            Else
                If q = 10 Then
                    AddaDigit((Digit+49))
                    If nines > 0 Then
                        For k = 1 To nines:AddaDigit(48):Next
                        End If
                        Digit = 0
                        nines = 0
                    Else
                        AddaDigit((Digit+48))
                        Digit = q
                        If nines <> 0 Then
                            For k = 1 To nines:AddaDigit(57):Next
                                nines = 0
                            End If
                        End If
                    End If
                Next j
                Return ("3."+Ltrim(pi,"03"))
            End Function
            
            
            Print pie(1000)
            Sleep 

[St_W]

Much easier to use a static GMP lib.

Just note the implications to your piece of software in case you decide to share it. Basically it means that you also have to use the LGPL license for your code and thus have to share it (not only binaries) if you link GMP statically. And again, this becomes only relevant as soon as you share your code.

[srvaldez]

When I run this i get:

Code: Select all

3.141592653589734

This isn't exactly enough digits when playing with Pi. Is there anyway to bypass this limitation and increase the digit amount?
Sorry if this is in the wrong area I'm new to this forum please redirect me to the correct area if so. Also if this already has an answer elsewhere please redirect me to where.
Thanks in advance,
James

in your program you only use the first term of the summation, the value you quote is correct when using only one term
to get better accuracy use two terms

Code: Select all

function factorial( byval i as integer ) as integer
    var answer = 1
    while( i > 1 )
        answer *= i
        i -= 1
    wend
    return answer
end function
dim pi as double
dim k as integer

k = 0
pi = (factorial(6*k) * (13591409 + (545140134)*k))/(factorial(3*k)*factorial(k)^3*(-640320)^(3*k))
k += 1
pi += (factorial(6*k) * (13591409 + (545140134)*k))/(factorial(3*k)*factorial(k)^3*(-640320)^(3*k))
pi*=(12/(640320^(3/2)))
print 1/pi

[bluatigro]

i was trying to to do this
but then whit a choise about how many digit's
viewtopic.php?f=3&t=26917
i tryed first the integer math to work
it is not jet good