rag@entropy.ms.washington.edu (David Ragozin) (05/02/87)
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!ragwoetzel@gmdzi.UUCP (Gerd Woetzel) (05/05/87)
The division errors of bc(1) are caused by the dc(1) command, which is used
indirectly by bc. Here are two fixes to dc.c made by Heinrich Schueth
(schue%gmdzi@unido.UUCP). The first one addresses an error in
the scale computation (line 60ff), the second is a correction to the
divide algorithm (line 665ff). Please send only credits to Heinrich, and
no flames to me. After the following diffs you'll find a little
shell script to test your implementation of bc/dc.
*** dc.c.OLD Tue May 5 13:27:28 1987
--- dc.c Tue May 5 14:13:14 1987
***************
*** 60,67
binop('*');
p = pop();
sunputc(p);
! savk = sk1+sk2;
! if(savk>k && savk>sk1 && savk>sk2){
sk = sk1;
if(sk<sk2)sk = sk2;
if(sk<k)sk = k;
--- 60,67 -----
binop('*');
p = pop();
sunputc(p);
! savk = n = sk1+sk2;
! if(n>k && n>sk1 && n>sk2){
sk = sk1;
if(sk<sk2)sk = sk2;
if(sk<k)sk = k;
***************
*** 65,71
sk = sk1;
if(sk<sk2)sk = sk2;
if(sk<k)sk = k;
! p = removc(p,savk-sk);
savk = sk;
}
sputc(p,savk);
--- 65,71 -----
sk = sk1;
if(sk<sk2)sk = sk2;
if(sk<k)sk = k;
! p = removc(p,n-sk);
savk = sk;
}
sputc(p,savk);
***************
*** 655,661
sbackc(p);
salterc(p,dig);
sbackc(p);
! if(--offset >= 0)divd->wt--;
}
if(divcarry != 0){
salterc(p,dig-1);
--- 655,668 -----
sbackc(p);
salterc(p,dig);
sbackc(p);
! if(--offset >= 0){
! if(d > 0){
! sbackc(divd);
! dd=sbackc(divd);
! salterc(divd,dd+100);
! }
! divd->wt--;
! }
}
if(divcarry != 0){
salterc(p,dig-1);
----------------- Test script --- Cut here -------------------
# Error in division. Found by schue%gmdzi@unido.UUCP (Heinrich Schueth)
# Fixed in dc.c line 665ff (see diffs above)"
echo "Test 1 (division): Result of expression is"
bc <<ENDSCRIPT1
scale=0
a=2799912455811601
b=99990001
(a/b)*b + (a%b) - a
ENDSCRIPT1
echo "Result should be"
echo "0"
echo ""
# Errors found by ronald@kulcs.UUCP (Ronald Cools)
# Same fix
echo "Test 2 (division): Results of divisions are"
bc <<ENDSCRIPT2
scale=45
a=999.999000000000000000000167001
1/a
10/a
100/a
1000/a
10000/a
ENDSCRIPT2
echo "Results should be"
echo ".001000001000001000001000000832999665998498997"
echo ".010000010000010000010000008329996659984989973"
echo "0.100000100000100000100000083299966599849899733"
echo "1.000001000001000001000000832999665998498997331"
echo "10.000010000010000010000008329996659984989973319"
echo ""
# Error at computation of scale. Found by Heinrich Schueth"
# Fix in dc.c line 60ff (see diffs above)"
echo "Test 3 (scale): Result of multiplication is"
bc <<ENDSCRIPT3
scale=64
a=1/9
a*a
ENDSCRIPT3
echo "Result should be"
echo ".0123456790123456790123456790123456790123456790123456790123456790"