[Search for users] [Overall Top Noters] [List of all Conferences] [Download this site]

Conference rusure::math

Title:Mathematics at DEC
Moderator:RUSURE::EDP
Created:Mon Feb 03 1986
Last Modified:Fri Jun 06 1997
Last Successful Update:Fri Jun 06 1997
Number of topics:2083
Total number of notes:14613

1061.0. "Just a llittle simple arithmetic" by NIZIAK::YARBROUGH (I PREFER PI) Tue Apr 18 1989 20:41

Here's an annoying problem from a recent Scientific Computing Text:

Evaluate, to six decimal places,

 c = 333.75*b^6 + a^2*(11*a^2*b^2-b^6-121*b^4-2) + 5.5*b^8 + a/(2*b)

at a = 77617, b = 33096.

Use your favorite programming language or calculator or whatever; you may 
rearrange the terms if you wish to make the calculation easier or more 
accurate. Use your own discretion in evaluating the exponentials.

The textbook reports that the IBM 4361, using 6, 14, or 28 hex digits
precision, in each case evaluates c to six places as 1.172604; how good is
that figure? 


Lynn Yarbrough 
T.RTitleUserPersonal
Name
DateLines
1061.1Wrong SignVAXRT::BRIDGEWATERTue Apr 18 1989 23:366
It's very inaccurate.  Spoiler follows


The correct answer is:    -54767/66192

Or, to 6 decimal places:  -0.827396
1061.2a/(2*b) for shortIOSG::CARLINDick Carlin IOSGWed Apr 19 1989 13:4810
    Interesting, the ALL-IN-1 desk calculator agrees with the 4361.
    However, even though I can't see the wicked glint in the eye of the
    person who orinally posed this, I assume that means they both fell into
    the same truncation trap.
    
    Or does everything except the last term really cancel out? It's not
    immediately obvious why.
    
    Dick
    
1061.3How high the moonNIZIAK::YARBROUGHI PREFER PIWed Apr 19 1989 14:2521
The correct result, -54767/66192, can be obtained either from MAPLE or VAX 
LISP (or other tools with extended arithmetic). The usual result is 
actually the value of the last term, a/2b, because everything else is 
astronomically large and cancels out within the precision available to 
machine floating-point arithmetic.

The expression was (first?) published in an IBM paper on ACRITH, the 
extended-precision arithmetic package; I saw it in the anthology 
*Scientific Computing with Automatic Result Verification*, publ. Springer-
Verlag 1988, Kerlisch & Stetter, Editors., in a paper by Fischer et al,
"Evaluation of Arithmetic Expressions with Guaranteed High Accuracy". A gem 
from this paper is an algorithm for evaluating such expressions in interval 
arithmetic, using the language PASCAL-SC (PASCAL for Sci. Comp.), which 
brackets the solution within +-10^-8 in three iterations.

Even as sophisticated a tool as MAPLE can crash badly on an ill-conditioned
problem like this, if not used wisely. MAPLE handles rational arithmetic
very well, but floating-point rather feebly, so typing the expression for c
as given results in .500000000*10^28! You have to express the constants as
33325/100 and 55/10, and a,b without decimal points, to get the full power
available. 
1061.4Roots of the problem.CADSYS::COOPERTopher CooperFri Apr 21 1989 16:5710
    Here's another one, not quite so spectacular but illustrating an
    important but little known caveat:
    
    Solve the following for X:
    
    	(X^2 + 1)*log((a+b)/(a+b^2)) + exp(a+b)*X = 0
    
    where a=11 and b = 1/2.
    
    					Topher
1061.5irational numbers will do .HANNAH::VESSALFri Apr 21 1989 17:333
    
    	Are we talking real numbers here ? the stated 
    	problem doesn't have a real number as an answer.
1061.6Where *do* those digits go?NIZIAK::YARBROUGHI PREFER PIFri Apr 21 1989 18:0714
>    Solve the following for X:
>    
>    	(X^2 + 1)*log((a+b)/(a+b^2)) + exp(a+b)*X = 0
>    
>    where a=11 and b = 1/2.

There are two solutions: approximately
	               -6
        X = -.222648+*10 
and	X = -4491386.777+

of which the first has only about 6 sig figs although based on 30-digit
arithmetic (thanks to MAPLE.) - VAX double precision yields 0. as the 
first answer. Big rounding errors.
1061.7No real answer...HANNAH::VESSALFri Apr 21 1989 18:579
    
          How did you test your answers fro correctness?
    
    	As X^2 + 1 >1 and (a+b)/(a+b^2) >1 which in turn makes
    	log  ((a+b)/(a+b^2)) > 0 and exp(a+b)*X >= 0 (always)
    	therefore no solution exist for the problem.
    
    			Hank
    
1061.8The answer may not be rational but its real.CADSYS::COOPERTopher CooperFri Apr 21 1989 20:1117
RE .7 (Hank)
    
    The requirement for there to be at least one real root for a quadratic
    is for:
    
    			B^2 >= 4*A*C
    
    In this case it translates to:
    
    			exp(a+b) >= 4*log((a+b)/(a+b^2))
    
    exp(a+b) is approximately 10^5, while log((a+b)/(a+b^2)) is about 10^-2
    (using common logrithms, which I used when I worked out the problem,
    but forgot to specify; although it doesn't really matter to the result)
    or about twice that using natural logs.  Obviously this qualifies.
    
    					Topher
1061.9sorry, but...NIZIAK::YARBROUGHI PREFER PIFri Apr 21 1989 20:163
... and exp(a+b)*X >= 0 (always) ...

Not if X<0. Perhaps you read that as exp((a+b)*X).
1061.10error on my part?HANNAH::VESSALFri Apr 21 1989 20:584
    My fault I read the expression incorrectly as Lynn stated in previous
    reply.
    
    
1061.11is it an integer or isn't it?CTCADM::ROTHIf you plant ice you'll harvest windFri Apr 21 1989 21:187
    Here's a couple of well-known examples:

	exp(pi*sqrt(163/9)) ~= 640320

	exp(pi*sqrt(163)) ~= 262537412640768744

    - Jim
1061.12AITG::DERAMODaniel V. {AITG,ZFC}:: D'EramoFri Apr 21 1989 23:0714
>>    Here's a couple of well-known examples:
>>
>>	exp(pi*sqrt(163/9)) ~= 640320

	In VAX LISP, (exp (* pi (sqrt (float 163/9 0.0l0)))) ==>
	640320.00000000060486373504901604L0.

>>	exp(pi*sqrt(163)) ~= 262537412640768744

	Also in VAX LISP, (exp (* pi (sqrt 163.0l0))) ==>
	2.62537412640768743999999999999251L17; subtracting same
	from 262537412640768744 left 7.49178497017055633477866649627686L-13.

	Dan