summaryrefslogtreecommitdiff
path: root/src/pkg/math/big/rat.go
diff options
context:
space:
mode:
Diffstat (limited to 'src/pkg/math/big/rat.go')
-rw-r--r--src/pkg/math/big/rat.go269
1 files changed, 209 insertions, 60 deletions
diff --git a/src/pkg/math/big/rat.go b/src/pkg/math/big/rat.go
index 7bd83fc0f..3e6473d92 100644
--- a/src/pkg/math/big/rat.go
+++ b/src/pkg/math/big/rat.go
@@ -10,14 +10,17 @@ import (
"encoding/binary"
"errors"
"fmt"
+ "math"
"strings"
)
// A Rat represents a quotient a/b of arbitrary precision.
// The zero value for a Rat represents the value 0.
type Rat struct {
- a Int
- b nat // len(b) == 0 acts like b == 1
+ // To make zero values for Rat work w/o initialization,
+ // a zero value of b (len(b) == 0) acts like b == 1.
+ // a.neg determines the sign of the Rat, b.neg is ignored.
+ a, b Int
}
// NewRat creates a new Rat with numerator a and denominator b.
@@ -25,6 +28,156 @@ func NewRat(a, b int64) *Rat {
return new(Rat).SetFrac64(a, b)
}
+// SetFloat64 sets z to exactly f and returns z.
+// If f is not finite, SetFloat returns nil.
+func (z *Rat) SetFloat64(f float64) *Rat {
+ const expMask = 1<<11 - 1
+ bits := math.Float64bits(f)
+ mantissa := bits & (1<<52 - 1)
+ exp := int((bits >> 52) & expMask)
+ switch exp {
+ case expMask: // non-finite
+ return nil
+ case 0: // denormal
+ exp -= 1022
+ default: // normal
+ mantissa |= 1 << 52
+ exp -= 1023
+ }
+
+ shift := 52 - exp
+
+ // Optimisation (?): partially pre-normalise.
+ for mantissa&1 == 0 && shift > 0 {
+ mantissa >>= 1
+ shift--
+ }
+
+ z.a.SetUint64(mantissa)
+ z.a.neg = f < 0
+ z.b.Set(intOne)
+ if shift > 0 {
+ z.b.Lsh(&z.b, uint(shift))
+ } else {
+ z.a.Lsh(&z.a, uint(-shift))
+ }
+ return z.norm()
+}
+
+// isFinite reports whether f represents a finite rational value.
+// It is equivalent to !math.IsNan(f) && !math.IsInf(f, 0).
+func isFinite(f float64) bool {
+ return math.Abs(f) <= math.MaxFloat64
+}
+
+// low64 returns the least significant 64 bits of natural number z.
+func low64(z nat) uint64 {
+ if len(z) == 0 {
+ return 0
+ }
+ if _W == 32 && len(z) > 1 {
+ return uint64(z[1])<<32 | uint64(z[0])
+ }
+ return uint64(z[0])
+}
+
+// quotToFloat returns the non-negative IEEE 754 double-precision
+// value nearest to the quotient a/b, using round-to-even in halfway
+// cases. It does not mutate its arguments.
+// Preconditions: b is non-zero; a and b have no common factors.
+func quotToFloat(a, b nat) (f float64, exact bool) {
+ // TODO(adonovan): specialize common degenerate cases: 1.0, integers.
+ alen := a.bitLen()
+ if alen == 0 {
+ return 0, true
+ }
+ blen := b.bitLen()
+ if blen == 0 {
+ panic("division by zero")
+ }
+
+ // 1. Left-shift A or B such that quotient A/B is in [1<<53, 1<<55).
+ // (54 bits if A<B when they are left-aligned, 55 bits if A>=B.)
+ // This is 2 or 3 more than the float64 mantissa field width of 52:
+ // - the optional extra bit is shifted away in step 3 below.
+ // - the high-order 1 is omitted in float64 "normal" representation;
+ // - the low-order 1 will be used during rounding then discarded.
+ exp := alen - blen
+ var a2, b2 nat
+ a2 = a2.set(a)
+ b2 = b2.set(b)
+ if shift := 54 - exp; shift > 0 {
+ a2 = a2.shl(a2, uint(shift))
+ } else if shift < 0 {
+ b2 = b2.shl(b2, uint(-shift))
+ }
+
+ // 2. Compute quotient and remainder (q, r). NB: due to the
+ // extra shift, the low-order bit of q is logically the
+ // high-order bit of r.
+ var q nat
+ q, r := q.div(a2, a2, b2) // (recycle a2)
+ mantissa := low64(q)
+ haveRem := len(r) > 0 // mantissa&1 && !haveRem => remainder is exactly half
+
+ // 3. If quotient didn't fit in 54 bits, re-do division by b2<<1
+ // (in effect---we accomplish this incrementally).
+ if mantissa>>54 == 1 {
+ if mantissa&1 == 1 {
+ haveRem = true
+ }
+ mantissa >>= 1
+ exp++
+ }
+ if mantissa>>53 != 1 {
+ panic("expected exactly 54 bits of result")
+ }
+
+ // 4. Rounding.
+ if -1022-52 <= exp && exp <= -1022 {
+ // Denormal case; lose 'shift' bits of precision.
+ shift := uint64(-1022 - (exp - 1)) // [1..53)
+ lostbits := mantissa & (1<<shift - 1)
+ haveRem = haveRem || lostbits != 0
+ mantissa >>= shift
+ exp = -1023 + 2
+ }
+ // Round q using round-half-to-even.
+ exact = !haveRem
+ if mantissa&1 != 0 {
+ exact = false
+ if haveRem || mantissa&2 != 0 {
+ if mantissa++; mantissa >= 1<<54 {
+ // Complete rollover 11...1 => 100...0, so shift is safe
+ mantissa >>= 1
+ exp++
+ }
+ }
+ }
+ mantissa >>= 1 // discard rounding bit. Mantissa now scaled by 2^53.
+
+ f = math.Ldexp(float64(mantissa), exp-53)
+ if math.IsInf(f, 0) {
+ exact = false
+ }
+ return
+}
+
+// Float64 returns the nearest float64 value to z.
+// If z is exactly representable as a float64, Float64 returns exact=true.
+// If z is negative, so too is f, even if f==0.
+func (z *Rat) Float64() (f float64, exact bool) {
+ b := z.b.abs
+ if len(b) == 0 {
+ b = b.set(natOne) // materialize denominator
+ }
+ f, exact = quotToFloat(z.a.abs, b)
+ if z.a.neg {
+ f = -f
+ }
+ return
+}
+
// SetFrac sets z to a/b and returns z.
func (z *Rat) SetFrac(a, b *Int) *Rat {
z.a.neg = a.neg != b.neg
@@ -36,7 +189,7 @@ func (z *Rat) SetFrac(a, b *Int) *Rat {
babs = nat(nil).set(babs) // make a copy
}
z.a.abs = z.a.abs.set(a.abs)
- z.b = z.b.set(babs)
+ z.b.abs = z.b.abs.set(babs)
return z.norm()
}
@@ -50,21 +203,21 @@ func (z *Rat) SetFrac64(a, b int64) *Rat {
b = -b
z.a.neg = !z.a.neg
}
- z.b = z.b.setUint64(uint64(b))
+ z.b.abs = z.b.abs.setUint64(uint64(b))
return z.norm()
}
// SetInt sets z to x (by making a copy of x) and returns z.
func (z *Rat) SetInt(x *Int) *Rat {
z.a.Set(x)
- z.b = z.b.make(0)
+ z.b.abs = z.b.abs.make(0)
return z
}
// SetInt64 sets z to x and returns z.
func (z *Rat) SetInt64(x int64) *Rat {
z.a.SetInt64(x)
- z.b = z.b.make(0)
+ z.b.abs = z.b.abs.make(0)
return z
}
@@ -72,7 +225,7 @@ func (z *Rat) SetInt64(x int64) *Rat {
func (z *Rat) Set(x *Rat) *Rat {
if z != x {
z.a.Set(&x.a)
- z.b = z.b.set(x.b)
+ z.b.Set(&x.b)
}
return z
}
@@ -97,15 +250,15 @@ func (z *Rat) Inv(x *Rat) *Rat {
panic("division by zero")
}
z.Set(x)
- a := z.b
+ a := z.b.abs
if len(a) == 0 {
- a = a.setWord(1) // materialize numerator
+ a = a.set(natOne) // materialize numerator
}
b := z.a.abs
if b.cmp(natOne) == 0 {
b = b.make(0) // normalize denominator
}
- z.a.abs, z.b = a, b // sign doesn't change
+ z.a.abs, z.b.abs = a, b // sign doesn't change
return z
}
@@ -121,38 +274,26 @@ func (x *Rat) Sign() int {
// IsInt returns true if the denominator of x is 1.
func (x *Rat) IsInt() bool {
- return len(x.b) == 0 || x.b.cmp(natOne) == 0
+ return len(x.b.abs) == 0 || x.b.abs.cmp(natOne) == 0
}
// Num returns the numerator of x; it may be <= 0.
// The result is a reference to x's numerator; it
-// may change if a new value is assigned to x.
+// may change if a new value is assigned to x, and vice versa.
+// The sign of the numerator corresponds to the sign of x.
func (x *Rat) Num() *Int {
return &x.a
}
// Denom returns the denominator of x; it is always > 0.
// The result is a reference to x's denominator; it
-// may change if a new value is assigned to x.
+// may change if a new value is assigned to x, and vice versa.
func (x *Rat) Denom() *Int {
- if len(x.b) == 0 {
- return &Int{abs: nat{1}}
+ x.b.neg = false // the result is always >= 0
+ if len(x.b.abs) == 0 {
+ x.b.abs = x.b.abs.set(natOne) // materialize denominator
}
- return &Int{abs: x.b}
-}
-
-func gcd(x, y nat) nat {
- // Euclidean algorithm.
- var a, b nat
- a = a.set(x)
- b = b.set(y)
- for len(b) != 0 {
- var q, r nat
- _, r = q.div(r, a, b)
- a = b
- b = r
- }
- return a
+ return &x.b
}
func (z *Rat) norm() *Rat {
@@ -160,17 +301,25 @@ func (z *Rat) norm() *Rat {
case len(z.a.abs) == 0:
// z == 0 - normalize sign and denominator
z.a.neg = false
- z.b = z.b.make(0)
- case len(z.b) == 0:
+ z.b.abs = z.b.abs.make(0)
+ case len(z.b.abs) == 0:
// z is normalized int - nothing to do
- case z.b.cmp(natOne) == 0:
+ case z.b.abs.cmp(natOne) == 0:
// z is int - normalize denominator
- z.b = z.b.make(0)
+ z.b.abs = z.b.abs.make(0)
default:
- if f := gcd(z.a.abs, z.b); f.cmp(natOne) != 0 {
- z.a.abs, _ = z.a.abs.div(nil, z.a.abs, f)
- z.b, _ = z.b.div(nil, z.b, f)
+ neg := z.a.neg
+ z.a.neg = false
+ z.b.neg = false
+ if f := NewInt(0).binaryGCD(&z.a, &z.b); f.Cmp(intOne) != 0 {
+ z.a.abs, _ = z.a.abs.div(nil, z.a.abs, f.abs)
+ z.b.abs, _ = z.b.abs.div(nil, z.b.abs, f.abs)
+ if z.b.abs.cmp(natOne) == 0 {
+ // z is int - normalize denominator
+ z.b.abs = z.b.abs.make(0)
+ }
}
+ z.a.neg = neg
}
return z
}
@@ -207,31 +356,31 @@ func scaleDenom(x *Int, f nat) *Int {
// +1 if x > y
//
func (x *Rat) Cmp(y *Rat) int {
- return scaleDenom(&x.a, y.b).Cmp(scaleDenom(&y.a, x.b))
+ return scaleDenom(&x.a, y.b.abs).Cmp(scaleDenom(&y.a, x.b.abs))
}
// Add sets z to the sum x+y and returns z.
func (z *Rat) Add(x, y *Rat) *Rat {
- a1 := scaleDenom(&x.a, y.b)
- a2 := scaleDenom(&y.a, x.b)
+ a1 := scaleDenom(&x.a, y.b.abs)
+ a2 := scaleDenom(&y.a, x.b.abs)
z.a.Add(a1, a2)
- z.b = mulDenom(z.b, x.b, y.b)
+ z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
return z.norm()
}
// Sub sets z to the difference x-y and returns z.
func (z *Rat) Sub(x, y *Rat) *Rat {
- a1 := scaleDenom(&x.a, y.b)
- a2 := scaleDenom(&y.a, x.b)
+ a1 := scaleDenom(&x.a, y.b.abs)
+ a2 := scaleDenom(&y.a, x.b.abs)
z.a.Sub(a1, a2)
- z.b = mulDenom(z.b, x.b, y.b)
+ z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
return z.norm()
}
// Mul sets z to the product x*y and returns z.
func (z *Rat) Mul(x, y *Rat) *Rat {
z.a.Mul(&x.a, &y.a)
- z.b = mulDenom(z.b, x.b, y.b)
+ z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
return z.norm()
}
@@ -241,10 +390,10 @@ func (z *Rat) Quo(x, y *Rat) *Rat {
if len(y.a.abs) == 0 {
panic("division by zero")
}
- a := scaleDenom(&x.a, y.b)
- b := scaleDenom(&y.a, x.b)
+ a := scaleDenom(&x.a, y.b.abs)
+ b := scaleDenom(&y.a, x.b.abs)
z.a.abs = a.abs
- z.b = b.abs
+ z.b.abs = b.abs
z.a.neg = a.neg != b.neg
return z.norm()
}
@@ -286,7 +435,7 @@ func (z *Rat) SetString(s string) (*Rat, bool) {
}
s = s[sep+1:]
var err error
- if z.b, _, err = z.b.scan(strings.NewReader(s), 10); err != nil {
+ if z.b.abs, _, err = z.b.abs.scan(strings.NewReader(s), 10); err != nil {
return nil, false
}
return z.norm(), true
@@ -317,11 +466,11 @@ func (z *Rat) SetString(s string) (*Rat, bool) {
}
powTen := nat(nil).expNN(natTen, exp.abs, nil)
if exp.neg {
- z.b = powTen
+ z.b.abs = powTen
z.norm()
} else {
z.a.abs = z.a.abs.mul(z.a.abs, powTen)
- z.b = z.b.make(0)
+ z.b.abs = z.b.abs.make(0)
}
return z, true
@@ -330,8 +479,8 @@ func (z *Rat) SetString(s string) (*Rat, bool) {
// String returns a string representation of z in the form "a/b" (even if b == 1).
func (x *Rat) String() string {
s := "/1"
- if len(x.b) != 0 {
- s = "/" + x.b.decimalString()
+ if len(x.b.abs) != 0 {
+ s = "/" + x.b.abs.decimalString()
}
return x.a.String() + s
}
@@ -355,9 +504,9 @@ func (x *Rat) FloatString(prec int) string {
}
return s
}
- // x.b != 0
+ // x.b.abs != 0
- q, r := nat(nil).div(nat(nil), x.a.abs, x.b)
+ q, r := nat(nil).div(nat(nil), x.a.abs, x.b.abs)
p := natOne
if prec > 0 {
@@ -365,11 +514,11 @@ func (x *Rat) FloatString(prec int) string {
}
r = r.mul(r, p)
- r, r2 := r.div(nat(nil), r, x.b)
+ r, r2 := r.div(nat(nil), r, x.b.abs)
// see if we need to round up
r2 = r2.add(r2, r2)
- if x.b.cmp(r2) <= 0 {
+ if x.b.abs.cmp(r2) <= 0 {
r = r.add(r, natOne)
if r.cmp(p) >= 0 {
q = nat(nil).add(q, natOne)
@@ -396,8 +545,8 @@ const ratGobVersion byte = 1
// GobEncode implements the gob.GobEncoder interface.
func (x *Rat) GobEncode() ([]byte, error) {
- buf := make([]byte, 1+4+(len(x.a.abs)+len(x.b))*_S) // extra bytes for version and sign bit (1), and numerator length (4)
- i := x.b.bytes(buf)
+ buf := make([]byte, 1+4+(len(x.a.abs)+len(x.b.abs))*_S) // extra bytes for version and sign bit (1), and numerator length (4)
+ i := x.b.abs.bytes(buf)
j := x.a.abs.bytes(buf[0:i])
n := i - j
if int(uint32(n)) != n {
@@ -427,6 +576,6 @@ func (z *Rat) GobDecode(buf []byte) error {
i := j + binary.BigEndian.Uint32(buf[j-4:j])
z.a.neg = b&1 != 0
z.a.abs = z.a.abs.setBytes(buf[j:i])
- z.b = z.b.setBytes(buf[i:])
+ z.b.abs = z.b.abs.setBytes(buf[i:])
return nil
}
generated by cgit v1.2.3 (git 2.39.1) at 2025-02-09 08:39:23 +0000