Created
February 22, 2012 06:59
-
-
Save ofan666/1882702 to your computer and use it in GitHub Desktop.
Arbitrary Precision using Python Mpmath part1 - http://ofan666.blogspot.com/2011/06/arbitrary-precision-using-python-mpmath.html
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import sys | |
# using mpmath for arbitrary precision, | |
# sympy for showing the equation in pretty way | |
if sys.platform == 'linux2': | |
import mpmath as mp, sympy as sm | |
elif sys.platform == 'win32': | |
import sympy as sm | |
mp = sm.mpmath | |
# just for write shorter code | |
f = mp.mpf; ms = mp.nstr; R = sm.Rational; pr = sm.pprint | |
sm.var('x y') | |
eqn = (333.75 - x**2)* y**6 + \ | |
x**2 *((11)* x**2 * y**2 - (121) * y**4 - 2)+\ | |
5.5 * y**8 + x/((2)*y) | |
print 'Equation for testing precision'; pr(eqn) | |
print 'with x=77617, y=33096' | |
print '\n\nSimplified equation using sympy or by hand by change to rational value' | |
sol = (R('333.75') - x**2)* y**6 + \ | |
x**2 *((11)* x**2 * y**2 - (121) * y**4 - 2)+\ | |
R('5.5') * y**8 + x/((2)*y) | |
pr(sol) | |
print 'Subtitute x:77617, y:33096 yield' | |
pr(sol.subs({x:77617, y:33096})) | |
# since the source(Scipy2008_proceedings.pdf) only told 6 decimal digits | |
# so i just evaluate to 6 decimal digits | |
pr(sol.subs({x:77617, y:33096}).evalf(6)) | |
print '\n\nDefault', mp.mp | |
# same function, but implemented with mpf from mpmath | |
def f_1(x,y): | |
return ((f('333.75') - x**2)* y**6 + \ | |
x**2 *((11)* x**2 * y**2 - (121) * y**4 - 2)+ \ | |
f('5.5') * y**8 + x/((2)*y)) | |
# calculate -54767/66192 using mpf from mpmath, by converting | |
# to str 6 decimal places for comparing to find the significant | |
# bit prec setting | |
exact = ms(-f(54767) / f(66192), n=6) # at this prec is equal 53 (default) | |
print '\nSearching prec value, to give result', exact | |
# calc -> calculate ; cpr -> for comparing with previosly calculated | |
calc = cpr = 0 | |
# loop for searching prec setting that matched with reference value -0.827396' | |
# by calculating f_1 function | |
while calc != exact: | |
mp.mp.prec +=1 | |
calc = mp.nstr(f_1(mp.mpf(77617), mp.mpf(33096)), n=6) | |
if cpr != calc: | |
cpr = calc | |
# uncomment below to find out when the calculated is changed | |
# than previous when prec is changed | |
#print cpr, mp.mp | |
print 'Found Significant', mp.mp | |
print 'Result using above setting =', calc | |
print f_1(mp.mpf(77617), mp.mpf(33096)) | |
print 'mpmath backend', mp.libmp.BACKEND #show mpmath backend |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment