dc(1): fix input of non-decimal fractional numbers

Inputting fractional non-decimal numbers has never worked correctly in our
OpenBSD-derived dc(1). It truncates the input to a number of decimal places
equal to the number of hexadecimal (or whatever base) places given on the
input. That's unacceptable, because many numbers require more precision to
represent in base 10 than in their original bases.

Fix this bug by using as many decimal places as needed to represent the
input, up to the maximum of the global scale factor.

This has one mildly surprising side effect: the scale of a number entered in
non-decimal mode will no longer necessarily equal the number of hexadecimal
(or whatever base) places given on the input. I think that's an acceptable
behavior change, given that inputting fractional non-decimal numbers never
worked in the first place, and the man page doesn't specify whether trailing
zeros on the input should affect a number's scale.

PR:		206230
Reported by:	nibbana@gmx.us
Reviewed by:	pfg
Differential Revision:	https://reviews.freebsd.org/D13336

etc/mtree/BSD.tests.dist

@@ -654,6 +654,8 @@
         ..
         cut
         ..
+        dc
+        ..
         diff
         ..
         dirname

usr.bin/dc/Makefile

@@ -1,8 +1,13 @@
 #	$FreeBSD$
 #	$OpenBSD: Makefile,v 1.2 2006/11/26 11:31:09 deraadt Exp $
 
+.include <src.opts.mk>
+
 PROG=	dc
 SRCS=	dc.c bcode.c inout.c mem.c stack.c
 LIBADD=	crypto
 
+HAS_TESTS=
+SUBDIR.${MK_TESTS}+= tests
+
 .include <bsd.prog.mk>

usr.bin/dc/bcode.c

@@ -58,7 +58,6 @@ static __inline void	 unreadch(void);
 static __inline char	*readline(void);
 static __inline void	 src_free(void);
 
-static __inline u_int	 max(u_int, u_int);
 static u_long		 get_ulong(struct number *);
 
 static __inline void	 push_number(struct number *);
@@ -326,18 +325,12 @@ pbn(const char *str, const BIGNUM *n)
 
 #endif
 
-static __inline u_int
-max(u_int a, u_int b)
-{
-
-	return (a > b ? a : b);
-}
-
 static unsigned long factors[] = {
 	0, 10, 100, 1000, 10000, 100000, 1000000, 10000000,
 	100000000, 1000000000
 };
 
+/* Multiply n by 10^s */
 void
 scale_number(BIGNUM *n, int s)
 {
@@ -411,6 +404,7 @@ split_number(const struct number *n, BIGNUM *i, BIGNUM *f)
 	}
 }
 
+/* Change the scale of n to s.  Reducing scale may truncate the mantissa */
 void
 normalize(struct number *n, u_int s)
 {
@@ -1067,8 +1061,6 @@ static void
 bdiv(void)
 {
 	struct number *a, *b, *r;
-	BN_CTX *ctx;
-	u_int scale;
 
 	a = pop_number();
 	if (a == NULL)
@@ -1079,21 +1071,8 @@ bdiv(void)
 		return;
 	}
 
-	r = new_number();
+	r = div_number(b, a, bmachine.scale);
-	r->scale = bmachine.scale;
-	scale = max(a->scale, b->scale);
 
-	if (BN_is_zero(a->number))
-		warnx("divide by zero");
-	else {
-		normalize(a, scale);
-		normalize(b, scale + r->scale);
-
-		ctx = BN_CTX_new();
-		bn_checkp(ctx);
-		bn_check(BN_div(r->number, NULL, b->number, a->number, ctx));
-		BN_CTX_free(ctx);
-	}
 	push_number(r);
 	free_number(a);
 	free_number(b);
@@ -1681,7 +1660,7 @@ parse_number(void)
 
 	unreadch();
 	push_number(readnumber(&bmachine.readstack[bmachine.readsp],
-	    bmachine.ibase));
+	    bmachine.ibase, bmachine.scale));
 }
 
 static void

usr.bin/dc/bcode.h

@@ -95,3 +95,10 @@ void			negate(struct number *);
 void			split_number(const struct number *, BIGNUM *, BIGNUM *);
 void			bmul_number(struct number *, struct number *,
 			    struct number *, u_int scale);
+	
+static __inline u_int
+max(u_int a, u_int b)
+{
+
+	return (a > b ? a : b);
+}

usr.bin/dc/extern.h

@@ -24,7 +24,7 @@
 /* inout.c */
 void		 src_setstream(struct source *, FILE *);
 void		 src_setstring(struct source *, char *);
-struct number	*readnumber(struct source *, u_int);
+struct number	*readnumber(struct source *, u_int, u_int);
 void		 printnumber(FILE *, const struct number *, u_int);
 char		*read_string(struct source *);
 void		 print_value(FILE *, const struct value *, const char *, u_int);
@@ -33,6 +33,7 @@ void		 print_ascii(FILE *, const struct number *);
 /* mem.c */
 struct number	*new_number(void);
 void		 free_number(struct number *);
+struct number	*div_number(struct number *, struct number *, u_int scale);
 struct number	*dup_number(const struct number *);
 void		*bmalloc(size_t);
 void		*breallocarray(void *, size_t, size_t);

usr.bin/dc/inout.c

@@ -183,12 +183,12 @@ printwrap(FILE *f, const char *p)
 }
 
 struct number *
-readnumber(struct source *src, u_int base)
+readnumber(struct source *src, u_int base, u_int bscale)
 {
 	struct number *n;
 	BN_ULONG v;
-	u_int i;
 	int ch;
+	u_int iscale = 0;
 	bool dot = false, sign = false;
 
 	n = new_number();
@@ -213,15 +213,48 @@ readnumber(struct source *src, u_int base)
 			break;
 		}
 		if (dot)
-			n->scale++;
+			iscale++;
 
 		bn_check(BN_mul_word(n->number, base));
 		bn_check(BN_add_word(n->number, v));
 	}
-	if (base != 10) {
+	if (base == 10) {
-		scale_number(n->number, n->scale);
+		n->scale = iscale;
-		for (i = 0; i < n->scale; i++)
+	} else {
-			BN_div_word(n->number, base);
+		/* At this point, the desired result is n->number / base^iscale*/
+		struct number *quotient, *divisor, *_n;
+		BIGNUM *base_n, *exponent;
+		BN_CTX *ctx;
+
+		ctx = BN_CTX_new();
+		base_n = BN_new();
+		exponent = BN_new();
+		divisor = new_number();
+		bn_check(BN_zero(base_n));
+		bn_check(BN_zero(exponent));
+
+		bn_check(BN_add_word(base_n, base));
+		bn_check(BN_add_word(exponent, iscale));
+		bn_check(BN_exp(divisor->number, base_n, exponent, ctx));
+		divisor->scale = 0;
+		quotient = div_number(n, divisor, bscale);
+		_n = n;
+		n = quotient;
+
+		/* 
+		 * Trim off trailing zeros to yield the smallest scale without
+		 * loss of accuracy
+		 */
+		while ( n->scale > 0 &&
+			BN_mod_word(n->number, 10) == 0) {
+			normalize(n, n->scale - 1);
+		}
+
+		free_number(_n);
+		free_number(divisor);
+		BN_CTX_free(ctx);
+		BN_free(base_n);
+		BN_free(exponent);
 	}
 	if (sign)
 		negate(n);

usr.bin/dc/mem.c

@@ -22,6 +22,7 @@ __FBSDID("$FreeBSD$");
 #include <openssl/err.h>
 
 #include <err.h>
+#include <limits.h>
 #include <stdlib.h>
 #include <string.h>
 
@@ -48,6 +49,37 @@ free_number(struct number *n)
 	free(n);
 }
 
+/*
+ * Divide dividend by divisor, returning the result.  Retain bscale places of
+ * precision.
+ * The result must be freed when no longer in use
+ */
+struct number *
+div_number(struct number *dividend, struct number *divisor, u_int bscale)
+{
+	struct number *quotient;
+	BN_CTX *ctx;
+	u_int scale;
+
+	quotient = new_number();
+	quotient->scale = bscale;
+	scale = max(divisor->scale, dividend->scale);
+
+	if (BN_is_zero(divisor->number))
+		warnx("divide by zero");
+	else {
+		normalize(divisor, scale);
+		normalize(dividend, scale + quotient->scale);
+
+		ctx = BN_CTX_new();
+		bn_checkp(ctx);
+		bn_check(BN_div(quotient->number, NULL, dividend->number,
+				divisor->number, ctx));
+		BN_CTX_free(ctx);
+	}
+	return (quotient);
+}
+
 struct number *
 dup_number(const struct number *a)
 {

usr.bin/dc/tests/Makefile

@@ -0,0 +1,7 @@
+# $FreeBSD$
+
+PACKAGE=	tests
+
+ATF_TESTS_SH=	inout
+
+.include <bsd.test.mk>

usr.bin/dc/tests/inout.sh

@@ -0,0 +1,102 @@
+# SPDX-License-Identifier: BSD-2-Clause-FreeBSD
+#
+# Copyright (c) 2017 Alan Somers
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions
+# are met:
+# 1. Redistributions of source code must retain the above copyright
+#    notice, this list of conditions and the following disclaimer.
+# 2. Redistributions in binary form must reproduce the above copyright
+#    notice, this list of conditions and the following disclaimer in the
+#    documentation and/or other materials provided with the distribution.
+#
+# THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+# OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+# OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+# SUCH DAMAGE.
+#
+# $FreeBSD$
+
+atf_test_case base16_input
+base16_input_head()
+{
+	atf_set "descr" "Input hexadecimal numbers"
+}
+base16_input_body()
+{
+	cat > input.dc << EOF
+4k	# set scale to 4 decimal places
+16i	# switch to base 16
+0   p
+10  p
+1   p
+1.  p	# The '.' should have no effect
+1.0 p	# Unlike with decimal, should not change the result's scale
+.8  p	# Can input fractions
+# Check that we can input fractions that need more scale in base 10 than in 16
+# See PR 206230
+.1  p
+.10 p	# Result should be .0625, with scale=4
+.01 p	# Result should be truncated to scale=4
+8k	# Increase scale to 8 places
+.01 p	# Result should be exact again
+0.1 p	# Leading zeros are ignored
+00.1 p	# Leading zeros are ignored
+EOF
+dc input.dc > output.txt
+cat > expect.txt << EOF
+0
+16
+1
+1
+1
+.5
+.0625
+.0625
+.0039
+.00390625
+.0625
+.0625
+EOF
+	atf_check cmp expect.txt output.txt
+}
+
+atf_test_case base3_input
+base3_input_head()
+{
+	atf_set "descr" "Input ternary numbers"
+}
+base3_input_body()
+{
+	cat > input.dc << EOF
+4k	# 4 digits of precision
+3i	# Base 3 input
+0 p
+1 p
+10 p
+.1 p	# Repeating fractions get truncated
+EOF
+dc input.dc > output.txt
+cat > expect.txt << EOF
+0
+1
+3
+.3333
+EOF
+	atf_check cmp expect.txt output.txt
+}
+
+atf_init_test_cases()
+{
+	atf_add_test_case base16_input
+	atf_add_test_case base3_input
+}