Various issues related to numerical evaluation

[GoogleCodeExporter]

I've been toying with various ways to implement a smart evalf function
using significance arithmetic, and I've come to the conclusion that it's
just too big a mess to get right. Instead of using significance arithmetic
to generate reliable numerical values (e.g. for proving inequalities), we
should use interval arithmetic right away; it is rigorous, and simpler to
implement.

Unfortunately, interval computations are typically 2-10 times slower than
direct float computations. Complex intervals might be even worse. In
addition, there are the standard pitfalls of interval arithmetic: x*x !=
x**2, many numerical algorithms don't converge, etc.

For these reasons, I think there should be one evalf function which just
uses Float arithmetic. It could use a settable extra working precision,
e.g. 10 digits, and that will cover nearly all rounding errors, but it
should not be required to check for cancellations that cause the result to
be wrong in extremely badly conditioned cases. There should be a second,
independent function for generating a rigorous interval or complex interval
for a given expression.

Third, there should be way to compile an expression to an optimized Python
function that uses only floats, for matrix and plotting performance. We
already have that in sympy (lambdify) -- what we should look at in
sympycore is how to integrate it better, and perhaps if it should be part
of a more general functionality for code generation in multiple languages.
It would be extremely cool if lambdify could check if a C compiler is
available and then dynamically both generate and load a C extension that
implements the function (I don't know if this is feasible).

A final detail. I think approximate numbers should be contagious like in
Mathematica (and even pure Python). That is, every numerical part of an
expression that comes in part with a Float should be evaluated to a Float,
e.g. 1.0*x*Sqrt(2) -> 1.414*x. Would this be difficult to implement with
our current implementation of the main arithmetic classes, and would it be
expensive?

Original issue reported on code.google.com by fredrik....@gmail.com on 3 Dec 2007 at 6:19

[GoogleCodeExporter]

> and perhaps if it should be part of a more general functionality for code
generation in multiple languages.

This part can and probably should be done in sympy, not sympycore (to make 
Ondrej
happy ;-)

Original comment by fredrik....@gmail.com on 3 Dec 2007 at 6:59

[GoogleCodeExporter]

> A final detail. I think approximate numbers should be contagious like in

Agree.

> already have that in sympy (lambdify) -- what we should look at in
> sympycore is how to integrate it better, and perhaps if it should be part
> of a more general functionality for code generation in multiple languages.

Agree. The code generation should probably go to a module. Pearu, is code 
generation
something you would like to have in sympycore? I think it belongs to IO.

> be wrong in extremely badly conditioned cases. There should be a second,
> independent function for generating a rigorous interval or complex interval
> for a given expression.

Does Mathematica have something like that? I think it'd be cool to have 
something
like that. But it belongs to mpmath, right? Then we'll just call it from 
sympy/sympycore.

Original comment by ondrej.c...@gmail.com on 8 Dec 2007 at 8:46

[GoogleCodeExporter]

Code generation has not high priority at the moment but
this is certainly something that I would need in projects
that will use python based CA, e.g. in g3f2py.

Original comment by pearu.peterson on 8 Dec 2007 at 9:19

[GoogleCodeExporter]

> Does Mathematica have something like that? I think it'd be cool to have 
something
> like that. But it belongs to mpmath, right? Then we'll just call it from
> sympy/sympycore.

Mathematica supports interval arithmetic, but it only works for simple 
functions (as
of version 5, at least).

Interval arithmetic can be implemented easily in mpmath; it's usually just a 
matter
of evaluating the function at both endpoints, both rounded down and up (some 
more
care is needed for non-monotone functions). See
http://code.google.com/p/mpmath/wiki/IntervalArithmetic for a usage demo (the 
code
used to generate the examples is not in SVN yet... still working on that).

Complex intervals are more difficult, and to the extent they can be supported, 
they
probably require a combination of symbolic and numerical evaluation.

Original comment by fredrik....@gmail.com on 8 Dec 2007 at 9:39