C Multiple Precision Variables in FreeBSD

I've been looking for a way to do this for several years [albeit not very hard]. Anyways, I saw a great example at mpfr.org. The libraries were already installed in my FreeBSD version [others may need to install them]. This example calculates Euler's number to an arbitrary degree of precision.

[The order of the header files is important, mpfr.h uses functions from gmp.h and stdio.h]


C:
//This example of multiple precision variables was taken from: https://www.mpfr.org/sample.html
//Make sure the gmp and mpfr libraries are installed
///header files should be in: "/usr/local/include/mpfr.h" and "/usr/local/include/gmp.h"
//To compile:
//gcc -o sample sample.c -lmpfr

#include <stdio.h>

#include <gmp.h>
#include <mpfr.h>

int main (void)
{
  unsigned int i;
  mpfr_t s, t, u;   //declares three floating-point variables s, t, u

  mpfr_init2 (t, 200);  //initializes the variable t with a 200-bit precision
  mpfr_set_d (t, 1.0, MPFR_RNDD);   //sets the value of t to the double-precision number 1.0 rounded toward minus infinity (here no rounding is involved since 1 is represented exactly as a double-precision number and also as a 200-bit MPFR number)
  mpfr_init2 (s, 200);
  mpfr_set_d (s, 1.0, MPFR_RNDD);
  mpfr_init2 (u, 200);
  for (i = 1; i <= 100; i++)
    {
      mpfr_mul_ui (t, t, i, MPFR_RNDU); //multiplies t in place by the unsigned integer i, where the result is rounded toward plus infinity
      mpfr_set_d (u, 1.0, MPFR_RNDD);
      mpfr_div (u, u, t, MPFR_RNDD);    //divides u by t, rounding the result toward minus infinity, and stores it in u
      mpfr_add (s, s, u, MPFR_RNDD);
    }
  printf ("Sum is ");
  mpfr_out_str (stdout, 10, 0, s, MPFR_RNDD); //prints the value of s in base 10, rounded toward minus infinity, where the third argument 0 means that the number of printed digits is automatically chosen from the precision of s (note: mpfr_printf could also be used instead of printf, mpfr_out_str and putchar)
  putchar ('\n');
  mpfr_clear (s); // The mpfr_clear and mpfr_free_cache calls free the space used by the MPFR variables and caches
  mpfr_clear (t);
  mpfr_clear (u);
  mpfr_free_cache ();
  return 0;
}
 
So with this I made a really precise estimate of Euler's number:

Code:
Euler's Number: 2.718281828459045235360287471352662497757247093699959574966967627724
076630353547594571382178525166427427466391932003059921817413596629043572900334295
260595630738132328627943490763233829880753195251019011573834187930702154089149934
884167509244761460668082264800168477411853742345442437107539077744992069551702761
838606261331384583000752044933826560297606737113200709328709127443747047230696977
209310141692836819025515108657463772111252389784425056953696770785449969967946864
454905987931636889230098793127736178215424999229576351482208269895193668033182528
869398496465105820939239829488793320362509443117301238197068416140397019837679320
683282376464804295311802328782509819455815301756717361332069811250996181881593041
690351598888519345807273866738589422879228499892086805825749279610484198444363463
244968487560233624827041978623209002160990235304369941849146314093431738143640546
253152096183690888707016768396424378140592714563549061303107208510383750510115747
704171898610687396965521267154688957035035402123407849819334321068170121005627880
235193033224745015853904730419957777093503660416997329725088687696640355570716226
844716256079882651787134195124665201030592123667719432527867539855894489697096409
754591856956380236370162112047742722836489613422516445078182442352948636372141740
238893441247963574370263755294448337998016125492278509257782562092622648326277933
386566481627725164019105900491644998289315056604725802778631864155195653244258698
294695930801915298721172556347546396447910145904090586298496791287406870504895858
671747985466775757320568128845920541334053922000113786300945560688166740016984205
580403363795376452030402432256613527836951177883863874439662532249850654995886234
281899707733276171783928034946501434558897071942586398772754710962953741521115136
835062752602326484728703920764310059584116612054529703023647254929666938115137322
753645098889031360205724817658511806303644281231496550704751025446501172721155519
486685080036853228183152196003735625279449515828418829478761085263981395599006737
648292244375287184624578036192981971399147564488262603903381441823262515097482798
777996437308997038886778227138360577297882412561190717663946507063304527954661855
096666185664709711344474016070462621568071748187784437143698821855967095910259686
200235371858874856965220005031173439207321139080329363447972735595527734907178379
342163701205005451326383544000186323991490705479778056697853358048966906295119432
473099587655236812859041383241160722602998330535370876138939639177957454016137223
618789365260538155841587186925538606164779834025435128439612946035291332594279490
433729908573158029095863138268329147711639633709240031689458636060645845925126994
655724839186564209752685082307544254599376917041977780085362730941710163434907696
423722294352366125572508814779223151974778060569672538017180776360346245927877846
585065605078084421152969752189087401966090665180351650179250461950136658543663271
254963990854914420001457476081930221206602433009641270489439039717719518069908699
860663658323227870937650226014929101151717763594460202324930028040186772391028809
786660565118326004368850881715723866984224220102495055188169480322100251542649463
981287367765892768816359831247788652014117411091360116499507662907794364600585194
199856016264790761532103872755712699251827568798930276176114616254935649590379804
583818232336861201624373656984670378585330527583333793990752166069238053369887956
513728559388349989470741618155012539706464817194670834819721448889879067650379590
366967249499254527903372963616265897603949857674139735944102374432970935547798262
961459144293645142861715858733974679189757121195618738578364475844842355558105002
561149239151889309946342841393608038309166281881150371528496705974162562823609216
807515017772538740256425347087908913729172282861151591568372524163077225440633787
593105982676094420326192428531701878177296023541306067213604600038966109364709514
141718577701418060644363681546444005331608778314317444081194942297559931401188868
331483280270655383300469329011574414756313999722170380461709289457909627166226074
071874997535921275608441473782330327033016823719364800217328573493594756433412994
302485023573221459784328264142168487872167336701061509424345698440187331281010794
512722373788612605816566805371439612788873252737389039289050686532413806279602593
038772769778379286840932536588073398845721874602100531148335132385004782716937621
800490479559795929059165547050577751430817511269898518840871856402603530558373783
242292418562564425502267215598027401261797192804713960068916382866527700975276706
977703643926022437284184088325184877047263844037953016690546593746161932384036389
313136432713768884102681121989127522305625675625470172508634976536728860596675274
086862740791285657699631378975303466061666980421826772456053066077389962421834085
988207186468262321508028828635974683965435885668550377313129658797581050121491620
765676995065971534476347032085321560367482860837865680307306265763346977429563464
371670939719306087696349532884683361303882943104080029687386911706666614680001512
114344225602387447432525076938707777519329994213727721125884360871583483562696166
198057252661220679754062106208064988291845439530152998209250300549825704339055357
016865312052649561485724925738620691740369521353373253166634546658859728665945113
644137033139367211856955395210845840724432383558606310680696492485123263269951460
359603729725319836842336390463213671011619282171115028280160448805880238203198149
309636959673583274202498824568494127386056649135252670604623445054922758115170931
492187959271800194096886698683703730220047531433818109270803001720593553052070070
607223399946399057131158709963577735902719628506114651483752620956534671329002599
439766311454590268589897911583709341937044115512192011716488056694593813118384376
562062784631049034629395002945834116482411496975832601180073169943739350696629571
241027323913874175492307186245454322203955273529524024590380574450289224688628533
654221381572213116328811205214648980518009202471939171055539011394331668151582884
368760696110250517100739276238555338617009e0
 
Here is a short program to calculate pi to an arbitrary degree of precision that I just wrote using the MPFR library.

I used the Gauss–Legendre iterative algorithm [this only requires 16 iterations to get a precise estimate at 150,515 decimal places utilizing 500,000 bit floating point variables].

{This is fairly memory intensive compared to a Machin-like algorithm, which uses arctans. Fortunately mpfr has an arctan function, which is the most difficult portion of this method to program, perhaps another day.}

Code:
// This is free code to calculate pi to an arbitrary degree of precision.
// There is no warranty or guarantee of any kind.
// The mpfr library has further restrictions.
// To Compile:
// gcc -o pi pi.c -lmpfr

#include <stdio.h>
#include <gmp.h>
#include <mpfr.h>
#include <stdlib.h>
#include <assert.h>


int pi(unsigned int *stop, unsigned int *bt){
    mpfr_t a, b, t, temp_A;
    mpfr_init2 (a, *bt);
    mpfr_init2 (b, *bt);
    mpfr_init2 (t, *bt);
    mpfr_init2 (temp_A, *bt);
    unsigned long p;

    //Initialize the values according to the Gauss–Legendre algorithm.
    mpfr_set_d (a, 1.0, MPFR_RNDD);

    mpfr_set_d (b, 2.0, MPFR_RNDD);
    mpfr_sqrt (b, b, MPFR_RNDD);
    mpfr_ui_div (b, 1, b, MPFR_RNDD);

    mpfr_set_d (t, 0.25, MPFR_RNDD);

    p = 1;

    //Perform iterations
    while (*stop>=1){
        (*stop)--;
    
        mpfr_add (temp_A, a, b, MPFR_RNDD);
        mpfr_div_ui (temp_A, temp_A, 2, MPFR_RNDD);
 
        mpfr_mul (b, a, b, MPFR_RNDU);
        mpfr_sqrt (b, b, MPFR_RNDD);
 
        mpfr_sub (a, a, temp_A, MPFR_RNDD);
        mpfr_sqr (a, a, MPFR_RNDD);
        mpfr_mul_ui (a, a, p, MPFR_RNDU);
        mpfr_sub (t, t, a, MPFR_RNDD);
 
        mpfr_set (a, temp_A, MPFR_RNDD);
 
        p = 2 * p;
    }

    //Perform final calculation
    mpfr_add (a, a, b, MPFR_RNDD);
    mpfr_sqr (a, a, MPFR_RNDD);

    mpfr_mul_ui (t, t, 4, MPFR_RNDU);

    mpfr_div (a, a, t, MPFR_RNDD);

    //Print Out Answer
    printf ("\n===================\nPI: ");
    mpfr_out_str (stdout, 10, 0, a, MPFR_RNDD);
    printf ("\n===================\n\n");

    mpfr_clear (a);
    mpfr_clear (b);
    mpfr_clear (t);
    mpfr_clear (temp_A);
    mpfr_free_cache ();
    return 0;
}


int main(int argc, char * argv[]){
  unsigned int i,b;

  if (argc <= 2){
    printf ("Usage: %s <number of iterations> <number of bits>\n", argv[0]);
    return 1;
  }

  i = atoi(argv[1]);
  b = atoi(argv[2]);

  assert( i >= 1);
  assert( b >= 1);
  pi(&i,&b);

  return 0;
}
 
Last edited:
Here's a fairly accurate estimate of PI, using the above program:


PI: 3.141592653589793238462643383279502884197169399375105820974
9445923078164062862089986280348253421170679821480865132823066
4709384460955058223172535940812848111745028410270193852110555
9644622948954930381964428810975665933446128475648233786783165
2712019091456485669234603486104543266482133936072602491412737
2458700660631558817488152092096282925409171536436789259036001
1330530548820466521384146951941511609433057270365759591953092
1861173819326117931051185480744623799627495673518857527248912
2793818301194912983367336244065664308602139494639522473719070
2179860943702770539217176293176752384674818467669405132000568
1271452635608277857713427577896091736371787214684409012249534
3014654958537105079227968925892354201995611212902196086403441
8159813629774771309960518707211349999998372978049951059731732
8160963185950244594553469083026425223082533446850352619311881
7101000313783875288658753320838142061717766914730359825349042
8755468731159562863882353787593751957781857780532171226806613
0019278766111959092164201989380952572010654858632788659361533
8182796823030195203530185296899577362259941389124972177528347
9131515574857242454150695950829533116861727855889075098381754
6374649393192550604009277016711390098488240128583616035637076
6010471018194295559619894676783744944825537977472684710404753
4646208046684259069491293313677028989152104752162056966024058
0381501935112533824300355876402474964732639141992726042699227
9678235478163600934172164121992458631503028618297455570674983
8505494588586926995690927210797509302955321231e0
 
If you modify that code you can use awk to determine the number of decimal places with this command:

awk '{ print length($0); }' pi_text

So 500,000 bits produces 150,519 characters or 150,515 decimal places [subtract the 3, the period, and the trailing e0, from the result. Does mpfr_out_str produce a newline character? I'm pretty sure it doesn't, that awk command doesn't reveal any blank lines.].

------
btw using 5 million bit floating point variables I yield a 1.5 mb text file with 1,505,150 decimal places [the awk command has a hard time with this size file].

It's possible to use the diff command with two different output files in order to discover if the algorithm has converged or not [both files should be the same].

In this case the algorithm converges somewhere around 32 iterations or less [21 iterations].
 
Last edited:
Here's a perl script which discovers how many iterations are required to converge to PI given an arbitrary number of bits [it assumes the only output from the pi program is the value of pi]. Perl appears to handle large strings better than awk.

Perl:
#!/usr/local/bin/perl -w
#---------------------------------#
# PROGRAM: pi_calc                #
#---------------------------------#

$na = $#ARGV + 1; #The number of arguments should equal 1
if($na != 1){
    print "\nSyntax Error: pi_calc <the number of bits>\n";
    exit(1);
}

$i = 1;
$OUT1 = `pi $i $ARGV[0]`;
while(1){
    $i++;
    $OUT2 = `pi $i $ARGV[0]`;
   
    if($OUT1 eq $OUT2){
        $i--;
        $length = length($OUT1) - 4;
        print "The Gauss–Legendre algorithm converges at $i iterations with a precision of $length decimal places using $ARGV[0] bit floating point variables.";
        exit(0);
    }
    $OUT1 = "$OUT2";
}
 
Here is a short program to calculate pi to an arbitrary degree of precision that I just wrote using the MPFR library.

I used the Gauss–Legendre iterative algorithm [this only requires 16 iterations to get a precise estimate at 150,515 decimal places utilizing 500,000 bit floating point variables].

{This is fairly memory intensive compared to a Machin-like algorithm, which uses arctans. Fortunately mpfr has an arctan function, which is the most difficult portion of this method to program, perhaps another day.}

Code:
// This is free code to calculate pi to an arbitrary degree of precision.
// There is no warranty or guarantee of any kind.
// The mpfr library has further restrictions.
// To Compile:
// gcc -o pi pi.c -lmpfr

#include <stdio.h>
#include <gmp.h>
#include <mpfr.h>
#include <stdlib.h>
#include <assert.h>


int pi(unsigned int *stop, unsigned int *bt){
    mpfr_t a, b, t, temp_A;
    mpfr_init2 (a, *bt);
    mpfr_init2 (b, *bt);
    mpfr_init2 (t, *bt);
    mpfr_init2 (temp_A, *bt);
    unsigned long p;

    //Initialize the values according to the Gauss–Legendre algorithm.
    mpfr_set_d (a, 1.0, MPFR_RNDD);

    mpfr_set_d (b, 2.0, MPFR_RNDD);
    mpfr_sqrt (b, b, MPFR_RNDD);
    mpfr_ui_div (b, 1, b, MPFR_RNDD);

    mpfr_set_d (t, 0.25, MPFR_RNDD);

    p = 1;

    //Perform iterations
    while (*stop>=1){
        (*stop)--;
   
        mpfr_add (temp_A, a, b, MPFR_RNDD);
        mpfr_div_ui (temp_A, temp_A, 2, MPFR_RNDD);
 
        mpfr_mul (b, a, b, MPFR_RNDU);
        mpfr_sqrt (b, b, MPFR_RNDD);
 
        mpfr_sub (a, a, temp_A, MPFR_RNDD);
        mpfr_sqr (a, a, MPFR_RNDD);
        mpfr_mul_ui (a, a, p, MPFR_RNDU);
        mpfr_sub (t, t, a, MPFR_RNDD);
 
        mpfr_set (a, temp_A, MPFR_RNDD);
 
        p = 2 * p;
    }

    //Perform final calculation
    mpfr_add (a, a, b, MPFR_RNDD);
    mpfr_sqr (a, a, MPFR_RNDD);

    mpfr_mul_ui (t, t, 4, MPFR_RNDU);

    mpfr_div (a, a, t, MPFR_RNDD);

    //Print Out Answer
    printf ("\n===================\nPI: ");
    mpfr_out_str (stdout, 10, 0, a, MPFR_RNDD);
    printf ("\n===================\n\n");

    mpfr_clear (a);
    mpfr_clear (b);
    mpfr_clear (t);
    mpfr_clear (temp_A);
    mpfr_free_cache ();
    return 0;
}


int main(int argc, char * argv[]){
  unsigned int i,b;

  if (argc <= 2){
    printf ("Usage: %s <number of iterations> <number of bits>\n", argv[0]);
    return 1;
  }

  i = atoi(argv[1]);
  b = atoi(argv[2]);

  assert( i >= 1);
  assert( b >= 1);
  pi(&i,&b);

  return 0;
}

Thanks a million BostonBSD Your example helped me a lot with understanding the MPFR library, and how to use specific functions. I have rewritten this procedure with use p variable and stop variable as buit-in mpfr library mpfr_t struct types. In this way the main loop in my version of this procedure may work for more steps.
 
I'm glad you liked it, page 13 of the manual also says
use mpfr_swap instead of mpfr_set whenever possible. This will avoid copying the signifi-
cands;

So I replaced:
Code:
mpfr_set (a, temp_A, MPFR_RNDD);

With
Code:
mpfr_swap (a, temp_A);

On line 48.

This works because the precision for both variables is the same and temp_A was already rounded.
 
Back
Top