Repository: sympy/sympy
Pull request: Polys: Refine dup_clear_denoms function and update tests
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->

Fixes: #26176

#### Brief description of what is fixed or changed

Refined the `dup_clear_denoms` function within the polys module, focusing on enhancing its efficiency and reliability when dealing with exact domains. The adjustments ensure that the function more accurately clears denominators, especially in scenarios involving symbolic expressions and rational numbers.

#### Other comments

Updated tests reflect these changes by focusing on:
- Ensuring the convert flag functions as intended across different scenarios.

#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->

* polys
  * Improved `dup_clear_denoms` function's efficiency in symbolic computations.

<!-- END RELEASE NOTES -->


Problem statement:
CoercionFailed when integrating exponential functions. But only for some values - why?
I am using sympy to solve and evaluate integrals. For some values, I get the error about coercion failed:

_raise CoercionFailed("Cannot convert %s of type %s from %s to %s" % (element, type(element), base, self))
sympy.polys.polyerrors.CoercionFailed: Cannot convert (8.48717*_Dummy_6413 + 2.999)/(8.0089*_Dummy_6413 + 2.83) of type <class 'sympy.polys.fields.FracElement'> from RR(_Dummy_6413) to RR[_Dummy_6413]_

Even debugging and going through other issues here in github couldn't help me to find a solution. The integral is bounded, and only using slightly different values works.
This is a minimal example:


`import sympy as sym
maximum_age = 1
a = sym.symbols('a')

# Example 1: these values work
[constant_1, constant_2] = [2.83, 2.99]

# Example 2: these values work
[constant_1, constant_2] = [1, 2.999]

# Example 3: these values don't work
[constant_1, constant_2] = [2.83, 2.999]

# This is the mathematical operation to do
factor_1 = 1 / (constant_1 * a + 1)
factor_2 = constant_2 * sym.exp(-a)
integral_term = sym.integrate(factor_1 * factor_2, (a, 0, maximum_age))`

Hints from issue discussion:
The problem is:
```python
In [11]: integrate(2.999*exp(-x)/(2.83*x + 1), (x, 0, 1))
---------------------------------------------------------------------------
CoercionFailed
```
More directly the problem is:
```python
In [18]: Poly([F(8.48717/(8.0089*x + 2.83)), F(0.0)], y, domain=F).clear_denoms()
Out[18]: (1.0⋅x + 0.353356890459364, Poly((8.48717*x + 2.999)/(8.0089*x + 2.83)*y, y, domain='RR(x)'))

In [19]: Poly([F(8.48717/(8.0089*x + 2.83)), F(0.0)], y, domain=F).clear_denoms(convert=True)
---------------------------------------------------------------------------
CoercionFailed
```
This is what happens if QQ is used instead of RR:
```python
In [20]: F2 = QQ.frac_field(x)

In [21]: Poly([F2(8.48717/(8.0089*x + 2.83)), F2(0.0)], y, domain=F2).clear_denoms()
Out[21]: 
⎛    5570104133810989                                                                               ↪
⎜x + ─────────────────, Poly(3900929712011235641763145417/3681104063018271642939543012*y, y, domain ↪
⎝    15763394698685100                                                                              ↪

↪          ⎞
↪ ='QQ(x)')⎟
↪          ⎠

In [22]: Poly([F2(8.48717/(8.0089*x + 2.83)), F2(0.0)], y, domain=F2).clear_denoms(convert=True)
Out[22]: 
⎛    5570104133810989                                                                               ↪
⎜x + ─────────────────, Poly(3900929712011235641763145417/3681104063018271642939543012*y, y, domain ↪
⎝    15763394698685100                                                                              ↪

↪          ⎞
↪ ='QQ[x]')⎟
↪          ⎠
```
So the problem seems to be in `clear_denoms`. When the domain is `RR(x)` is does not extract the polynomial denominator properly.

The problem seems to be this code:
https://github.com/sympy/sympy/blob/0fef9fb5da57ab3ae58a7cbd4f63981f243e2747/sympy/polys/densetools.py#L1199-L1208
It assumes that multiplying by the lcm of the denominator will give a rational function without polynomia denominator but for an inexact domain like `RR` it is possible that the denominators do not cancel completely due to rounding error. Actually the error looks much bigger than a typical rounding error:
```
1201 	        common = K1.lcm(common, K0.denom(c))
1202 	
1203 	    if not K1.is_one(common):
1204 	        f = dup_mul_ground(f, common, K0)
1205 	
1206 ->	    if not convert:
1207 	        return common, f
1208 	    else:
1209 	        return common, dup_convert(f, K0, K1)
1210 	
1211 	
(Pdb) p f
[(8.48717*x + 2.999)/(8.0089*x + 2.83), 0.0]
(Pdb) p common
x + 0.353356890459364
```
It would be worth trying to understand how the error is so large but in any case `dup_clear_denoms` should be written differently to handle inexact domains. Rather than just multiplying by the LCM of the denominators and hoping that it cancels it should separate the numerator explicitly and multiply each numerator by the lcm divided by the denominator.

Actually that is possibly a more efficient approach for exact fields including `QQ` and `QQ(x)` etc as well.

This seems to be the fix:
```diff
diff --git a/sympy/polys/densetools.py b/sympy/polys/densetools.py
index b56a4a773b..ecd1c98c3b 100644
--- a/sympy/polys/densetools.py
+++ b/sympy/polys/densetools.py
@@ -1199,13 +1199,20 @@ def dup_clear_denoms(f, K0, K1=None, convert=False):
     for c in f:
         common = K1.lcm(common, K0.denom(c))
 
-    if not K1.is_one(common):
-        f = dup_mul_ground(f, common, K0)
+    if K1.is_one(common):
+        if not convert:
+            return common, f
+        else:
+            return common, dup_convert(f, K0, K1)
+
+    # Use quo rather than exquo to handle inexact domains by discarding the
+    # remainder.
+    f = [K0.numer(c)*K1.quo(common, K0.denom(c)) for c in f]
 
     if not convert:
-        return common, f
+        return common, dup_convert(f, K1, K0)
     else:
-        return common, dup_convert(f, K0, K1)
+        return common, f
 
 
 def _rec_clear_denoms(g, v, K0, K1):
```
> It would be worth trying to understand how the error is so large

Actually the error is not as big as I thought:
```python
In [5]: monic((8.48717*x + 2.999))
Out[5]: 1.0⋅x + 0.353356890459364

In [6]: monic(8.0089*x + 2.83)
Out[6]: 1.0⋅x + 0.353356890459364

In [7]: cancel((8.48717*x + 2.999)/(8.0089*x + 2.83))
Out[7]: 
8.48717⋅x + 2.999
─────────────────
 8.0089⋅x + 2.83 
```
Interestingly with the fix above the integral evaluates but only for floats and not for rationals:
```python
In [9]: integrate(3*exp(-x)/(2*x + 1), (x, 0, 1))
Out[9]: 
  1               
  ⌠               
  ⎮      1        
3⋅⎮ ─────────── dx
  ⎮      x    x   
  ⎮ 2⋅x⋅ℯ  + ℯ    
  ⌡               
  0               

In [10]: integrate(3.0*exp(-x)/(2.0*x + 1), (x, 0, 1))
Out[10]: 
                     ⎛     ⅈ⋅π⎞                      ⎛     ⅈ⋅π⎞
- 2.47308190605019⋅Ei⎝0.5⋅ℯ   ⎠ + 2.47308190605019⋅Ei⎝1.5⋅ℯ   ⎠
```
In any case the issue can be closed with the diff above and a test for `clear_denoms`.
@oscarbenjamin Can you please tell me the reason for interchanging the return statments in the `dup_clear_denoms` ? 
There are two different domains `K0` and `K1`. The changed way of computing `f` gives a result over `K1 = RR[x]` rather than `K0 = RR(x)`.

See:
https://docs.sympy.org/latest/modules/polys/domainsintro.html

Commit message: This commit refines the `dup_clear_denoms` function within the polys module, focusing on enhancing its efficiency and reliability when dealing with exact domains. The adjustments ensure that the function more accurately clears denominators, especially in scenarios involving symbolic expressions and rational numbers.

Changes made include:
- Fine-tuning the logic for handling denominators to ensure robust performance across a variety of use cases.
- Streamlining the conversion process to guarantee that expressions are correctly converted to the target domain when requested.

The primary aim of these modifications was to address feedback and improve the function's handling of symbolic and rational expressions, thereby increasing the utility and reliability of polynomial operations within Sympy.

Updated tests reflect these changes by focusing on:
- Verifying the accurate clearing of denominators in exact domains.
- Ensuring the `convert` flag functions as intended across different scenarios.

While initial considerations included enhancements for inexact domain handling, this version concentrates on exact domains to solidify the foundation of `dup_clear_denoms`. Future work may revisit inexact domain improvements.

Fixes issue: #26176 (https://github.com/sympy/sympy/issues/26176)

improved code quality

Updated PR with comments

Updated code quality

Updated test cases for domain handling

Updated code quality

Added test with simplified expression

Updated test

Improved code quality

Improved code quality

improved code quality

Updated PR with comments

Updated code quality

Updated test cases for domain handling

Updated code quality

Added test with simplified expression

Updated test

Improved code quality

Improved code quality

Squashed commits

Updated PR with comments

Updated code quality

Updated test cases for domain handling

Updated code quality

Added test with simplified expression

Updated test

Improved code quality

Improved code quality

Patch to review:
diff --git a/sympy/polys/densetools.py b/sympy/polys/densetools.py
index b56a4a773b6b..ecd1c98c3be3 100644
--- a/sympy/polys/densetools.py
+++ b/sympy/polys/densetools.py
@@ -1199,13 +1199,20 @@ def dup_clear_denoms(f, K0, K1=None, convert=False):
     for c in f:
         common = K1.lcm(common, K0.denom(c))
 
-    if not K1.is_one(common):
-        f = dup_mul_ground(f, common, K0)
+    if K1.is_one(common):
+        if not convert:
+            return common, f
+        else:
+            return common, dup_convert(f, K0, K1)
+
+    # Use quo rather than exquo to handle inexact domains by discarding the
+    # remainder.
+    f = [K0.numer(c)*K1.quo(common, K0.denom(c)) for c in f]
 
     if not convert:
-        return common, f
+        return common, dup_convert(f, K1, K0)
     else:
-        return common, dup_convert(f, K0, K1)
+        return common, f
 
 
 def _rec_clear_denoms(g, v, K0, K1):
diff --git a/sympy/polys/tests/test_densetools.py b/sympy/polys/tests/test_densetools.py
index e8efb52f5230..4bd2f8029e1c 100644
--- a/sympy/polys/tests/test_densetools.py
+++ b/sympy/polys/tests/test_densetools.py
@@ -1,5 +1,6 @@
 """Tests for dense recursive polynomials' tools. """
 
+from sympy.core.symbol import symbols
 from sympy.polys.densebasic import (
     dup_normal, dmp_normal,
     dup_from_raw_dict,
@@ -28,7 +29,6 @@
     dup_sign_variations,
     dup_revert, dmp_revert,
 )
-
 from sympy.polys.polyclasses import ANP
 
 from sympy.polys.polyerrors import (
@@ -40,7 +40,7 @@
 
 from sympy.polys.specialpolys import f_polys
 
-from sympy.polys.domains import FF, ZZ, QQ, EX
+from sympy.polys.domains import FF, ZZ, QQ, EX, RR
 from sympy.polys.rings import ring
 
 from sympy.core.numbers import I
@@ -48,7 +48,6 @@
 from sympy.functions.elementary.trigonometric import sin
 
 from sympy.abc import x
-
 from sympy.testing.pytest import raises
 
 f_0, f_1, f_2, f_3, f_4, f_5, f_6 = [ f.to_dense() for f in f_polys() ]
@@ -613,6 +612,9 @@ def test_dup_sign_variations():
 
 
 def test_dup_clear_denoms():
+
+    x = symbols('x')
+
     assert dup_clear_denoms([], QQ, ZZ) == (ZZ(1), [])
 
     assert dup_clear_denoms([QQ(1)], QQ, ZZ) == (ZZ(1), [QQ(1)])
@@ -637,6 +639,9 @@ def test_dup_clear_denoms():
     assert dup_clear_denoms([EX(7)], EX) == (EX(1), [EX(7)])
     assert dup_clear_denoms([EX(sin(x)/x), EX(0)], EX) == (EX(x), [EX(sin(x)), EX(0)])
 
+    F = RR.frac_field(x)
+    result = dup_clear_denoms([F(8.48717/(8.0089*x + 2.83)), F(0.0)], F)
+    assert str(result) == "(x + 0.353356890459364, [1.05971731448763, 0.0])"
 
 def test_dmp_clear_denoms():
     assert dmp_clear_denoms([[]], 1, QQ, ZZ) == (ZZ(1), [[]])


Code review comments:
diff_hunk: @@ -613,6 +612,9 @@ def test_dup_sign_variations():
 
 
 def test_dup_clear_denoms():
+
+    x = symbols('x')
+
path: sympy/polys/tests/test_densetools.py
text: @user1:
```suggestion
```