How do I use sympy.solve() with greater precision than 15?

[Damon Turney]

Hi Forum,

Help!  sympy.solve() always uses precision of 15 no matter what I try.

Recreated problem results:

In [164]: import sympy                                                                                                                                                                                           

In [165]: w,x,y,z=sympy.symbols('w,x,y,z')                                                                                                                                                                      

In [166]: eqn=sympy.Eq(sympy.Float('1/55500',45)*w,sympy.Float('96500.0')*y)                                                                                                                                    

In [167]: eqn.evalf(45)                                                                                                                                                                                         

Out[167]: Eq(0.000018018018018018018018018018018018018018018018*w, 96500.0*y)

In [168]: sympy.solve(eqn,y)[0].evalf(45)                                                                                                                                                                       

Out[168]: 0.000000000186715212621948369903318027445619613213789734*w

In [169]: sympy.solve(eqn.evalf(45),y)[0].evalf(45)                                                                                                                                                             

Out[169]: 0.000000000186715212621948369903318027445619613213789734*w

The more accurate coefficient (to 45 digit accuracy) is:

In [152]: sympy.Float('1/55500',45)/sympy.Float('96500.0',45).evalf(45)                                                                                                                                          

Out[152]: 0.000000000186715212621948373243710031274798114176352518

In [174]: sympy.Float('1',45)/sympy.Float('55500',45)/sympy.Float('96500.0',45).evalf(45)                                                                                                                        

Out[174]: 0.000000000186715212621948373243710031274798114176352518

So you see, the coefficient that results from sympy.solve: 

0.000000000186715212621948369903318027445619613213789734

is accurate to only 15 digits of precision.

Clearly sympy.solve is using only 15 digit precision, even though I'm specifying 45 digit precision. 

How can I force sympy.solve to use 45 digit precision when it solves for y?

[Damon Turney]

Update:  An angel on GitHub told me that I can overcome this problem by using solveset() instead of solve().

But why does solve() not work as expected?

[Aaron Meurer]

solve() seems to have a bug here.

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 view this discussion on the web visit https://groups.google.com/d/msgid/sympy/996dada4-fc41-49a7-ae08-cd8577bca974%40googlegroups.com.