Opened 3 years ago
Closed 3 years ago
#19822 closed enhancement (fixed)
Fast polynomial evaluation fmpz_poly/ZZX with mpfr/mpfi input
Reported by:  vdelecroix  Owned by:  

Priority:  major  Milestone:  sage7.1 
Component:  algebra  Keywords:  
Cc:  Merged in:  
Authors:  Vincent Delecroix  Reviewers:  Jeroen Demeyer 
Report Upstream:  N/A  Work issues:  
Branch:  232a575 (Commits)  Commit:  232a575c4933fc2e02d61fda511dca65b7a101e2 
Dependencies:  Stopgaps: 
Description (last modified by )
We implement functions that allow fast evaluations of integer polynomials with real input:
fmpz_poly_evaluation_mpfr(mpfr_t res, const fmpz_poly_t poly, const mpfr_t a)
fmpz_poly_evaluation_mpfi(mpfi_t res, const fmpz_poly_t poly, const mpfi_t a)
ZZX_evaluation_mpfr(mpfr_t res, ZZX_c poly, const mpfr_t a)
ZZX_evaluation_mpfi(mpfi_t res, ZZX_c poly, const mpfi_t a)
These functions are integrated in the polynomial code within the methods _eval_mpfr_
and _eval_mpfi_
.
See also the following discussion for upstream integrations:
And an algorithm for accuracte computations:
Attachments (10)
Change History (61)
comment:1 Changed 3 years ago by
 Branch set to u/vdelecroix/19822
 Commit set to 683cd0de16503dc913750fd7c417129409b010c6
 Status changed from new to needs_review
comment:2 Changed 3 years ago by
 Description modified (diff)
comment:3 Changed 3 years ago by
 Status changed from needs_review to needs_work
Some quick comments:
only do
>only does
 What's the point of the
mpfr_rnd_t rnd
parameter, given that you don't make any guarantees about exact rounding or rounding direction? Given the code, the only sensible rounding mode isMPFR_RNDN
.
 Did you investigate specialcasing small degrees (say, degree <= 1)?
 Both
fmpz_poly_degree(poly)
andZZX_deg
returnlong
, so the loop variable must belong
too.
comment:4 Changed 3 years ago by
 Please move the changes to
src/sage/rings/number_field/number_field_element.pyx
to a new ticket such that we can focus on the actual issue from the ticket here.
comment:5 Changed 3 years ago by
 Commit changed from 683cd0de16503dc913750fd7c417129409b010c6 to 131cd594f7a5b4977040b9f681dc8b02197e4247
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
131cd59  Trac 19822: review 1

comment:6 followup: ↓ 10 Changed 3 years ago by
 Description modified (diff)
 Status changed from needs_work to needs_review
For 2. I just forward the rounding mode to the mpfr functions. Having exact rounding for polynomials might actually be tricky... and also useful. The main reason for having this option is to handle simply polynomial evaluations with different RealField input.
For 3. I did not consider it on purpose. Having special case for degree 0,1 will slow down the code of higher degree cases which are indeed the interesting ones.
comment:7 Changed 3 years ago by
To be more complient with flint calling convention I will modify the names of the function eval > evaluation
.
comment:8 Changed 3 years ago by
 Commit changed from 131cd594f7a5b4977040b9f681dc8b02197e4247 to ad0721694614fe3546cdace8b3ec27e415539b74
Branch pushed to git repo; I updated commit sha1. New commits:
ad07216  Trac 19822: change method names + warning in doc

comment:9 Changed 3 years ago by
 Description modified (diff)
comment:10 in reply to: ↑ 6 Changed 3 years ago by
Replying to vdelecroix:
For 2. I just forward the rounding mode to the mpfr functions.
I know that but it doesn't make sense to do that!
What would be the point of doing every single operation with MPFR_RNDU
? That doesn't guarantee anything about the rounding of the result (if I take a number obtained with MPFR_RNDU
and multiply it with 1, it effectively becomes rounded with MPFR_RNDD
). Since polynomials use a mixture of additions and multiplications, it becomes difficult to avoid this problem.
comment:11 Changed 3 years ago by
 Status changed from needs_review to needs_work
Tiny detail: please keep the alphabetic order in src/module_list.py
.
comment:12 Changed 3 years ago by
Typo: Hornder
> Horner
.
comment:13 Changed 3 years ago by
Avoid sage/libs/ntl/decl.pxi
: just cimport
what you need from the .pxd
files in NTL.
comment:14 followup: ↓ 15 Changed 3 years ago by
Replying to jdemeyer:
Replying to vdelecroix:
For 2. I just forward the rounding mode to the mpfr functions.
I know that but it doesn't make sense to do that!
I see and you are definitely right. Then what should with p(x)
when p
is an integer polynomial and x
belongs to RealField(prec, rnd=MPFR_RNDU)
?
comment:15 in reply to: ↑ 14 Changed 3 years ago by
Just use MPFR_RNDN
for polynomial evaluation in all cases, it's the only sensible thing if you cannot guarantee anything about the rounding.
comment:16 Changed 3 years ago by
 Commit changed from ad0721694614fe3546cdace8b3ec27e415539b74 to 2453496c0c37e647e8a052c576dc975ff91e4ac0
comment:17 Changed 3 years ago by
 Description modified (diff)
 Status changed from needs_work to needs_review
comment:18 Changed 3 years ago by
 Description modified (diff)
comment:19 Changed 3 years ago by
remark: in flint they have two polynomial evaluation algorithms: Horner
and divide_and_conquer
and they hardcoded the threshold to 50.
comment:20 Changed 3 years ago by
 Commit changed from 2453496c0c37e647e8a052c576dc975ff91e4ac0 to 82f121bfac062934185f741610c7f97ddd504882
Branch pushed to git repo; I updated commit sha1. New commits:
82f121b  Trac 19822: add doctest for corner cases

comment:21 Changed 3 years ago by
 Description modified (diff)
comment:22 Changed 3 years ago by
 Commit changed from 82f121bfac062934185f741610c7f97ddd504882 to 5562a78b592125e698021743d545f176cc0604bb
Branch pushed to git repo; I updated commit sha1. New commits:
5562a78  Trac 19822: remove two Cython directives that belong to #19841

comment:23 Changed 3 years ago by
ready for review again.
comment:24 Changed 3 years ago by
ping?
comment:25 Changed 3 years ago by
The code looks good but I need to benchmark it.
comment:26 Changed 3 years ago by
You can also command it if you have precise requests ;)
BTW, I added special cases for degree 0 and 1 and I have seen no difference.
comment:27 followup: ↓ 28 Changed 3 years ago by
Benchmarks indicate that the bound degree < 10 makes no sense at all. The ratio in speed between the old and the new code does not depend on the degree at all.
Something which is strange: it looks like the precision of the real number/interval matters: for small precisions, the new code is faster but for large precisions (> 50000 bits) the old code is faster. This makes no sense to me! I would guess that both implementations should scale the same way if the precision is increased. Any clues?
comment:28 in reply to: ↑ 27 ; followup: ↓ 29 Changed 3 years ago by
Replying to jdemeyer:
Benchmarks indicate that the bound degree < 10 makes no sense at all. The ratio in speed between the old and the new code does not depend on the degree at all.
Something which is strange: it looks like the precision of the real number/interval matters: for small precisions, the new code is faster but for large precisions (> 50000 bits) the old code is faster. This makes no sense to me! I would guess that both implementations should scale the same way if the precision is increased. Any clues?
What kind of polynomials did you use? The old code tries to minimize the number of operation in evaluation. For example (x+1)^3
is infinitely smarter than 1 + x(3 + x(3 + x))
.
comment:29 in reply to: ↑ 28 Changed 3 years ago by
Replying to vdelecroix:
What kind of polynomials did you use?
R.<x> = ZZ[] pol = R.random_element(deg)
The old code tries to minimize the number of operation in evaluation.
I see. I guess that those "random" polynomials are somewhat sparse, so less operations are needed. I will try again with very dense polynomials.
comment:30 Changed 3 years ago by
 Reviewers set to Jeroen Demeyer
 Status changed from needs_review to needs_work
With dense polynomials and large precision, the timings of the old and new code are essentially the same. That's expected since the MPFR/MPFI operations dominate the runtime.
Can you remove the "degree < 10" check?
comment:31 followup: ↓ 32 Changed 3 years ago by
On second thought, what if the polynomial is x^1000
? Then we will lose a lot with your code.
comment:32 in reply to: ↑ 31 Changed 3 years ago by
Replying to jdemeyer:
On second thought, what if the polynomial is
x^1000
? Then we will lose a lot with your code.
Indeed. It is possible to implement a less stupid Horner with grouping. In other words x^6 + 3x^2 + 1
would be 1 + x^2(3 + x^4)
. There would be a small overhead for dense polynomials.
comment:33 Changed 3 years ago by
Or just skip the addition if the term is 0. It would be nice to benchmark that.
comment:34 Changed 3 years ago by
And there is no mpfi_pow
nor mpfi_pow_ui
!!
comment:35 Changed 3 years ago by
mpfr_pow
is *more* reliable than trial multiplication! Hence we should use it
sage: a = RR(1.1) sage: a*a*a*a == a**4 # understandable False sage: b = RealField(256)(1.1) sage: ans=b*b*b*b sage: (ansa**4) 4.44089209850063e16 sage: (ansa*a*a*a) 6.66133814775094e16
(EDIT: I was completely wrong in my first version)
comment:36 followup: ↓ 39 Changed 3 years ago by
I attached files which compare the various possibilities for Horner scheme:
eval1
: the current version which is straightforward Horner expansion (red)eval2
: the same but skip addition if the term is 0 (green)eval3
: usempfr_pow
(blue)
There are files degree_{deg}_{prec}.png
for fixed degrees where the number of terms vary and poly_{deg}_{nb_terms}.png
where the percision vary.
eval2
seems superior over eval1
. But eval3
seems only interesting for very sparse polynomials or very large precision. Indeed, for half dense polynomials (degree=20, nb_terms=13), eval3
gets better than the other two around prec=2^{8}
comment:37 Changed 3 years ago by
 Description modified (diff)
comment:38 Changed 3 years ago by
 Description modified (diff)
comment:39 in reply to: ↑ 36 Changed 3 years ago by
Do you have code for eval2
? I think that might be the best in the end, it is sort of a compromise.
comment:40 followup: ↓ 43 Changed 3 years ago by
 Commit changed from 5562a78b592125e698021743d545f176cc0604bb to 6970134df31c259484fd3d5fda0a2e6516a7db70
comment:41 Changed 3 years ago by
 Commit changed from 6970134df31c259484fd3d5fda0a2e6516a7db70 to 568f0d4c39d64f00a89b1a4293346b93317fcd9b
Branch pushed to git repo; I updated commit sha1. New commits:
568f0d4  Trac 19822: remove condition on degree for evaluation

comment:42 Changed 3 years ago by
 Status changed from needs_work to needs_review
comment:43 in reply to: ↑ 40 Changed 3 years ago by
Replying to git:
6970134 Trac 19822: ignore zero coefficients + copyright
Could you at least use the standard copyright template instead?
comment:44 Changed 3 years ago by
In doctests, use the x^n
notation instead of x**n
.
comment:45 Changed 3 years ago by
 Milestone changed from sage7.0 to sage7.1
 Status changed from needs_review to needs_work
comment:46 Changed 3 years ago by
Remove the mention of "complex" values:
This file provides fast evaluation of integer polynomials with a real or complex value
comment:47 Changed 3 years ago by
The doctests with the various rounding modes are obsolete, since the rounding mode isn't used. I suggest to remove those tests.
comment:48 Changed 3 years ago by
 Commit changed from 568f0d4c39d64f00a89b1a4293346b93317fcd9b to 232a575c4933fc2e02d61fda511dca65b7a101e2
comment:49 Changed 3 years ago by
 Status changed from needs_work to needs_review
comment:50 Changed 3 years ago by
 Status changed from needs_review to positive_review
comment:51 Changed 3 years ago by
 Branch changed from u/vdelecroix/19822 to 232a575c4933fc2e02d61fda511dca65b7a101e2
 Resolution set to fixed
 Status changed from positive_review to closed
New commits:
Trac 19822: polynomial evaluation
Trac 19822: call from python
Trac 19822: faster number field comparisons