Answer exact results past the range of a double, such as 2^1200 (closes #16)

`2^1200` was refused although the calculator computed it exactly: writing a result went through a float64 and refused anything outside a double's range. An exact result under the 4096-bit limit is now written from its exact value, rounded to 17 significant digits past that range, and a fractional power of such a number is computed by first bringing its base into range with square roots taken from the exact value. Results computed in float64 must still be normal doubles; past the 4096-bit limit is still refused.

Disclosures: `1e-310` is answered again, reversing a call made in PR 11; a fractional power carries `math.Pow`'s rounding, which for a large base can reach the last digits shown.

Model: opus-5-5

internal/calc/calc.go

@@ -48,6 +48,18 @@ const (
 	plainLower = 1e-6
 )
 
+// A result past the normal range of a double is written to
+// significantDigits significant digits, the most the shortest form of a
+// double takes. It is rounded to them from a float of floatPrecision
+// bits, the bits a numerator or denominator can hold and 64 more for the
+// digits, so that the float rounds as the exact result would. The square
+// roots of a power's base past that range are taken in such a float too:
+// see nonNegativePower.
+const (
+	significantDigits = 17
+	floatPrecision    = bitLimit + 64
+)
+
 // The precedence of the binary operators: the higher, the tighter the
 // operator binds.
 const (
@@ -82,8 +94,7 @@ var (
 )
 
 // Evaluate computes an arithmetic expression and returns its result as
-// text: whole numbers without a decimal point, fractions in the
-// shortest form that reads back as the same float64.
+// text, written as format describes.
 func Evaluate(input string) (string, error) {
 	s := strings.TrimSpace(input)
 	if len(s) > MaxInputLength {
@@ -106,7 +117,7 @@ func Evaluate(input string) (string, error) {
 		return "", ErrNotArithmetic
 	}
 
-	return format(v)
+	return format(v), nil
 }
 
 // tokenize splits an expression into operators, parentheses and
@@ -395,11 +406,25 @@ func nonNegativePower(x, y, n constant.Value) (constant.Value, error) {
 
 	xf, _ := constant.Float64Val(x)
 	yf, _ := constant.Float64Val(y)
-	f := math.Pow(xf, yf)
 
-	// 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) {
+	// x^y is (√x)^(2y). An x outside the normal range of a double, such
+	// as 2^1200, would lose digits as a double, or all of them, so square
+	// roots taken from its exact value bring it into that range first. As
+	// x is between 2^-4096 and 2^4096 (see exact), three at most are
+	// needed.
+	r, _ := constant.Val(x).(*big.Rat)
+	root := new(big.Float).SetPrec(floatPrecision).SetRat(r)
+
+	for !normal(xf) {
+		root.Sqrt(root)
+		xf, _ = root.Float64()
+		yf *= 2
+	}
+
+	// x^y is not zero. If it is not a normal double, it has lost digits,
+	// or all of them.
+	f := math.Pow(xf, yf)
+	if !normal(f) {
 		return nil, ErrOutOfRange
 	}
 
@@ -459,26 +484,31 @@ func normal(f float64) bool {
 }
 
 // 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. 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) {
+// ordinary size is written exactly, digit for digit. Any other result in
+// the normal range of a double 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. Past that
+// range a double keeps fewer digits, or none (2^-1074 would be written
+// 5e-324, and 2^1024 is infinite), so such a result is written from its
+// exact value, to significantDigits.
+func format(v constant.Value) string {
 	f, _ := constant.Float64Val(v)
 	if constant.Sign(v) != 0 && !normal(f) {
-		return "", ErrOutOfRange
+		// Every number here is exact: see exact.
+		r, _ := constant.Val(v).(*big.Rat)
+
+		return new(big.Float).SetPrec(floatPrecision).SetRat(r).Text('g', significantDigits)
 	}
 
 	abs := math.Abs(f)
 
 	if i := constant.ToInt(v); i.Kind() == constant.Int && abs < plainUpper {
-		return i.ExactString(), nil
+		return i.ExactString()
 	}
 
 	if abs >= plainUpper || abs < plainLower {
-		return strconv.FormatFloat(f, 'g', -1, 64), nil
+		return strconv.FormatFloat(f, 'g', -1, 64)
 	}
 
-	return strconv.FormatFloat(f, 'f', -1, 64), nil
+	return strconv.FormatFloat(f, 'f', -1, 64)
 }