Opened 3 years ago
Closed 3 years ago
#23140 closed enhancement (fixed)
GaussLegendre Integrator
Reported by:  sasha.zotine  Owned by:  

Priority:  major  Milestone:  sage8.0 
Component:  numerical  Keywords:  GaussLegendre Integration sd86.5 
Cc:  Merged in:  
Authors:  Nils Bruin, Alexandre Zotine  Reviewers:  Aly Deines 
Report Upstream:  N/A  Work issues:  
Branch:  a455f3e (Commits)  Commit:  a455f3e4532f4b642543aaa88aa701b1b3b8d490 
Dependencies:  Stopgaps: 
Description
This code implements a method for doing GaussLegendre integration on a fastcallable polynomial. Designed with multiprecision in mind.
Change History (15)
comment:1 Changed 3 years ago by
 Branch set to u/nbruin/gauss_legendre_integrator
comment:2 Changed 3 years ago by
 Commit set to fa4495a1ce973563a586538f773296035d9fbc83
 Keywords sd86.5 added
comment:3 Changed 3 years ago by
 Commit changed from fa4495a1ce973563a586538f773296035d9fbc83 to fec8f0026f21bbc32216c440cd89a418a3d81f33
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
fec8f00  Creation of a vectorized gauss_legendre integrator

comment:4 Changed 3 years ago by
 Commit changed from fec8f0026f21bbc32216c440cd89a418a3d81f33 to f1287234af088ddf7d631e75aa6a9581c09b860c
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
f128723  Creation of a vectorized gauss_legendre integrator

comment:5 Changed 3 years ago by
Ready for review. Code is basically a translation of what is in mpmath. See for instance
https://github.com/fredrikjohansson/mpmath/blob/master/mpmath/calculus/quadrature.py#L203
for comparison.
(most computations are straightforward. It's particularly the error estimator that is tricky, and that is basically copied literally; modulo some logbase changes)
comment:6 Changed 3 years ago by
 Commit changed from f1287234af088ddf7d631e75aa6a9581c09b860c to d68683278e55e155745e9330e62f0a1617eb7e0d
Branch pushed to git repo; I updated commit sha1. New commits:
d686832  Change to node for degree=1, based on bugfix in mpmath:

comment:7 Changed 3 years ago by
 Priority changed from minor to major
 Status changed from new to needs_review
For documentation completeness, here is an issue tracker for an issue concerning log bases in the error estimation routine: https://github.com/fredrikjohansson/mpmath/pull/329
comment:8 Changed 3 years ago by
Minor comments:
 In Cython, you can use the libc math library instead of the Python math library. It suffices to replace
import math
bycimport libc.math as math
.
 To ensure Python 3 compatibility, use
from __future__ import absolute_import, division, print_function
.
 Keep in mind PEP 8 which says that assignments should have spaces: replace
Rout=RealField(prec)
byRout = RealField(prec)
.
comment:9 followup: ↓ 12 Changed 3 years ago by
 You use
cdef int
in various places. Only do this when you know that the variable must always be small. If there any tiny chance that this will ever be larger than2^31
, usecdef long
or some other type instead. Even if you cannot support GaussLegendre integration in "degree">= 31
now, there is no reason to a priori assume that it will never be needed.
 In documentation, use
^
instead of**
.
(*) I find "degree" a confusing name. Is this standard in Gauss–Legendre integration theory? I would expect the degree to be one less than the number of nodes.
comment:10 Changed 3 years ago by
 Commit changed from d68683278e55e155745e9330e62f0a1617eb7e0d to 8c823208f6305337954d5c1b6f90f3ee655bcb6a
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
8c82320  Change to node for degree=1, based on bugfix in mpmath:

comment:11 Changed 3 years ago by
(Pushed an amended branch because git had filled in a bogus author. cocalc projects and git mix badly)
comment:12 in reply to: ↑ 9 Changed 3 years ago by
Replying to jdemeyer: Replying to jdemeyer:
Minor comments:
 In Cython, you can use the libc math library instead of the Python math library. It suffices to replace
import math
bycimport libc.math as math
.
I tried and it generated errors for me (it got confused with complex numbers in a way that python's math doesn't). It's only used in a noncritical part (computing the initial value for newton iteration), so it doesn't need to be fast. The python library seems to work fine.
 To ensure Python 3 compatibility, use
from __future__ import absolute_import, division, print_function
.
Thanks
 Keep in mind PEP 8 which says that assignments should have spaces: replace
Rout=RealField(prec)
byRout = RealField(prec)
.
Done
 You use
cdef int
in various places. Only do this when you know that the variable must always be small. If there any tiny chance that this will ever be larger than2^31
, usecdef long
or some other type instead. Even if you cannot support GaussLegendre integration in "degree">= 31
now, there is no reason to a priori assume that it will never be needed.
Good point.
 In documentation, use
^
instead of**
.
Good point
(*) I find "degree" a confusing name. Is this standard in Gauss–Legendre integration theory? I would expect the degree to be one less than the number of nodes.
Frederic uses it in mpmath. He scales the parameter so that an increase roughly doubles the precision. I agree, though: it's much easier to document if you just give the degree of the legendre polynomials. It would be nice to use this routine as a dropin replacement in sage's mpmath too.
comment:13 Changed 3 years ago by
 Commit changed from 8c823208f6305337954d5c1b6f90f3ee655bcb6a to a455f3e4532f4b642543aaa88aa701b1b3b8d490
Branch pushed to git repo; I updated commit sha1. New commits:
a455f3e  some PEP8, reviewer comments; rationalized "degree" parameter for "nodes"

comment:14 Changed 3 years ago by
 Reviewers set to Aly Deines
 Status changed from needs_review to positive_review
If you really really wanted to you could clean up a tiny bit to meet all PEP 8 standards, but nothing else needs work or to be changed.
comment:15 Changed 3 years ago by
 Branch changed from u/nbruin/gauss_legendre_integrator to a455f3e4532f4b642543aaa88aa701b1b3b8d490
 Resolution set to fixed
 Status changed from positive_review to closed
New commits:
initial checkin of gausslegendre integrator