Re: Making a polynomial and evaluating it many times from C++

[Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>]

On Thu, Jul 22, 2021 at 06:07:47PM +0000, Neven Sajko wrote:
> Hello,
> 
> I want to make a certain polynomial in Pari (see script below) and
> evaluate it millions of times (plus some other bits) from my C or C++
> code.

Hello Neven,

Note that there are several dedicated algorithms for simultaneous
evaluation of a polynomial at multiple points that you could use,
depending of your choise of evaluation points, that would be
much faster than separate evaluations.

> I looked into gp2c, but it seems aimed at producing C code which will
> again be used from gp, instead of from other C or C++ code.

This is not the case, not sure why you get this
impression. If you need an example of a main() function, you
can look at examples/extgcd.c

> default(parisize, 2147483648);
> default(realbitprecision, 131072);
> default(seriesprecision, 2048);
> default(strictargs, 1);
> 
> default(format, "g.20");
> 
> truncated_power_series = simplify(sum(n = 0, 767, sqr(binomial(1/2, n)) * h^n));
> 
> ellipse_perimeter(a, b) = {
>         my(sum = a + b);
>         my(x = sqr((a - b) / sum));
>         Pi * sum * subst(truncated_power_series, h, x);
> };
> 
> ellipse_perimeter(0.5, 0.5)
> ellipse_perimeter(0.75, 0.25)
> ...
> 
> 
> The arguments to the ellipse_perimeter function must be IEEE 754
> binary64 numbers (known as 'double' in the C/C++ world), but they
> should then be converted to arbitrary precision floats, as the point
> of the computation is to get an accurate result.

You can call the libpari function 'dbltor' to convert C double to
PARI t_REAL.

> Another question, perhaps subjective: would it perhaps make more sense
> to just export the coefficients of the polynomial to my C++ code, and
> then use MPFR directly to evaluate it? That way Pari does only the
> symbolic manipulation, which it is more suited for, I guess?

PARI t_REAL format is not compatible with MPFR, so it is not that
simple.

For what is worth there are better formula for computing the ellipse
perimeter at high precision, see
http://www.ams.org/notices/201208/rtx120801094p.pdf
An Eloquent Formula for the Perimeter of an Ellipse
Semjon Adlaj
Notices of the AMS

and the simple-minded GP implementation below (function ellperimeter).

Cheers,
Bill

/*
http://www.ams.org/notices/201208/rtx120801094p.pdf
An Eloquent Formula for the Perimeter of an Ellipse 
Semjon Adlaj
Notices of the AMS
*/

magm(a,b)=
{
  my(c=0,eps=10^-(default(realprecision)-5));
  while(abs(a-b)>eps,
    my(u = sqrt((a-c)*(b-c))); 
    [a,b,c] = [(a+b)/2,c+u,c-u]);
  (a+b)/2
}

numE(g)=my(g2=g^2);intnum(x=0,[1,-1/2],my(x2=x^2);sqrt((1-g2*x2)/(1-x2)))

numK(g)=my(g2=g^2);intnum(x=0,[1,-1/2],my(x2=x^2);1/sqrt((1-g2*x2)*(1-x2)))

agmE(g)=my(b2=1-g^2);Pi/2/agm(1,sqrt(b2))*magm(1,b2)

agmK(g)=my(b2=1-g^2);Pi/2/agm(1,sqrt(b2))

agmPi(g)=my(b=sqrt(1-g^2));2*agm(1,b)*agm(1,g)/(magm(1,b^2)+magm(1,g^2)-1)

ellperimeter(a, b)=4*a*agmE(1-(b/a)^2)