golang/go

Switch branches/tags
Nothing to show
Fetching contributors…
Cannot retrieve contributors at this time
525 lines (475 sloc) 12.9 KB
 // Copyright 2010 The Go Authors. All rights reserved. // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. // This file implements multi-precision rational numbers. package big import ( "fmt" "math" ) // A Rat represents a quotient a/b of arbitrary precision. // The zero value for a Rat represents the value 0. // // Operations always take pointer arguments (*Rat) rather // than Rat values, and each unique Rat value requires // its own unique *Rat pointer. To "copy" a Rat value, // an existing (or newly allocated) Rat must be set to // a new value using the Rat.Set method; shallow copies // of Rats are not supported and may lead to errors. type Rat struct { // 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. 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 // Optimization (?): 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() } // quotToFloat32 returns the non-negative float32 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 quotToFloat32(a, b nat) (f float32, exact bool) { const ( // float size in bits Fsize = 32 // mantissa Msize = 23 Msize1 = Msize + 1 // incl. implicit 1 Msize2 = Msize1 + 1 // exponent Esize = Fsize - Msize1 Ebias = 1<<(Esize-1) - 1 Emin = 1 - Ebias Emax = Ebias ) // 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<= B). // This is 2 or 3 more than the float32 mantissa field width of Msize: // - the optional extra bit is shifted away in step 3 below. // - the high-order 1 is omitted in "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 := Msize2 - 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 := low32(q) haveRem := len(r) > 0 // mantissa&1 && !haveRem => remainder is exactly half // 3. If quotient didn't fit in Msize2 bits, redo division by b2<<1 // (in effect---we accomplish this incrementally). if mantissa>>Msize2 == 1 { if mantissa&1 == 1 { haveRem = true } mantissa >>= 1 exp++ } if mantissa>>Msize1 != 1 { panic(fmt.Sprintf("expected exactly %d bits of result", Msize2)) } // 4. Rounding. if Emin-Msize <= exp && exp <= Emin { // Denormal case; lose 'shift' bits of precision. shift := uint(Emin - (exp - 1)) // [1..Esize1) lostbits := mantissa & (1<>= shift exp = 2 - Ebias // == exp + shift } // Round q using round-half-to-even. exact = !haveRem if mantissa&1 != 0 { exact = false if haveRem || mantissa&2 != 0 { if mantissa++; mantissa >= 1< 100...0, so shift is safe mantissa >>= 1 exp++ } } } mantissa >>= 1 // discard rounding bit. Mantissa now scaled by 1<= B). // This is 2 or 3 more than the float64 mantissa field width of Msize: // - the optional extra bit is shifted away in step 3 below. // - the high-order 1 is omitted in "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 := Msize2 - 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 Msize2 bits, redo division by b2<<1 // (in effect---we accomplish this incrementally). if mantissa>>Msize2 == 1 { if mantissa&1 == 1 { haveRem = true } mantissa >>= 1 exp++ } if mantissa>>Msize1 != 1 { panic(fmt.Sprintf("expected exactly %d bits of result", Msize2)) } // 4. Rounding. if Emin-Msize <= exp && exp <= Emin { // Denormal case; lose 'shift' bits of precision. shift := uint(Emin - (exp - 1)) // [1..Esize1) lostbits := mantissa & (1<>= shift exp = 2 - Ebias // == exp + shift } // Round q using round-half-to-even. exact = !haveRem if mantissa&1 != 0 { exact = false if haveRem || mantissa&2 != 0 { if mantissa++; mantissa >= 1< 100...0, so shift is safe mantissa >>= 1 exp++ } } } mantissa >>= 1 // discard rounding bit. Mantissa now scaled by 1< 0 && !z.a.neg // 0 has no sign return z } // Inv sets z to 1/x and returns z. func (z *Rat) Inv(x *Rat) *Rat { if len(x.a.abs) == 0 { panic("division by zero") } z.Set(x) a := z.b.abs if len(a) == 0 { a = a.set(natOne) // materialize numerator } b := z.a.abs if b.cmp(natOne) == 0 { b = b[:0] // normalize denominator } z.a.abs, z.b.abs = a, b // sign doesn't change return z } // Sign returns: // // -1 if x < 0 // 0 if x == 0 // +1 if x > 0 // func (x *Rat) Sign() int { return x.a.Sign() } // IsInt reports whether the denominator of x is 1. func (x *Rat) IsInt() bool { 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, 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, and vice versa. func (x *Rat) Denom() *Int { 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 &x.b } func (z *Rat) norm() *Rat { switch { case len(z.a.abs) == 0: // z == 0 - normalize sign and denominator z.a.neg = false z.b.abs = z.b.abs[:0] case len(z.b.abs) == 0: // z is normalized int - nothing to do case z.b.abs.cmp(natOne) == 0: // z is int - normalize denominator z.b.abs = z.b.abs[:0] default: neg := z.a.neg z.a.neg = false z.b.neg = false if f := NewInt(0).lehmerGCD(nil, nil, &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[:0] } } z.a.neg = neg } return z } // mulDenom sets z to the denominator product x*y (by taking into // account that 0 values for x or y must be interpreted as 1) and // returns z. func mulDenom(z, x, y nat) nat { switch { case len(x) == 0: return z.set(y) case len(y) == 0: return z.set(x) } return z.mul(x, y) } // scaleDenom computes x*f. // If f == 0 (zero value of denominator), the result is (a copy of) x. func scaleDenom(x *Int, f nat) *Int { var z Int if len(f) == 0 { return z.Set(x) } z.abs = z.abs.mul(x.abs, f) z.neg = x.neg return &z } // Cmp compares x and y and returns: // // -1 if x < y // 0 if x == y // +1 if x > y // func (x *Rat) Cmp(y *Rat) int { 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.abs) a2 := scaleDenom(&y.a, x.b.abs) z.a.Add(a1, a2) 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.abs) a2 := scaleDenom(&y.a, x.b.abs) z.a.Sub(a1, a2) 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 { if x == y { // a squared Rat is positive and can't be reduced z.a.neg = false z.a.abs = z.a.abs.sqr(x.a.abs) z.b.abs = z.b.abs.sqr(x.b.abs) return z } z.a.Mul(&x.a, &y.a) z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs) return z.norm() } // Quo sets z to the quotient x/y and returns z. // If y == 0, a division-by-zero run-time panic occurs. func (z *Rat) Quo(x, y *Rat) *Rat { if len(y.a.abs) == 0 { panic("division by zero") } a := scaleDenom(&x.a, y.b.abs) b := scaleDenom(&y.a, x.b.abs) z.a.abs = a.abs z.b.abs = b.abs z.a.neg = a.neg != b.neg return z.norm() }