Skip to content

Instantly share code, notes, and snippets.

@ofan666
Created February 22, 2012 06:59
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save ofan666/1882702 to your computer and use it in GitHub Desktop.
Save ofan666/1882702 to your computer and use it in GitHub Desktop.
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