Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP
Path: utzoo!mnetor!seismo!gatech!mit-eddie!uw-beaver!uw-june!uw-entropy!rag
From: rag@entropy.ms.washington.edu (David Ragozin)
Newsgroups: sci.math.symbolic,comp.bugs.4bsd,comp.bugs.sys5
Subject: Re: bug in bc and dc division and sqrt
Message-ID: <586@entropy.ms.washington.edu>
Date: Sat, 2-May-87 05:38:41 EDT
Article-I.D.: entropy.586
Posted: Sat May  2 05:38:41 1987
Date-Received: Sun, 3-May-87 05:39:18 EDT
References: <778@kulcs.UUCP>
Reply-To: rag@entropy.ms.washington.edu (David Ragozin)
Distribution: world
Organization: UW MathStat, Seattle
Lines: 54
Summary: bc and dc division,sqrt unreliable. arctan on BSD always wrong with fix
Xref: mnetor sci.math.symbolic:72 comp.bugs.4bsd:335 comp.bugs.sys5:113

Several notes about bc (dc) bugs.
A) A simplified version of someones earlier sqrt bug program also shows
   the division problem, with integer arithmetic:

	a=10^18-10^12+1      #(12 can be replaced by 9 through 15)
	b=a*a
	b/a
	a
	sqrt(b)

    All three displayed numbers should be equal, but are not. First
    two should be equal even with scale=0 as the division should be
    exact. Failure of third to be correct could be attributed to
    scale problems, but see B) below.

B)  Since sqrt is computed using newton's method, it involves divisions,
    so can be expected to fail if division fails.  Now for the
    "strange" part: If you prefix the previous program with 

	scale=2n+1   (i.e. any odd integer)
   
    a-sqrt(b) computes as 0, while if 

	scale=2n  (even scale)

    then  sqrt(b) computes as different from a!

    (b/a computes as unequal to a no matter what scale is set.)

C) Above results hold on Vax11/750,785's under BSD4.3 and on 
   3b2/400 under SysVRel3 and ATT UNIX PC under 3.5.

D) A TOTALLY different (trivial) bug occurs on all machines using
   /usr/lib/lib.b distributed by BSD labelled 4.1, dated 83/04/02.
   The code for a(x) (arctan of x) is missing braces, causing it
   to return pi/4 (up to 51 digits) independent of actual value of x.
   (This bug does not occur in the lib.b distributed with Sys V).
   A correct version of the offending lines:

	if(x==1) { 
		if(scale<52) {
return(.7853981633974483096156608458198757210492923498437764/1) }}

   simply adds the two sets of braces.

Unfortunately, until the division bug is fixed, all bc library
computations and sqrt's are suspect.  Does anyone know of a
domain of correct operation?

David Ragozin, Department of Mathematics, GN-50,
University of Washington, Seattle, WA. 98195. (206-543-1148)

rag@uw.ms.entropy.edu or
...!uw-beaver!uw-entropy!rag