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

[Neven Sajko <nsajko@gmail.com>]

On Thu, 22 Jul 2021 at 23:27, Bill Allombert
<Bill.Allombert@math.u-bordeaux.fr> wrote:
>
> 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.

Do these come with Pari?

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

Yes, I figured out that gp2c is useful to me after all. Also, thank
you for the example of initializing Pari.

> > 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.

Thanks!

> > 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.

Noted.

> 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).

FTR, your function for the perimeter of an ellipse is not accurate.

? ellipse_perimeter(0.75, 0.25)
%9 = 3.3412233051388145575
? ellperimeter(0.75, 0.25)
%19 = 3.5553025205101225533

An independent implementation also agrees with my answer (3.34122...).

Thanks,
Neven