Closed jdemeyer closed 8 years ago
Branch: u/jdemeyer/ticket/19797
Isn't this sacrificing a lot of numerical performance? At least IEEE 754-2008 contains fused multiply-add.
The doctest is
sage: x = var('x', domain='real')
sage: s = abs((1+I*x)^4); s
(I*x + 1)^2*(-I*x + 1)^2
sage: s._plot_fast_callable(x)
<sage.ext.interpreters.wrapper_py.Wrapper_py object at ...>
sage: s._plot_fast_callable(x)(10)
10201
sage: abs((I*10+1)^4)
10201
sage: plot(s)
Graphics object consisting of 1 graphics primitive
Isn't the real problem that abs does not default to hold=True
sage: ((1+I*x)^4).abs()
(I*x + 1)^2*(-I*x + 1)^2
so something that ought to be manifestly real is not. Compare:
sage: abs((1+I*x)^4)._plot_fast_callable(x)(RDF(10))
10201.00000000001
sage: type(_)
<type 'sage.rings.complex_double.ComplexDoubleElement'>
vs.
sage: ((1+I*x)^4).abs(hold=True)._plot_fast_callable(x)(RDF(10))
10201.000000000005
sage: type(_)
<type 'sage.rings.real_double.RealDoubleElement'>
New commits:
e468673 | Simplify build of interpreters by skipping header files |
694feb2 | Build GSL in IEEE 754 compliant mode |
Replying to @vbraun:
Isn't this sacrificing a lot of numerical performance?
I have no idea but I doubt it. For example, Intel architectures added this instruction only very recently.
Yes but CPU's haven't seen any general speedups recently, either. Only more specialized instructions, and if you don't use them then you are not participating in progress.
Replying to @vbraun:
Yes but CPU's haven't seen any general speedups recently, either. Only more specialized instructions, and if you don't use them then you are not participating in progress.
I would say that at least in the specific case of complex multiplication, the use of FMA instructions is wrong.
Alternatively, we could probably play some tricks with volatile
to fix complex multiplication.
So the question is: Should we guarantee that x*x.conj()
has no numerical noise as imaginary part.
I'm not sure.
What I am pretty sure about is: multiplying x*x.conj()
as complex numbers is a bad way of evaluating the absolute value. Fast callables using abs should return real and not complex values. And for that we need to hold abs. At which point the issue of the numerical noise in the imaginary part is pretty meaningless.
Replying to @vbraun:
sage: x = var('x', domain='real') sage: s = abs((1+I*x)^4); s (I*x + 1)^2*(-I*x + 1)^2
Isn't the real problem that abs does not default to
hold=True
Maybe it's better to expand automatically if the result contains I
. Automatic expansion also resolved a similar problem in #18952
sage: x = var('x', domain='real')
sage: ((1+I*x)^4).abs().expand()
x^4 + 2*x^2 + 1
so something that ought to be manifestly real is not.
It is but only when expanded.
Maybe it's better to expand automatically if the result contains
I
.
Also only if the result is not of form abs(...)
And not of the form 1+abs(...)
and so on
Automatically expanding isn't good enough either in general:
sage: x = var('x', domain='real')
sage: ((exp(I*x)+exp(-I*x))^4).abs().expand()
e^(4*I*x) + 4*e^(2*I*x) + 4*e^(-2*I*x) + e^(-4*I*x) + 6
Upstream: Reported upstream. Developers acknowledge bug.
Description changed:
---
+++
@@ -20,3 +20,4 @@
Graphics object consisting of 0 graphics primitives
**********************************************************************
+See https://github.com/pynac/pynac/issues/117
Replying to @vbraun:
- Even besides FMA there is always the possibility of 80-bit x87 instructions messing things up
With x87, it is reasonable to assume that the whole imaginary part will be evaluated with 80 bits, so there won't be numerical noise either (see also the successful i386 buildbot reports).
Replying to @vbraun:
And not of the form
1+abs(...)
and so on
This can be fixed in GiNaC::abs_eval()
Automatically expanding isn't good enough either in general:
sage: x = var('x', domain='real') sage: ((exp(I*x)+exp(-I*x))^4).abs().expand() e^(4*I*x) + 4*e^(2*I*x) + 4*e^(-2*I*x) + e^(-4*I*x) + 6
Why not? This is real as well and could be rewritten as a sum of trig function expressions.
Replying to @rwst:
This is real as well and could be rewritten as a sum of trig function expressions.
Volker's point is about being "manifestly" real, which means: an expression which will always give a real result when evaluated, even if sub-expressions introduce numerical noise.
An expression like
e^(4*I*x) + e^(-4*I*x)
is not manifestly real, since there is no guarantee that the imaginary parts of the two terms will exactly cancel out.
Well then, exclude functions in the treewalk in the search for I
. Treewalks are very fast in comparison to other Pynac operations, and those are only a small fraction when compared to what Python costs.
Its not just that there is no guarantee against numerical noise in the imaginary part of exp(I*x)+exp(-I*x)
, the return type of the fast callable also changes from real to complex (comment:2). The symbolic expression containing abs() should only do automatic simplifications that are again manifestly real.
Expanding polynomials over QQ[I] is safe, and that is the only safe case I can think of
Transcendental functions are not safe
Expanding polynomials whose coefficients are radical expressions over QQ[I] are not safe
Thoughts?
Replying to @vbraun:
- Expanding polynomials whose coefficients are radical expressions over QQ[I] are not safe
Is this decidable without expansion? Simply weeding out rational exponents will make (I<sup>(5/4))*(I</sup>(3/4))
a false negative. OTOH first expanding then checking would be costly.
If so, I no longer think it useful.
So, you're proposing to always hold abs()
and probably have a keyword for x*conjugate(x)
output?
Even with expansion I think it can be difficult to tell if a radical is real or not; Cardano's formula which can't be written in real radicals even if the roots are real comes to mind.
Of course it is possible to decide if a radical is real or not, for example using QQbar.
Can't we just abs(hold=True)
by default? You can always simplify or set hold=False
if you don't want that.
Replying to @vbraun:
Can't we just
abs(hold=True)
by default?
We could certainly do that as a stopgap-style solution. However, this is certainly going to break some doctests...
Replying to @jdemeyer:
Replying to @vbraun:
Can't we just
abs(hold=True)
by default?We could certainly do that as a stopgap-style solution. However, this is certainly going to break some doctests...
I get only two fails in functions/
, symbolic/
, calculus/
, doc/
if I outcomment this snippet in Pynac which is responsible for the specific behaviour:
https://github.com/pynac/pynac/blob/master/ginac/inifcns.cpp#L253
So the result is arrived at by two code parts: first abs
of real power (or power of positive) is rewritten as power of abs
and, second, even power of abs(x)
is rewritten as x*x.conj()
to the half power.
So I'll just disable the abs(ex<sup>real)-->abs(ex)</sup>real
part in that Pynac ticket if there are no objections.
Replying to @rwst:
second, even power of
abs(x)
is rewritten asx*x.conj()
to the half power.
It's really this second part which is problematic.
Description changed:
---
+++
@@ -1,12 +1,13 @@
-The following can happen on some hardware architectures due to fused multiply-add instructions:
+Calling `abs()` on certain expressions can yield something which is not guaranteed to be evaluated as a real number:
-sage: x = CDF(0.99, 0.2) -sage: x x.conj() -1.0201 + 1.1102230246251575e-19I +sage: abs(x^2) +x*conjugate(x)
-It's annoying because it causes doctest failures, for example
+In the presence of numerical noise, the expression `x*conjugate(x)` can actually a non-zero imaginary part. And even if there is no noise, the type would be wrong: complex instead of real.
+
+This causes the following doctest failure:
sage -t --long --warn-long 62.3 src/sage/symbolic/expression.pyx @@ -20,4 +21,5 @@ Graphics object consisting of 0 graphics primitives
+
See https://github.com/pynac/pynac/issues/117
Changed dependencies from #19796 to none
Description changed:
---
+++
@@ -5,7 +5,7 @@
x*conjugate(x)
-In the presence of numerical noise, the expression x*conjugate(x)
can actually a non-zero imaginary part. And even if there is no noise, the type would be wrong: complex instead of real.
+In the presence of numerical noise, the expression x*conjugate(x)
can have a non-zero imaginary part. And even if there is no noise, the type would be wrong: complex instead of real.
This causes the following doctest failure:
The simplification abs(ex<sup>real)-->abs(ex)</sup>real
is safe since the result is still manifestly real.
Description changed:
---
+++
@@ -2,6 +2,8 @@
sage: abs(x^2) +xconjugate(x) +sage: abs(x)^2 xconjugate(x)
This was introduced with pynac-0.3.9.2 from GiNaC. I'll revert it for pynac-0.5.4.
Reviewer: Jeroen Demeyer
"Duplicate" of #19819.
Changed author from Jeroen Demeyer to none
Changed branch from u/jdemeyer/ticket/19797 to none
Calling
abs()
on certain expressions can yield something which is not guaranteed to be evaluated as a real number:In the presence of numerical noise, the expression
x*conjugate(x)
can have a non-zero imaginary part. And even if there is no noise, the type would be wrong: complex instead of real.This causes the following doctest failure:
See https://github.com/pynac/pynac/issues/117
Upstream: Reported upstream. Developers acknowledge bug.
CC: @rwst
Component: packages: standard
Reviewer: Jeroen Demeyer
Issue created by migration from https://trac.sagemath.org/ticket/19797