restart: interface(rtablesize=10):
#
# Specify required number of digits numDig
# so the "best case error" will be
#
# 1.0*10^(-numDig)
#
# However the worst case error may be
# substantially greater than this, due to
# rounding at every stage of the calculation
# so define the "desired" error to be a couple
# of digits better, as in
#
numDig:= 10:
eps:= 1.0*10^(-numDig-2):
f:= 2*sin(x^2/2)-sin(1*x/2)^2:
df:= diff(f,x):
dff:= diff(f,x,x):
#
# Define a function whih computes the required
# quantity.Use the 'fulldigits' option to ensure
# that fsolve() is always using the current
# setting of Digits internally
#
do_dff:= ()-> eval
( dff,
x = fsolve
( df,
x = 1..2,
fulldigits
)
):
#
# Compute the desired answer using default
# Digits setting.
#
old:= do_dff():
#
# Increment 'Digits' and repeat the calculation
# until the difference between the value on the
# current setting of Digits and that on the previous
# setting is less than the desired error
#
for j from 1 do
Digits:= Digits+1:
new:= do_dff();
if abs(new-old) < eps
then break
else old:= new:
fi:
od:
#
# Out of idle curiosity, check the value of Digits
# where the error criterion was met, and the value
# obtained.
#
Digits;
new;
#
# Round the obtained value to the required number
# of Digits
#
evalf[numDig](new);