Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP
Path: utzoo!mnetor!seismo!mcvax!unido!gmdzi!woetzel
From: woetzel@gmdzi.UUCP (Gerd Woetzel)
Newsgroups: sci.math.symbolic,comp.bugs.4bsd,comp.bugs.sys5
Subject: Re: bug in bc and dc division and sqrt (and fix, hopefully)
Message-ID: <336@gmdzi.UUCP>
Date: Tue, 5-May-87 16:58:57 EDT
Article-I.D.: gmdzi.336
Posted: Tue May  5 16:58:57 1987
Date-Received: Thu, 7-May-87 07:14:08 EDT
References: <586@entropy.ms.washington.edu>
Distribution: world
Organization: GMD, Sankt Augustin, F.R.Germany
Lines: 119
Xref: mnetor sci.math.symbolic:73 comp.bugs.4bsd:341 comp.bugs.sys5:114

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"