C/C++ Multiple Precision Variables in FreeBSD

BostonBSD

Member

Reaction score: 16
Messages: 49

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]


Code:
//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;
}
 
OP
B

BostonBSD

Member

Reaction score: 16
Messages: 49

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
 
OP
B

BostonBSD

Member

Reaction score: 16
Messages: 49

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:
OP
B

BostonBSD

Member

Reaction score: 16
Messages: 49

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
 
OP
B

BostonBSD

Member

Reaction score: 16
Messages: 49

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:
OP
B

BostonBSD

Member

Reaction score: 16
Messages: 49

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.

Code:
#!/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";
}
 
Top