precision loss in mpmath

[Tobias Hartung]

Hi everyone,

I am trying to integrate polynomials on an n-sphere and implemented (see attached) the algorithm published in Alan Genz "Fully Symmetric Interpolatory Rules for Multiple Integrals over Hyper-Spherical Surfaces" Journal of Computational and Applied Mathematics 157: 187-195, 2003. On double precision, everything is fine; however, I lose precision in some integrals using more than double prec.

I am using randomized samples on the sphere which (by construction) should integrate the polynomials (that I tested) exactly. However, integrating all polynomials up to degree 6 on a 5-dimensional sphere with points that should integrate all polynomials up to degree 7 exactly, I obtain the following test results :

mp.prec = 1024 (aka mp.dps = 307)

polynomial  |  degree of polynomial  |  log_10 ( difference of two randomized integrations )  |  log_10 ( relative error of weights )  |  analytic value of integral

[0, 0, 0, 0, 0, 0]      0      -inf                     -306.962262      31.0062767

[1, 0, 0, 0, 0, 0]      1      -308.111506      -306.962262      0.0

[2, 0, 0, 0, 0, 0]      2      -14.519354        -306.962262      5.16771278

[1, 1, 0, 0, 0, 0]      2      -14.8787358      -306.962262      0.0

[3, 0, 0, 0, 0, 0]      3      -308.583548      -306.962262      0.0

[2, 1, 0, 0, 0, 0]      3      -309.157806      -306.962262      0.0

[1, 1, 1, 0, 0, 0]      3      -309.759866      -306.962262      0.0

[4, 0, 0, 0, 0, 0]      4      -14.6442927      -306.962262      1.93789229

[3, 1, 0, 0, 0, 0]      4      -15.3047045      -306.962262      0.0

[2, 2, 0, 0, 0, 0]      4      -15.0673319      -306.962262      0.645964098

[2, 1, 1, 0, 0, 0]      4      -16.0594155      -306.962262      0.0

[1, 1, 1, 1, 0, 0]      4      -31.8146065      -306.962262      0.0

[5, 0, 0, 0, 0, 0]      5      -309.37189        -306.962262      0.0

[4, 1, 0, 0, 0, 0]      5      -309.481766      -306.962262      0.0

[3, 2, 0, 0, 0, 0]      5      -309.739262      -306.962262      0.0

[3, 1, 1, 0, 0, 0]      5      -310.095075      -306.962262      0.0

[2, 2, 1, 0, 0, 0]      5      -310.474598      -306.962262      0.0

[2, 1, 1, 1, 0, 0]      5      -310.714154      -306.962262      0.0

[1, 1, 1, 1, 1, 0]      5      -311.46643        -306.962262      0.0

[6, 0, 0, 0, 0, 0]      6      -14.7692314      -306.962262      0.968946146

[5, 1, 0, 0, 0, 0]      6      -15.6057345      -306.962262      0.0

[4, 2, 0, 0, 0, 0]      6      -15.4314091      -306.962262      0.193789229

[3, 3, 0, 0, 0, 0]      6      -15.8275833      -306.962262      0.0

[4, 1, 1, 0, 0, 0]      6      -16.5822942      -306.962262      0.0

[3, 2, 1, 0, 0, 0]      6      -16.4753651      -306.962262      0.0

[2, 2, 2, 0, 0, 0]      6      -16.0997045      -306.962262      0.0645964098

[3, 1, 1, 1, 0, 0]      6      -32.3374853      -306.962262      0.0

[2, 2, 1, 1, 0, 0]      6      -17.3552924      -306.962262      0.0

[2, 1, 1, 1, 1, 0]      6      -33.4100353      -306.962262      0.0

[1, 1, 1, 1, 1, 1]      6      -49.1751847      -306.962262      0.0

It should be noted that the polynomial are written in the following form: a list p of length n represents the polynomial  

z_0^{p[0]}  z_1^{p[1]}  z_2^{p[2]}  ...  z_{n-1}^{p[n-1]}.

If we look at the second to last column, we can see that the weights add are exact with respect to mp.prec=1024 but do not add up to 1.0 exactly. Hence, the loss in precision is unlikely to be caused by the weights being only on double prec (and even if it is, then I don't understand why because they are calculated with mp.prec=1024, as well). 

The really interesting column is the middle one. Here, I took two randomized set of integration points (both should integrate all polynomials exactly up to machine error) and printed the decadic logarithm of the absolute value of their difference. In other words, we are expecting the middle column to be populated with numbers in the vicinity of -307 (just like the fourth column). However, we only get these values for odd degree polynomials (and the constant 1 but that is just the sum of weights, i.e., the randomization doesn't do anything). For even degree polynomials, we have significantly reduced precision. 

Does anyone have an idea what might be the reason for such behavior?

Thank you very much,

Tobias

[Aaron Meurer]

You may want to crosspost this to the mpmath list.

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/cf76f223-96a0-4789-8530-db224622f6bb%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Tobias Hartung]

Thank you and sorry about the wrong group post. Do you happen to know, how I can crosspost it?

Tobias

[Aaron Meurer]

Here is the mpmath mailing list https://groups.google.com/forum/#!forum/mpmath.

Aaron Meurer

On Thu, Jan 14, 2016 at 11:17 AM, Tobias Hartung

> https://groups.google.com/d/msgid/sympy/015dd6ad-d2e8-495c-82f7-1c157b0c7ab0%40googlegroups.com.

[Tobias Hartung]

Yes, I found it, but it won't let me re-post.

Tobias

[Tobias Hartung]

Ok. Apparently the post had to be validated on the mpmath list first.