Answer fractional powers of exact numbers past a double's range (closes #16)

A fractional power turned its base into a double first and refused a
base outside the normal range of a double, so (2^1200)^0.5, (2^1024)^0.5,
(2^-1200)^0.5 and 1e400^0.5 were refused although each answer is an
ordinary double. Such a base is now brought into that range by square
roots taken from its exact value in a big.Float, at most three under the
4096-bit limit, with the exponent doubled for each, and only a result
that is not a normal double is refused. The tests give these powers and
the edges of the range their values, and the bounded-work test covers a
base that needs three roots. The README's refusal sentence and example
are updated.

Model: opus-5-5

internal/calc/calc.go

@@ -52,7 +52,9 @@ const (
 // 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.
+// 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
@@ -404,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
 	}