Refuse numbers a double cannot hold; exact powers up to the limit (closes #3)

A power computed in float64 is refused unless its base and its result
are normal doubles, and so is a result other than zero that is not one:
below about 2.2e-308 a double keeps fewer digits, and at 0 or infinity
none. A literal that go/constant reads as 0 though it is not, such as
1e-999999999, is refused too.

A whole exponent up to 4096 either way is computed exactly and kept
when go/constant holds the result exactly, so 2^-1400 is exact.

ErrTooLarge becomes ErrOutOfRange, and its reply, "That needs a number
too large or too small for me.", is true of both.

Model: opus-5-5

internal/calc/calc.go

@@ -6,8 +6,8 @@
 // go/constant, which does exact rational arithmetic: 5 * 5/2 is exactly
 // 12.5, and 0.1 + 0.2 is exactly 0.3, so a result carries no binary
 // floating point noise until the moment it is formatted. A power is the
-// exception: one with a fractional exponent, or too large to compute
-// exactly, is computed in float64.
+// exception: one with a fractional exponent, or whose result go/constant
+// cannot hold exactly, is computed in float64.
 package calc
 
 import (
@@ -21,16 +21,20 @@ import (
 	"strings"
 )
 
-// MaxInputLength caps an expression, in bytes, and maxExactPowerBits
-// caps a power, so that a message cannot make the bot do unbounded work.
+// MaxInputLength caps an expression, in bytes, and maxExactExponent caps
+// a power computed exactly, so that a message cannot make the bot do
+// unbounded work.
 const MaxInputLength = 256
 
-// maxExactPowerBits caps a power computed exactly: its numerator and
-// denominator together have at most this many bits, estimated before
-// multiplying as the exponent times the bits in the base's numerator and
-// denominator. A larger power is computed in float64, whose cost does not
-// grow with it.
-const maxExactPowerBits = 4096
+// maxExactExponent is the largest exponent, either way, of a power
+// computed exactly. Past it, x^n has a numerator or denominator of more
+// than 4096 bits, which go/constant holds only rounded, unless x is 0 or
+// 1, and float64 computes those exactly.
+const maxExactExponent = 4096
+
+// smallestNormal is the smallest positive normal double, about 2.2e-308.
+// Below it a double keeps fewer digits, down to one.
+const smallestNormal = 0x1p-1022
 
 // Results of magnitude plainUpper or more are written in exponent form
 // (1e+21 rather than twenty-two digits), and so are fractions smaller
@@ -53,7 +57,7 @@ var (
 	ErrTooLong        = errors.New("expression too long")
 	ErrNotArithmetic  = errors.New("not an arithmetic expression")
 	ErrDivisionByZero = errors.New("division by zero")
-	ErrTooLarge       = errors.New("result too large")
+	ErrOutOfRange     = errors.New("number too large or too small")
 	ErrNoRealResult   = errors.New("no real result")
 )
 
@@ -245,10 +249,16 @@ func number(tok string) (constant.Value, error) {
 	// leading zero would make it octal.
 	v := constant.MakeFromLiteral(tok, token.FLOAT, 0)
 
-	// The syntax was checked above, so Unknown here means the exponent
-	// overflowed. A literal such as 1e1300 is held rounded: see rounded.
-	if v.Kind() == constant.Unknown || rounded(v) {
-		return nil, ErrTooLarge
+	// A literal such as 1e1300 or 1e-1300 is held rounded: see exact.
+	if !exact(v) {
+		return nil, ErrOutOfRange
+	}
+
+	// One too small even to be held rounded, such as 1e-999999999, is
+	// read as 0.
+	mantissa, _, _ := strings.Cut(strings.ToLower(tok), "e")
+	if constant.Sign(v) == 0 && strings.ContainsAny(mantissa, "123456789") {
+		return nil, ErrOutOfRange
 	}
 
 	return v, nil
@@ -282,10 +292,8 @@ func apply(x constant.Value, op string, y constant.Value) (constant.Value, error
 		return nil, err
 	}
 
-	// go/constant represents an overflow to infinity as Unknown. A
-	// rounded number is refused too: see rounded.
-	if v.Kind() == constant.Unknown || rounded(v) {
-		return nil, ErrTooLarge
+	if !exact(v) {
+		return nil, ErrOutOfRange
 	}
 
 	return v, nil
@@ -312,8 +320,8 @@ func modulo(x, y constant.Value) (constant.Value, error) {
 
 	// The fractional part of a rounded quotient, and so the remainder,
 	// would be wrong.
-	if rounded(q) {
-		return nil, ErrTooLarge
+	if !exact(q) {
+		return nil, ErrOutOfRange
 	}
 
 	// x % y is y times the fractional part of x/y, which is at least 0
@@ -345,13 +353,16 @@ func power(x, y constant.Value) (constant.Value, error) {
 	case constant.Sign(x) == 0 && constant.Sign(y) < 0:
 		return nil, ErrDivisionByZero
 	case constant.Sign(x) >= 0:
-		return nonNegativePower(x, y, n), nil
+		return nonNegativePower(x, y, n)
 	case n.Kind() != constant.Int:
 		return nil, ErrNoRealResult
 	}
 
 	// x is negative and n whole: x^n is (-x)^n, negated if n is odd.
-	v := nonNegativePower(constant.UnaryOp(token.SUB, x, 0), y, n)
+	v, err := nonNegativePower(constant.UnaryOp(token.SUB, x, 0), y, n)
+	if err != nil {
+		return nil, err
+	}
 
 	odd := constant.BinaryOp(n, token.AND, constant.MakeInt64(1))
 	if constant.Sign(odd) != 0 {
@@ -362,36 +373,38 @@ func power(x, y constant.Value) (constant.Value, error) {
 }
 
 // nonNegativePower computes x^y for x of at least zero, and y not below
-// zero if x is zero: exactly if y is a whole number n and the result
-// fits in maxExactPowerBits, otherwise in float64.
-func nonNegativePower(x, y, n constant.Value) constant.Value {
+// zero if x is zero: exactly if y is a whole number n and go/constant
+// holds the result exactly, otherwise in float64.
+func nonNegativePower(x, y, n constant.Value) (constant.Value, error) {
 	e, ok := constant.Int64Val(n)
-	if ok && exactPowerFits(x, e) {
-		return exactPower(x, e)
+	if ok && -maxExactExponent <= e && e <= maxExactExponent {
+		v := exactPower(x, e)
+		if exact(v) {
+			return v, nil
+		}
+	}
+
+	// y is above zero here if x is zero.
+	if constant.Sign(x) == 0 {
+		return x, nil
 	}
 
 	xf, _ := constant.Float64Val(x)
 	yf, _ := constant.Float64Val(y)
+	f := math.Pow(xf, yf)
 
-	// An infinite result becomes Unknown, which apply refuses as too
-	// large.
-	return constant.MakeFloat64(math.Pow(xf, yf))
-}
-
-// exactPowerFits reports whether x^e fits in maxExactPowerBits.
-func exactPowerFits(x constant.Value, e int64) bool {
-	// Checked first so that the product below cannot overflow.
-	if e < -maxExactPowerBits || e > maxExactPowerBits {
-		return false
+	// Neither x nor x^y is zero. If either is not a normal double, it
+	// has lost digits, or all of them.
+	if !normal(xf) || !normal(f) {
+		return nil, ErrOutOfRange
 	}
 
-	bits := int64(constant.BitLen(constant.Num(x)) + constant.BitLen(constant.Denom(x)))
-
-	return bits*max(e, -e) <= maxExactPowerBits
+	return constant.MakeFloat64(f), nil
 }
 
 // exactPower computes x^e by repeated squaring. x is not zero if e is
-// negative.
+// negative. Each step's numbers stay small: go/constant holds one whose
+// numerator or denominator reaches 4096 bits as a 512-bit float.
 func exactPower(x constant.Value, e int64) constant.Value {
 	result := constant.MakeInt64(1)
 
@@ -410,28 +423,41 @@ func exactPower(x constant.Value, e int64) constant.Value {
 	return result
 }
 
-// rounded reports whether go/constant holds v rounded. It holds a
-// number exactly, as a fraction, only while the numerator and the
-// denominator each stay under 4096 bits; past that, and for a literal of
-// that size, it holds a 512-bit float. Such a number is refused as too
-// large wherever it appears: a sum can lose the answer entirely
-// (7^1000*7^1000 + 5 - 7^1000*7^1000 would be 0), and a remainder, or
-// whether an exponent is whole or odd, cannot be read from one.
-func rounded(v constant.Value) bool {
-	_, isFloat := constant.Val(v).(*big.Float)
+// exact reports whether go/constant holds v exactly. It holds a number
+// as a fraction until its numerator or denominator reaches 4096 bits,
+// then as a 512-bit float, and past that float's range as Unknown. A
+// number not held exactly is refused wherever it appears: a sum can lose
+// the answer entirely (7^1000*7^1000 + 5 - 7^1000*7^1000 would be 0), and
+// a remainder, or whether an exponent is whole or odd, cannot be read
+// from one.
+func exact(v constant.Value) bool {
+	switch constant.Val(v).(type) {
+	case int64, *big.Int, *big.Rat:
+		return true
+	default:
+		return false
+	}
+}
 
-	return isFloat
+// normal reports whether f is a normal double, finite and at least
+// smallestNormal in magnitude: a number other than zero keeps all of a
+// double's digits only as one.
+func normal(f float64) bool {
+	abs := math.Abs(f)
+
+	return abs >= smallestNormal && abs <= math.MaxFloat64
 }
 
 // format writes a result for a person to read. A whole number of
 // ordinary size is written exactly, digit for digit; anything else goes
 // through float64, whose shortest round-trip form is free of the noise
 // (0.30000000000000004) that printing a binary fraction to a fixed
-// precision produces.
+// precision produces. A result that is not zero must therefore be a
+// normal double: 2^-1074 would be written 5e-324.
 func format(v constant.Value) (string, error) {
 	f, _ := constant.Float64Val(v)
-	if math.IsInf(f, 0) || math.IsNaN(f) {
-		return "", ErrTooLarge
+	if constant.Sign(v) != 0 && !normal(f) {
+		return "", ErrOutOfRange
 	}
 
 	abs := math.Abs(f)