# Undeterministic numerical approximation

Let the following expression of x :

 f=1/152444172305856930250752000000*x^28 + 1/10888869450418352160768000000*x^27 + 1/15511210043330985984000000*x^25 + 1/310224200866619719680000*x^24 + 1/25852016738884976640000*x^23 + 1/51090942171709440000*x^21 + 1/1216451004088320000*x^20 + 1/121645100408832000*x^19 + 1/355687428096000*x^17 + 1/10461394944000*x^16 + 1/1307674368000*x^15 + 1/6227020800*x^13 + 1/239500800*x^12 + 1/39916800*x^11 + 1/362880*x^9 + 1/20160*x^8 + 1/5040*x^7 + 1/120*x^5 + 1/12*x^4 + 1/6*x^3 + x - cos(x) + 2 -exp(x) 

I noticed that evaluating that for the value x=numerical_approximation(10) leads to non-deterministic results (from the 12th decimal).

I guess that this is a "cancellation error" : the last term (-exp(x)) and the whole stuff before are both of the order of magnitude 22026.

I'm only interested in the first 3 or 4 decimals.

• QUESTION 1 How can I "force" sage to provide deterministic results ? Doing numerical_approx(...,prec=30) makes the job, but the "30" is quite heuristic; I'm not sure that this is a good approach.

• QUESTION 2 (not really sage-related) How should I unit-test a function that has undeterministic bahaviour like that ?

edit retag close merge delete

Sort by » oldest newest most voted

The following is a possibility to make the result predictible:

def f(x):
return \
( 1/152444172305856930250752000000*x^28
+ 1/10888869450418352160768000000*x^27
+ 1/15511210043330985984000000*x^25
+ 1/310224200866619719680000*x^24
+ 1/25852016738884976640000*x^23
+ 1/51090942171709440000*x^21
+ 1/1216451004088320000*x^20
+ 1/121645100408832000*x^19
+ 1/355687428096000*x^17
+ 1/10461394944000*x^16
+ 1/1307674368000*x^15
+ 1/6227020800*x^13
+ 1/239500800*x^12
+ 1/39916800*x^11
+ 1/362880*x^9
+ 1/20160*x^8
+ 1/5040*x^7
+ 1/120*x^5
+ 1/12*x^4
+ 1/6*x^3
+ x
- cos(x) + 2 -exp(x) )

R = RealField( 200 )
C = ComplexField( 200 )

f( C(10) )

f( R(10) )
f( R(10) )
f( R(10) )


Now the last four calls of f( C(10) ) and f( R(10) ) lead to...

-0.013417186790728922927885279844822282173830250128225204250732
-0.013417186790728922927885279844822282173830250128225204250732
-0.013417186790728922927885279844822282173830250128225204250732
-0.013417186790728922927885279844822282173830250128225204250732


We can also compare with...

sage: f(x) = exp(x) + cos(x) - 2
sage: p = taylor( f(x), x, 0, 28 )
sage: g(x) = sum( [ c*x^pow for c,pow in p.coefficients() ] ) - f(x)
sage: g(10)
-cos(10) - e^10 + 185630709328907872984/8427947351530707
sage: g(10.)
-0.0134171867909140
sage: g( RealField(100)(10) )
-0.013417186790728922927885322847
sage: g( RealField(200)(10) )
-0.013417186790728922927885279844822282173830250128225204254272
sage: g( RealField(400)(10) )
-0.0134171867907289229278852798448222821738302501282252042588478039342715326239015599771064232777651085304454608204788898295
sage: gp( g(10) )
-0.01341718679072892292788521428
sage: g(10).n()
-0.0134171867909140
sage: g(10).n( digits=100 )
-0.01341718679072892292788527984482228217383025012822520425884780393427153262390155997710642327776499594

more

Thanks for the answer. Since I am only interested in the firsts decimals (these are coordinates of of points on a A4 paper, in centimeters), I didn't though to increase the precision.

I'l have to tests that.