On Tue, 8 Feb 2022 at 20:07, Luke Hartman <luke.hart...@tmlt.io> wrote:
Hi Luke,
> When you use evalf, you can provide an argument n, which is "the number of
> digits of precision" that you want
>
> One straightforward interpretation of this is that Sympy's evalf interface
> always returns a number whose nth significant digit is within 1 of the true
> answer (in either direction).
Be careful here as there are multiple things that are called evalf.
You are using what I would consider to be an internal function which
has a different interface from the public .evalf() method. Also you're
asking for a single binary digit of precision which is extreme. I
don't think that the public `.evalf()` method would ever call the
internal `evalf` function with such a low precision request.
> The only examples I have seen this break are around zero. In particular, if
> you have the wrong sign, then tweaking the nth significant digit is not going
> to help
>
> Here is a simple example where sp.core.evalf.evalf gives a different sign
> with different precision (so at least one of them has the wrong sign)
> ```
> expr = sp.fibonacci(27) - (sp.GoldenRatio**27)/sp.sqrt(5)
> re, _, re_acc, _ = sp.core.evalf.evalf(expr, prec=1, options={"strict": True})
> sign, man, exp, bc = re
> print([1, -1][sign]*man*2**exp)
> print(float(expr))
> ```
> prints
> ```
> -128
> 1.0182366178305287e-06
> ```
If you use the public .evalf method then the same sign is returned in all cases:
In [11]: for i in range(20):
...: print(i, expr.evalf(i))
0 0.e-6
1 1.e-6
2 1.0e-6
3 1.02e-6
4 1.018e-6
5 1.0182e-6
6 1.01824e-6
7 1.018237e-6
8 1.0182366e-6
9 1.01823662e-6
10 1.018236618e-6
11 1.0182366178e-6
12 1.01823661783e-6
13 1.018236617831e-6
14 1.0182366178305e-6
15 1.01823661783053e-6
16 1.018236617830529e-6
17 1.0182366178305288e-6
18 1.01823661783052880e-6
19 1.018236617830528800e-6
> I also dug into the code a bit to see what I could find, when you use
> "strict=True", behind the scenes, Sympy basically just checks that
> sp.core.evalf.evalf's returned re_acc (only concerned with) is at least as
> high as the desired precision. (Where re_acc seems to be the number of
> accurate bits in the number that re is a rounded version of)
>
> So it isn't really clear to me what the interpretation of n/prec/re_acc is
> for these cases around zero
The concept of precision is not well defined if the number is smaller
than its error bound.
> In a perfect world, I would like to have a function which takes the
> input/output of evalf and returns an upper and lower bound for the true value
> of the expression.
This is what evalf tries to do but firstly it does not return the
results in that form and secondly it is heuristic in some places so it
usually gives a good bound but there will be obscure cases where it
might not.
> As a related question, I would like to know how rigorous Sympy tries to be
> with adhering to a specific precision guarantee (for example mpmath has [this
> page](https://mpmath.org/doc/current/technical.html#correctness-guarantees))
Note that the goal in mpmath is somewhat different because it tries to
give a guarantee for individual operations rather than for the overall
result of evaluating a complex expression. Under the covers SymPy uses
mpmath for pretty much all of this and evalf does aim to satisfy
accuracy for the final output. There are definitely cases where it
could be improved though and overall the implemented methods have a
number of heuristics that might not be formally justifiable.
I think that it would be better in general to work with upper and
lower bounds using interval arithmetic or ball arithmetic. The
original author of the evalf module is Fredrik Johansson who has since
moved on and created the Arb library which does exactly that:
https://arblib.org/
The Arb library is included in python_flint:
https://fredrikj.net/python-flint/arb.html
I would definitely like to make python_flint/Arb something that SymPy
could make use of. Ultimately pretty much all of evalf should be
reworked along the lines of Arb.
--
Oscar
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to sympy+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/CAHVvXxRB%2B%2B%2BaY4YUtbXt0nzBgwM7NDy8TytrADDhYWTAXtWPsA%40mail.gmail.com.
