Sympy vs Numpy, better accuracy in precision?

[Amy Valhausen]

Sympy vs Numpy, better accuracy in precision?

I've been trying to solve a problem with numpy and other code routines
to raise a base to a large power and then take the modulus.

Precision accuracy is very important, speed isnt as much - although it would
be convenient if I didnt have to wait a long long time for processing.

Constraints Im under is that Im working on a winxp system, Im using
python 3.4 and numpy version 1.1.

When using numpy Ive been using the code lines;

import numpy
(np.longdouble(1.4142)** 6000 )%400

I am not sure how accurate the result is, and I have tried using other methods too
but recently I found a post comparing sympy to numpy and it seems someone
is claiming that sympy can return superior precision results, can anyone confirm
that this is true and do you know if this would be a good solution to run on my system?

Link with info comparing the below is shown ;

http://stackoverflow.com/questions/25181642/how-set-numpy-floating-point-accuracy

In normal numpy use, the numbers are double. Which means that the accuracy will be less than 16 digits. Here is a solved subject that contains the same problematic ...

If you need to increase the accuracy, you can use symbolic computation .... The library mpmath ... is a quiet good one. The advantage is that you can use limitless precision. However, calculations are slower than what numpy can do.

Here is an example:

# The library mpmath is a good solution
>>> import sympy as smp
>>> mp = smp.mpmath

>>> mp.mp.dps = 50  # Computation precision is 50 digits
# Notice that the difference between x and y is in the digit before last (47th)
>>> x = smp.mpmath.mpf("0.910221324013388510820732335560023784637451171875")
>>> y = smp.mpmath.mpf("0.910221324013388510820732335560023784637451171865")
>>> x - y  # Must be equal to 1e-47 as the difference is on the 47th digit
mpf('1.000014916280995001003481719184726944958705912691304e-47')

You can't do better with numpy. You can calculate exponentials with a better accuracy.

smp.exp(x).evalf(20) = 2.4848724344693696167

http://stackoverflow.com/questions/25181642/how-set-numpy-floating-point-accuracy

https://github.com/sympy/sympy/releases
http://docs.sympy.org/latest/index.html

[Aaron Meurer]

SymPy uses mpmath under the hood to get arbitrary precision. With
SymPy, you should either create a symbolic expression and evaluate it
with evalf and a precision, or start with something like
Float("1.4142", 100) (for 100 digits of precision) and manipulate it.

I'm actually seeing a lot of issues around this.

It seems there is a SymPy issue with creating this expression
symbolically https://github.com/sympy/sympy/issues/10963.

I should also point out that if you are taking a large power and the a
modulus, it's better to use pow(x, a, b) instead of x**a % b, since
the former is computed more efficiently using modular arithmetic.

With NumPy, I get wildly different answers, both of which are likely
wrong due to precision issues:

In [19]: (np.longdouble(1.4142)** 6000 )%400
Out[19]: 32.0

In [20]: pow(np.longdouble(1.4142), 6000, 400)
Out[20]: 1.1614417884357119756e+903

With SymPy's Float, pow() doesn't work
(https://github.com/sympy/sympy/issues/5715), but I get

In [21]: x = Float("1.4142", 100)

In [22]: x**6000 % 400
Out[22]: 272.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

mpmath unsurprisingly gives the same thing (by the way, sympy.mpmath
no longer works with SymPy 1.0, just use "import mpmath").

In [28]: x = mpmath.mpf("1.4142")

In [29]: x**6000 % 400
Out[29]: mpf('272.0')

(also it seems pow() doesn't work with mpf)

But, you can also use SymPy to convert 1.4142 to a rational number:

In [30]: nsimplify("1.4142")
Out[30]:
7071
────
5000

In [31]: 7071/5000
Out[31]: 1.4142

In [32]: x = nsimplify("1.4142")

In [33]: a = x**6000 % 400

In [34]: a.evalf()
Out[34]: 271.048181008631

I didn't print the value of a because it's a very large rational
number. Unless there's a bug somewhere in the rational number
calculations, or in the conversion of large rational numbers to
floats, I would wager this is the correct answer (at any rate, a would
have evaluated to 272 exactly if the answer were an exact integer).

Aaron Meurer

> --
> 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+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/e54f1042-193a-44db-88ec-8a4e13a8e44b%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Oscar Benjamin]

On 5 April 2016 at 01:56, Amy Valhausen <amy.vau...@gmail.com> wrote:
>
> import numpy
> (np.longdouble(1.4142)** 6000 )%400

...

>
> # The library mpmath is a good solution
>>>> import sympy as smp
>>>> mp = smp.mpmath
>
>>>> mp.mp.dps = 50 # Computation precision is 50 digits

50 digits is nowhere near enough. Think about what the number
1.4142**6000 is: it's a big number. If I wrote it out in digits it
might look something like:

5432321453425234...234234234.123235345345

The digits that contain the precision needed to get the remainder
modulo 400 begin just before the decimal point. So how many digits are
there before that? If d is the number of digits preceding the decimal
point then loosely:

10**d = 1.4142**6000

Which gives that

d = log10(1.4142**6000) = 6000*log10(1.4142) ~= 903

So if you want an answer that's good for m digits you'll need to use
about 900+m digits for the exponentiation:

In [1]: from sympy import mpmath

In [2]: mpmath.mp.dps = 950

In [3]: mpmath.mpf('1.4142') ** 6000 % 400
Out[3]: mpf('271.048181008630921815939109488389292518324580362344398848121124779167483584534976647550313880646779627157871790825164801629065802757168057723902165889739990234375')

BTW it's also possible to do this particular calculation easily enough
just with Python stdlib (and without numpy or sympy):

In [1]: from fractions import Fraction

In [2]: float(Fraction('1.4142')**6000 % 400)
Out[2]: 271.04818100863093


I though that it should be possible to easily do this with sympy
Floats but it doesn't seem to work:

In [1]: x = S(1.4142)

In [2]: x
Out[2]: 1.41420000000000

In [3]: x**6000
Out[3]: 1.16144178843571e+903

In [4]: x**6000 % 400
Out[4]: 32.0000000000000

This doesn't work because auto-evaluation of x**6000 destroyed all the
precision so I tried:

In [5]: Pow(x, 6000, evaluate=False)
Out[5]:
6000
1.4142

In [6]: Pow(x, 6000, evaluate=False) % 400
Out[6]: Mod(1.16144178843571e+903, 400)

Well that's annoying. It's gone ahead and auto-evaluated the x**6000
but created a Mod object instead of evaluating the modular division.

So let's disable evaluation of at the mod step as well:

In [10]: Mod(Pow(x, 6000, evaluate=False), 400, evaluate=False).evalf(50)
Out[10]:
⎛ 6000 ⎞
Mod⎝1.4142 , 400⎠

How do I get that to actually evaluate?


--
Oscar

[Aaron Meurer]

No, it's because the default precision is 15. I didn't realize you
needed a precision of 950. If you do what I did above with that you
get the right answer

In [6]: x = Float("1.4142", 950)

In [7]: x**6000%400
Out[7]:
271.04818100863092181593910948838929251832458036234439884812112477916748358453497664755031388064677962715787179082516480162906580275716805772390216588973999
023437500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000

>
> In [5]: Pow(x, 6000, evaluate=False)
> Out[5]:
> 6000
> 1.4142
>
> In [6]: Pow(x, 6000, evaluate=False) % 400
> Out[6]: Mod(1.16144178843571e+903, 400)
>
> Well that's annoying. It's gone ahead and auto-evaluated the x**6000
> but created a Mod object instead of evaluating the modular division.
>
> So let's disable evaluation of at the mod step as well:
>
> In [10]: Mod(Pow(x, 6000, evaluate=False), 400, evaluate=False).evalf(50)
> Out[10]:
> ⎛ 6000 ⎞
> Mod⎝1.4142 , 400⎠
>
> How do I get that to actually evaluate?

expr.doit() will evaluate the expression, but it evaluates it using
the prevision of x. For your x, that's 15, so you get the wrong answer
32. If you set x = Float("1.4142", 950) you get the right answer.

By the way, it seems the issue
https://github.com/sympy/sympy/issues/10963 is fixed in master (but
not in 1.0), so there you can do

In [3]: x = Symbol('x')

In [4]: expr = x**6000%400

In [7]: expr.evalf(950, subs={x:Float('1.4142', 950)})
Out[7]:
271.04818100863092181593910948838929251832458036234491201894546248073154256753666973207104639844775979650788048282365959816813533178603279338858556002378463
745117187500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000

But note even there you get the wrong answer unless you create 1.4142
as a Float with precision 950.

In [6]: expr.evalf(950, subs={x:1.4142})
Out[6]:
397.92898251746155210046958579278755870619830529784805636903899631078936617612319719713188467826854090314235414634273414074584762634145818083197809755802154
541015625000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000

I don't know if it should be considered a bug, but it's worth noting
that if you want SymPy to give the right precision in general you have
to start with Float objects that are set with the precision you need.
To me it feels like a bug because it negates the purpose of the evalf
precision argument.

Aaron Meurer

>
>
> --
> 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+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxRWe%2BDykG15%3Dks%2Bs04mF%2BgS2LmtRJ_kmAWQdwXNHAg5Bw%40mail.gmail.com.

[Oscar Benjamin]

On 5 April 2016 at 17:15, Aaron Meurer <asme...@gmail.com> wrote:
> On Tue, Apr 5, 2016 at 6:19 AM, Oscar Benjamin
> <oscar.j....@gmail.com> wrote:
>>
>> I though that it should be possible to easily do this with sympy
>> Floats but it doesn't seem to work:
>>
>> In [1]: x = S(1.4142)
>>
>> In [2]: x
>> Out[2]: 1.41420000000000
>>
>> In [3]: x**6000
>> Out[3]: 1.16144178843571e+903
>>
>> In [4]: x**6000 % 400
>> Out[4]: 32.0000000000000
>>
>> This doesn't work because auto-evaluation of x**6000 destroyed all the
>> precision so I tried:
>
> No, it's because the default precision is 15. I didn't realize you
> needed a precision of 950. If you do what I did above with that you
> get the right answer
>
> In [6]: x = Float("1.4142", 950)

...

It is because of the auto-evaluation in the sense that there were
enough digits to exactly represent all the input numbers but the
auto-evaluation cheated me from specifying how many digits to use for
the exponentiation with evalf later. I think that inexact
auto-evaluation should be disabled by default. Explicitly calling
.evalf() serves as a useful reminder of the numeric approximations
that have occurred in the calculation and should be required for this.

>> So let's disable evaluation of at the mod step as well:
>>
>> In [10]: Mod(Pow(x, 6000, evaluate=False), 400, evaluate=False).evalf(50)
>> Out[10]:
>> ⎛ 6000 ⎞
>> Mod⎝1.4142 , 400⎠
>>
>> How do I get that to actually evaluate?
>
> expr.doit() will evaluate the expression, but it evaluates it using
> the prevision of x. For your x, that's 15, so you get the wrong answer
> 32. If you set x = Float("1.4142", 950) you get the right answer.

Is there a reason for not evaluating this when evalf is called?

...

>
> I don't know if it should be considered a bug, but it's worth noting
> that if you want SymPy to give the right precision in general you have
> to start with Float objects that are set with the precision you need.
> To me it feels like a bug because it negates the purpose of the evalf
> precision argument.

Is there a coherent policy on float-handling in sympy?

My ideal would be:

1) All float objects are created exact (having the exact value of the
object passed in).
2) No inexact auto-evaluation.
3) .evalf() can be used to fully evaluate the expression with desired precision.
4) Ideally the precision argument to .evalf would be the precision of
the *output* rather than the internal precision of the intermediate
calculations

Currently 1) already occurs for decimal strings but Float(float)
rounds to 15 digits and you can explicitly force something impossible
as a ratio string: Float("1/3"). I think Float should be more like
decimal.Decimal here: all input arguments are treated as exact
regardless of precision etc. (and I don't see any good reason for
allowing Float("1/3"))

Without 2) it is impossible to achieve 3). If approximate
auto-evaluation can occur before calling .evalf then there's no way
for evalf to set the precision to be used for the auto-evaluation.

Obviously 4) is harder than current behaviour and perhaps impossible
in general but it is achievable for simple cases like in this thread.
From a user perspective it is definitely what is actually wanted and
much easier to understand.

--
Oscar

[Aaron Meurer]

No, that's also a bug. Unevaluated objects in SymPy tend to be pretty buggy.

>
> ...
>>
>> I don't know if it should be considered a bug, but it's worth noting
>> that if you want SymPy to give the right precision in general you have
>> to start with Float objects that are set with the precision you need.
>> To me it feels like a bug because it negates the purpose of the evalf
>> precision argument.
>
> Is there a coherent policy on float-handling in sympy?
>
> My ideal would be:
>
> 1) All float objects are created exact (having the exact value of the
> object passed in).
> 2) No inexact auto-evaluation.
> 3) .evalf() can be used to fully evaluate the expression with desired precision.
> 4) Ideally the precision argument to .evalf would be the precision of
> the *output* rather than the internal precision of the intermediate
> calculations

Can you clarify what you mean by "exact" here?

Note that there's no way to know what the input value of a float is.
That is, there's no way to write Float(0.2) (with no quotes) and have
it be treated as Float(2/10). The 0.2 object is converted to a Python
floating point by the time that Float sees it, and it's not a decimal:

In [49]: (0.2).as_integer_ratio()
Out[49]: (3602879701896397, 18014398509481984)

That's why Float allows string input (and it's the recommended way of
creating them).

With that being said, I don't think the fact that
(1.4142).as_integer_ratio() isn't (7071, 5000) is the problem here.
Float(1.4142) is indeed inexact compared to Float('1.4142'), but the
wrong answers from x**6000%400 come from lack of computing precision,
not lack of input accuracy.

>
> Currently 1) already occurs for decimal strings but Float(float)
> rounds to 15 digits and you can explicitly force something impossible
> as a ratio string: Float("1/3"). I think Float should be more like
> decimal.Decimal here: all input arguments are treated as exact
> regardless of precision etc. (and I don't see any good reason for
> allowing Float("1/3"))
>
> Without 2) it is impossible to achieve 3). If approximate
> auto-evaluation can occur before calling .evalf then there's no way
> for evalf to set the precision to be used for the auto-evaluation.
>
> Obviously 4) is harder than current behaviour and perhaps impossible
> in general but it is achievable for simple cases like in this thread.
> From a user perspective it is definitely what is actually wanted and
> much easier to understand.

I'm unclear how this works, because if you take my example above with
x = nsimplify("1.4142"), evalf() gave the right answer with the
default precision (15). That is, when everything in the expression is
a non-float, it gives the right answer. However, it seems that as soon
as an expression contains a Float, that Float must have whatever
precision you need set on it.

My knowledge of how this works internally is rather limited. Most of
what I know is based on my observations on how things work. I've CC'd
Fredrik Johansson. Hopefully he can give better insight into how
things should be working here.

Aaron Meurer

>
> --
> 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+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxRxmU9wrh__88Jw0ERkht3BPzFjo_cu3_iuvRdJdy_BrA%40mail.gmail.com.

[Isuru Fernando]

Hi,

I think the current way of representing Floats is reasonable.

Float internally keeps it in binary representation, so any non-terminating number in base 2 is truncated when stored as a Float. That's why there is a string constructor which can take in an exact decimal value and create a number in binary representation with precision mentioned.

Note that in this way only numbers representable as a sum of different powers of two have an exact value. Others can also be represented in an exact form if the base of the float is different than 2, but I think mpmath and many other libraries use only base 2.

If you want an exact form, best way would be to define it as a Rational and then you can do 1/3 as well which you can't if it's in decimal form. One other good thing about Rational is that `sin(Rational)` won't be evaluated unless SymPy knows the exact value, whereas `sin(Float)` would evaluate leading to error propagation even if it's something like `0.5` which represented exactly in binary representation.

Isuru Fernando

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6KhqcKoJgoK013m%2BotSRa1cAC%3DOxmtPVh%2B_6C_jGHbmzw%40mail.gmail.com.

[Oscar Benjamin]

On 5 April 2016 at 18:08, Aaron Meurer <asme...@gmail.com> wrote:
> On Tue, Apr 5, 2016 at 12:54 PM, Oscar Benjamin
> <oscar.j....@gmail.com> wrote:
>>
>>>
>>> I don't know if it should be considered a bug, but it's worth noting
>>> that if you want SymPy to give the right precision in general you have
>>> to start with Float objects that are set with the precision you need.
>>> To me it feels like a bug because it negates the purpose of the evalf
>>> precision argument.
>>
>> Is there a coherent policy on float-handling in sympy?
>>
>> My ideal would be:
>>
>> 1) All float objects are created exact (having the exact value of the
>> object passed in).
>> 2) No inexact auto-evaluation.
>> 3) .evalf() can be used to fully evaluate the expression with desired precision.
>> 4) Ideally the precision argument to .evalf would be the precision of
>> the *output* rather than the internal precision of the intermediate
>> calculations
>
> Can you clarify what you mean by "exact" here?
>
> Note that there's no way to know what the input value of a float is.
> That is, there's no way to write Float(0.2) (with no quotes) and have
> it be treated as Float(2/10). The 0.2 object is converted to a Python
> floating point by the time that Float sees it, and it's not a decimal:
>
> In [49]: (0.2).as_integer_ratio()
> Out[49]: (3602879701896397, 18014398509481984)

Passing a Python float into a sympy expression (I know I accidentally
did it earlier in the thread) is usually not going to do what is
wanted e.g. 0.1*x creates a number not truly equal to 0.1 and passes
it to x.__rmul__. The good fix for this is as you say to use a string.
However there are times when it would be good to pass in a float that
you have obtained from some other calculation and have sympy work with
it.

Currently Sympy will round an input float to 15d.p. and as a result
S(0.1) will result in an mpf which really does have a true value of
0.1. This is useful for novices but IMO just hides the binary float
problem a bit. The right solution is for users to understand that they
should be using S("0.1") or something.

If OTOH I received my float from some other source (rather than trying
to make a simple number like 0.1) then sympy is rounding the number
rather than taking it in exactly. I would prefer it in this case if
sympy would retain the exact value of the input number the same way
that Fraction and Decimal do:

In [19]: from fractions import Fraction

In [20]: Fraction(0.1)
Out[20]: Fraction(3602879701896397, 36028797018963968)

In [21]: from decimal import Decimal

In [22]: Decimal(0.1)
Out[22]:

Decimal('0.1000000000000000055511151231257827021181583404541015625')
In [23]: from sympy import Float

In [24]: Float(0.1)
Out[24]: 0.100000000000000

Both Fraction and Decimal always retain the exact input value. If you
want to round it then you need to do that explicitly. This is useful
because once you're numbers are converted to say Decimal then you can
use the calculation contexts to control exactly how the rounding
occurs (or to prevent any rounding) etc. If the rounding occurs
straight away at input then you cannot.

> That's why Float allows string input (and it's the recommended way of
> creating them).
>
> With that being said, I don't think the fact that
> (1.4142).as_integer_ratio() isn't (7071, 5000) is the problem here.
> Float(1.4142) is indeed inexact compared to Float('1.4142'), but the
> wrong answers from x**6000%400 come from lack of computing precision,
> not lack of input accuracy.

It is a separate but related issue to this. In this particular case
S(1.4142) does what a user may be intending it to do and creates
mpf("1.4142") so the issue occurs later.

>> Currently 1) already occurs for decimal strings but Float(float)
>> rounds to 15 digits and you can explicitly force something impossible
>> as a ratio string: Float("1/3"). I think Float should be more like
>> decimal.Decimal here: all input arguments are treated as exact
>> regardless of precision etc. (and I don't see any good reason for
>> allowing Float("1/3"))
>>
>> Without 2) it is impossible to achieve 3). If approximate
>> auto-evaluation can occur before calling .evalf then there's no way
>> for evalf to set the precision to be used for the auto-evaluation.
>>
>> Obviously 4) is harder than current behaviour and perhaps impossible
>> in general but it is achievable for simple cases like in this thread.
>> From a user perspective it is definitely what is actually wanted and
>> much easier to understand.
>
> I'm unclear how this works, because if you take my example above with
> x = nsimplify("1.4142"), evalf() gave the right answer with the
> default precision (15). That is, when everything in the expression is
> a non-float, it gives the right answer. However, it seems that as soon
> as an expression contains a Float, that Float must have whatever
> precision you need set on it.

I assume that you mean this:

In [4]: ((nsimplify("1.4142") ** 6000) % 400).evalf()
Out[4]: 271.048181008631

So let's pull that apart:

When you do nsimplify("1.4142") you get back a Rational.
Exponentiation that with an integer gives another Rational. The
expression is auto-evaluated to get the new rational but it is done
exactly. Likewise the modulo 400 is auto-evaluated but again it gives
an exact Rational.

Print this to see the massive numerator/denominator that are created:
nsimplify("1.4142") ** 6000 % 400

Since the final .evalf() is only called on a Rational which was itself
computed exactly and can easily give as many correct digits as you
desire.

However when you do this:

In [13]: Float("1.4142") ** 6000 % 400
Out[13]: 32.0000000000000

It works out differently. So what happens here is the expression
Float("1.4142") ** 6000 is auto-evaluated but only using the default
precision that is attached to the mpf Float (or maybe the context):

In [15]: Float("1.4142") ** 6000
Out[15]: 1.16144178843571e+903

So here we have auto-evaluation that has changed the numeric value of
the expression. This is not the same as simplifying 2**3 -> 8 or some
other simplification that is *exactly* correct because it has
introduced a numerical error. I think that this kind of
auto-evaluation should not occur by default.

So in the ideal scenario any auto-evaluation/simplification that would
not be exactly correct should not be applied. Then at the end when
calling .evalf it is possible to specify exactly the precision used
for *all* of the approximate steps in the calculation.

However as I said earlier a better general approach (although harder)
would be to have an API that allows to specify the accuracy of the
desired answer rather than the precision used for calculations. So in
this case it works back recursively. I have an expression:

s = Mod(Pow(Float("1.4142"), 6000), 400)

and I call

s.evalf(correct_digits=10)

and then this recursively calls evalf along the tree asking each time
for the number of correct digits needed in order to ensure that the
*end result* is correct to 10 digits.

--
Oscar

[Oscar Benjamin]

On 5 April 2016 at 18:21, Isuru Fernando <isu...@gmail.com> wrote:
>
> I think the current way of representing Floats is reasonable.
>
> Float internally keeps it in binary representation, so any non-terminating
> number in base 2 is truncated when stored as a Float. That's why there is a
> string constructor which can take in an exact decimal value and create a
> number in binary representation with precision mentioned.
>
> Note that in this way only numbers representable as a sum of different
> powers of two have an exact value. Others can also be represented in an
> exact form if the base of the float is different than 2, but I think mpmath
> and many other libraries use only base 2.

Some do and some use decimal. The advantage of decimal is being able
to exactly represent numbers input as decimal. Maybe the key here is
not to have Float automatically turn into an mpf in the first place if
it can't be done losslessly. Maybe a Float can remember the exact
input it was given e.g. "1.4142" and not try to turn that into a
binary representation until the precision is specified later.

> If you want an exact form, best way would be to define it as a Rational and
> then you can do 1/3 as well which you can't if it's in decimal form. One
> other good thing about Rational is that `sin(Rational)` won't be evaluated
> unless SymPy knows the exact value, whereas `sin(Float)` would evaluate
> leading to error propagation even if it's something like `0.5` which
> represented exactly in binary representation.

Exactly sin(Rational) won't be auto-evaluated because it's not clear
how many digits to use until someone later calls .evalf(). I think it
should be the same for sin(Float). It shouldn't really matter whether
the user does S("0.1") or S("1/10").

--
Oscar

[Oscar Benjamin]

I just realised I was confusing how mpmath works (thinking it was decimal rather than binary) so a few of the things I said were incorrect. When I'm next at a computer I'll start a new thread and clarify what I mean.

[Aaron Meurer]

You can see how mpmath works here
http://mpmath.org/doc/current/technical.html#representation-of-numbers.
It basically uses man*2**exp binary representation. For instance

In [249]: Float("0.1")._mpf_
Out[249]: (0, 3602879701896397, -55, 52)

The first entry is the sign (positive). The next two entries say that
0.1 is represented as

In [250]: 3602879701896397*2**-55
Out[250]: 0.1

and the last entry is the number of bits in 3602879701896397 (which is
redundant). Obviously, this number isn't exactly 1/10, which can't be
represented as <integer>/<power of 2>

In [251]: 3602879701896397*2**-S(55)
Out[251]:
3602879701896397
─────────────────
36028797018963968

You could also use a decimal representation man*10**exp, but using
base 2 is faster because you can use bit shifts for certain basic
operations (like addition). Even for Python longs, from what I can
tell, bit shifts are faster than multiplying by 2 or 10.

> Exactly sin(Rational) won't be auto-evaluated because it's not clear
how many digits to use until someone later calls .evalf(). I think it
should be the same for sin(Float). It shouldn't really matter whether
the user does S("0.1") or S("1/10").

I think most users don't care about precision beyond the default 15.
For them, auto-evaluating functions like sin at floating point values
is a useful feature. It's also easy to work around this, either
provide the precision when you create the float (like Float('0.1',
100)), or use a symbol, and substitute it at the end with evalf(100,
subs={x: '0.1'}).

The thing that bugs me is that, for instance in this example,
(x**6000%400).evalf(950, subs={x:"1.4142"}) doesn't work. You have to
provide the precision twice, like (x**6000%400).evalf(950,
subs={x:Float("1.4142", 950)}). So perhaps the default precision
shouldn't be 15, but rather some kind of "auto" flag that uses the
precision that's provided elsewhere (or 15 if no other precision is
provided). Or maybe evalf with the subs flag should just override the
precision for the float inputs with the precision argument of evalf.

Aaron Meurer

On Thu, Apr 7, 2016 at 11:19 AM, Oscar Benjamin

> --
> 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+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.
> To view this discussion on the web visit

> https://groups.google.com/d/msgid/sympy/CAHVvXxQNqEcLXYcFcVn0bW_9BNmO0O1xqLDrHmgfJURzdwv_DA%40mail.gmail.com.

[Oscar Benjamin]

On 7 April 2016 at 17:42, Aaron Meurer <asme...@gmail.com> wrote:
> You can see how mpmath works here
> http://mpmath.org/doc/current/technical.html#representation-of-numbers.
> It basically uses man*2**exp binary representation. For instance
>
> In [249]: Float("0.1")._mpf_
> Out[249]: (0, 3602879701896397, -55, 52)

Yeah this is what I found yesterday...

>> Exactly sin(Rational) won't be auto-evaluated because it's not clear
> how many digits to use until someone later calls .evalf(). I think it
> should be the same for sin(Float). It shouldn't really matter whether
> the user does S("0.1") or S("1/10").
>
> I think most users don't care about precision beyond the default 15.

Well actually this thread was started by a user who ended up with zero
correct digits because many more than 15 digits may be needed for
intermediate calculations in order to ensure that the end result is
correct to 15 digits.

Consider this:

>>> S("1e20 + 0.1 - 1e20")
0.0937500000000000000000

So we're using at least 15 d.p. for each quantity and calculation but
we get almost one correct digit in the output. I think this is
reasonable behaviour for a floating point library but not really from
a CAS like sympy (at least that's not why I use sympy!). I don't have
another CAS to hand but I just noticed that Wolfram Alpha apparently
gives both a wrong and a right answer for this
https://www.wolframalpha.com/input/?i=%281e20%2B0.1%29-1e20

Really when you do expr.evalf(15) what you want is not for sympy to
use 15 digits of precision internally and compute a result that has an
unspecified number of correct digits. What you want is for sympy to go
do what it takes to get an answer that you can trust to 15 digits.
Making that happen generally means going to a higher precision
internally.

For example if c = a*b and I want c to 15 digits I should compute a
and b to 16 digits, multiply them using 16 digits, and then round to
15 digits. Likewise if I want to compute a%b to 15 digits then I need
to compute a and b to something like ceil(log10(a/b))+15+1 digits
calculate the modulo and then round to 15 digits *at the end*. It is
possible to make something that will recursively compute the
expression tree in such a way that the end result of the computation
is trustable to the specified number of digits.

> For them, auto-evaluating functions like sin at floating point values
> is a useful feature. It's also easy to work around this, either
> provide the precision when you create the float (like Float('0.1',
> 100)), or use a symbol, and substitute it at the end with evalf(100,
> subs={x: '0.1'}).

I think I'll avoid using Floats altogether. I don't like the
possibility that any approximations occur implicitly and I don't
really feel like I know enough about how sympy works with floats to
trust it. Even the way floats are displayed is misleading:

>>> a = S("1e20")
>>> a
100000000000000000000.
>>> b = S("0.1")
>>> a
100000000000000000000.
>>> (a+b)
100000000000000000000.
>>> a == (a+b)
False
>>> str(a) == str(a+b)
True

> The thing that bugs me is that, for instance in this example,
> (x**6000%400).evalf(950, subs={x:"1.4142"}) doesn't work. You have to
> provide the precision twice, like (x**6000%400).evalf(950,
> subs={x:Float("1.4142", 950)}). So perhaps the default precision
> shouldn't be 15, but rather some kind of "auto" flag that uses the
> precision that's provided elsewhere (or 15 if no other precision is
> provided). Or maybe evalf with the subs flag should just override the
> precision for the float inputs with the precision argument of evalf.

The problem is that S("1.4142") throws away the string and stores a 53
bit binary approximation. The only way to increase the precision
retrospectively is to store the number in decimal/rational somewhere
(or just store the string itself).

--
Oscar

[Amy Valhausen]

Hi Aaron

I tried ;

import mpmath

x = mpmath.mpf("1.4142")

x**6000 % 400

but this returned "32" for me, not 272 like you got?

Also when I tried ;

>>> x = Float("1.4142", 100)

I got this error message ?

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
NameError: name 'Float' is not defined

I tried this ;
 

x = mpmath.mpf("1.4142")

x**6000 % 400

but again got "32" not 272 like you did?

You mentioned ;
 

But, you can also use SymPy to convert 1.4142 to a rational number:

 In [30]: nsimplify("1.4142")

I tried this but got the error message ;

>>> nsimplify("1.4142")
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
NameError: name 'nsimplify' is not defined

 
Could you also explain the benefit here for using 'nsimplify' ?

[Amy Valhausen]

[Amy Valhausen]

Oscar you mentioned below finding the # of digits using log
and the code ;

In [1]: from sympy import mpmath
In [2]: mpmath.mp.dps = 950
In [3]: mpmath.mpf('1.4142') ** 6000 % 400

How accurate of a return do you feel this will give?  Are we losing
any data, will the return be very accurate?


 

  10**d = 1.4142**6000
Which gives that
    d = log10(1.4142**6000) = 6000*log10(1.4142) ~= 903

So if you want an answer that's good for m digits you'll need to use
about 900+m digits for the exponentiation:

In [1]: from sympy import mpmath
In [2]: mpmath.mp.dps = 950
In [3]: mpmath.mpf('1.4142') ** 6000 % 400
Out[3]: mpf('271.048181008630921815939109488389

 

BTW it's also possible to do this particular calculation easily enough

[Amy Valhausen]

Hi Aaron,

In tried ;

import mpmath

x = Float("1.4142", 950)

x**6000%400

But got the same sort of error message again ;

>>> x = Float("1.4142", 950)

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
NameError: name 'Float' is not defined

Also, I thought I read recently that the float term is too small?
Some one said there are different kinds of floats ;
float, float96, float128??? 

I tried using "expr.doit()" but also got an error message?

[Amy Valhausen]

Hi, regarding the below, is it possible to increase precision for this past
(15) ?  Im a little bit confused.... are you saying that a float set to the default
precision will not return the accurate result but specifying the precision will?
Also you mentioned a non float expression will return the correct value, but
when I try ; x = nsimplify("1.4142"), evalf()

It returns only an error message?

 

[Oscar Benjamin]

On 11 April 2016 at 03:55, Amy Valhausen <amy.vau...@gmail.com> wrote:
> In tried ;
>
> import mpmath
> x = Float("1.4142", 950)
> x**6000%400
>
> But got the same sort of error message again ;
>
>>>> x = Float("1.4142", 950)
> Traceback (most recent call last):
> File "<stdin>", line 1, in <module>
> NameError: name 'Float' is not defined

You need to import Float and nsimplify from sympy with:

from sympy import Float, nsimplify

Or for these simple examples you could just do

from sympy import *

to import everything from sympy.

For your other questions about float accuracy see the responses from
Casevh and myself in the new thread that you started.

--
Oscar

[Chris Smith]

I see that this is an old post but for posterity...

The exact result is

 >>> exact = pow(Rational('1.4142'),6000,400)

It is exact because mpmath results *for Rationals* have unlimited precision.

To calculate to an arbitrary number of digits, evaluate this to a Float:

 >>> exact.n(17)
 271.04818100863092

Notice that the digit after the 309 rounds to 2 here but rounded to 3 in what you showed. This can be confirmed by evaluating to more digits

 >>> exact.n(22)
 271.0481810086309218159

You can also see how many digits for the input would have been needed to compute the output:

>>> [pow(Rational('1.4142').n(2**i),6000,400).n(17) for i in range(12)]
[1.0, 336.0, 224.0, 48.0, 368.0, 80.0, 192.0, 288.0, 192.0, 160.0, 271.04818100863092, 271.04818100863092]

As mentioned above, over 900 were needed.

Note also the failure of using evalf with subs to compute this sort of value:

>>> pow(x,6000,400).n(subs={x:'1.4142'})
272.000000000000

Note, too, the sensitivity to the input

>>> pow(Rational('1.41421'),6000,400).n(17)
7.7284031843133987

/c