Bill Allombert on Wed, 14 Oct 2009 15:52:41 +0200


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: using PARI library with arbitrary precision in C++


On Wed, Oct 14, 2009 at 12:12:41PM +0200, Michael Koehn wrote:
> Dear PARi users,
>
> please help me on my problem. I need to compute the complex gamma  
> function to arbitrary precision and wanted to use PARI's built-in  
> function. However I don't know how to convert the GEN value back to some 
> arbitrary precision type that I can use for my program in C++. (I am 
> using the arprec package.)

I see two ways:
1) you convert PARI t_REAL to strings (using GENtostr) and the strings to
arprec format. You can convert a string back to t_REALGEN using strtoGEN.

2) you can write a function that rebuild a arprec real from PARI t_REAL
sign, exponent and mantissa, and back.

> I am only able to convert my GEN with t_REAL back to double with  
> rtodbl(GEN x), but then I lose the high precision. Can anybody give me  
> please a hint on this?
>
> Also, I am wondering why I can only compute values for gamma(x) in  
> libpari to the 18th digit, if I print them with brute(). (See below.)

The issue is that dbltor return a real number with only 18 digits, so the
resulting gamma value will also have 18 digit. 

The precision passed to ggamma is only used for exact input (e.g. gamma(5))
and is in word+2 not in digits (so 50 si actually 463 digits in 32 bit and
926 in 64bits).

You can augment the precision of a number (by padding it with 0) using rtor,
e.g. rtor(dbltor(50.),8), but of course this assume the input is exactly equal
to the double.

> This is my example program:

> -----------
>
> #include </usr/local/include/pari/pari.h>
> #include <iostream>
> #include <iomanip>
> using namespace std;
>
> static GEN s;
>
> int main()	
> {
> 	pari_init(500000,0);
> 	
> //	output(ggamma(gadd(dbltor(200.),gmul(gi,dbltor(10.))),50));
> 	brute(ggamma(dbltor(50.),50),'f',100);
> 	cout << "\n" << setprecision(50) << rtodbl(ggamma(dbltor(50.),50)) <<  
> "\n";
> }
>
> ------------
>
> The returned values of gamma(50)
> -- in gp:(19:32) gp > gamma(50)
> %11 =  
> 608281864034267560872252163321295376887552831379210240000000000.0000000000000000000000000000000000000
>
> -- from libpari with brute(ggamma(dbltor(50.),50),'f',100):
> 6.082818640342675609E62

This result is correct.

> -- from libpari with brute(ggamma(dbltor(50.),50),'f',100) after  
> converting to double with rtodbl(GEN x):
> 6.0828186403426752248860160811673162316877754210242e+62

This result is wrong: you are creating artificial digits by expanding the double
outside its precision (53bits ~ 18digits).

Below is the correct result:

  6.0828186403426756087225216332129537688755283137921e62

Cheers,
Bill.