123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544 |
- // 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. At
- // the earliest opportunity (when an assignment to the Rat
- // is made), such uninitialized denominators are set to 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<<Msize1, 1<<(Msize2+1)
- // (Msize2 bits if A < B when they are left-aligned, Msize2+1 bits if A >= 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 - 1)
- haveRem = haveRem || lostbits != 0
- mantissa >>= 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<<Msize2 {
- // Complete rollover 11...1 => 100...0, so shift is safe
- mantissa >>= 1
- exp++
- }
- }
- }
- mantissa >>= 1 // discard rounding bit. Mantissa now scaled by 1<<Msize1.
- f = float32(math.Ldexp(float64(mantissa), exp-Msize1))
- if math.IsInf(float64(f), 0) {
- exact = false
- }
- return
- }
- // quotToFloat64 returns the non-negative float64 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 quotToFloat64(a, b nat) (f float64, exact bool) {
- const (
- // float size in bits
- Fsize = 64
- // mantissa
- Msize = 52
- 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<<Msize1, 1<<(Msize2+1)
- // (Msize2 bits if A < B when they are left-aligned, Msize2+1 bits if A >= 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 - 1)
- haveRem = haveRem || lostbits != 0
- mantissa >>= 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<<Msize2 {
- // Complete rollover 11...1 => 100...0, so shift is safe
- mantissa >>= 1
- exp++
- }
- }
- }
- mantissa >>= 1 // discard rounding bit. Mantissa now scaled by 1<<Msize1.
- f = math.Ldexp(float64(mantissa), exp-Msize1)
- if math.IsInf(f, 0) {
- exact = false
- }
- return
- }
- // Float32 returns the nearest float32 value for x and a bool indicating
- // whether f represents x exactly. If the magnitude of x is too large to
- // be represented by a float32, f is an infinity and exact is false.
- // The sign of f always matches the sign of x, even if f == 0.
- func (x *Rat) Float32() (f float32, exact bool) {
- b := x.b.abs
- if len(b) == 0 {
- b = natOne
- }
- f, exact = quotToFloat32(x.a.abs, b)
- if x.a.neg {
- f = -f
- }
- return
- }
- // Float64 returns the nearest float64 value for x and a bool indicating
- // whether f represents x exactly. If the magnitude of x is too large to
- // be represented by a float64, f is an infinity and exact is false.
- // The sign of f always matches the sign of x, even if f == 0.
- func (x *Rat) Float64() (f float64, exact bool) {
- b := x.b.abs
- if len(b) == 0 {
- b = natOne
- }
- f, exact = quotToFloat64(x.a.abs, b)
- if x.a.neg {
- f = -f
- }
- return
- }
- // SetFrac sets z to a/b and returns z.
- // If b == 0, SetFrac panics.
- func (z *Rat) SetFrac(a, b *Int) *Rat {
- z.a.neg = a.neg != b.neg
- babs := b.abs
- if len(babs) == 0 {
- panic("division by zero")
- }
- if &z.a == b || alias(z.a.abs, babs) {
- babs = nat(nil).set(babs) // make a copy
- }
- z.a.abs = z.a.abs.set(a.abs)
- z.b.abs = z.b.abs.set(babs)
- return z.norm()
- }
- // SetFrac64 sets z to a/b and returns z.
- // If b == 0, SetFrac64 panics.
- func (z *Rat) SetFrac64(a, b int64) *Rat {
- if b == 0 {
- panic("division by zero")
- }
- z.a.SetInt64(a)
- if b < 0 {
- b = -b
- z.a.neg = !z.a.neg
- }
- 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.abs = z.b.abs.setWord(1)
- return z
- }
- // SetInt64 sets z to x and returns z.
- func (z *Rat) SetInt64(x int64) *Rat {
- z.a.SetInt64(x)
- z.b.abs = z.b.abs.setWord(1)
- return z
- }
- // SetUint64 sets z to x and returns z.
- func (z *Rat) SetUint64(x uint64) *Rat {
- z.a.SetUint64(x)
- z.b.abs = z.b.abs.setWord(1)
- return z
- }
- // Set sets z to x (by making a copy of x) and returns z.
- func (z *Rat) Set(x *Rat) *Rat {
- if z != x {
- z.a.Set(&x.a)
- z.b.Set(&x.b)
- }
- if len(z.b.abs) == 0 {
- z.b.abs = z.b.abs.setWord(1)
- }
- return z
- }
- // Abs sets z to |x| (the absolute value of x) and returns z.
- func (z *Rat) Abs(x *Rat) *Rat {
- z.Set(x)
- z.a.neg = false
- return z
- }
- // Neg sets z to -x and returns z.
- func (z *Rat) Neg(x *Rat) *Rat {
- z.Set(x)
- z.a.neg = len(z.a.abs) > 0 && !z.a.neg // 0 has no sign
- return z
- }
- // Inv sets z to 1/x and returns z.
- // If x == 0, Inv panics.
- func (z *Rat) Inv(x *Rat) *Rat {
- if len(x.a.abs) == 0 {
- panic("division by zero")
- }
- z.Set(x)
- z.a.abs, z.b.abs = z.b.abs, z.a.abs
- 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, unless
- // x is an uninitialized (zero value) Rat, in which case
- // the result is a new Int of value 1. (To initialize x,
- // any operation that sets x will do, including x.Set(x).)
- // If 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 {
- // Note that x.b.neg is guaranteed false.
- if len(x.b.abs) == 0 {
- // Note: If this proves problematic, we could
- // panic instead and require the Rat to
- // be explicitly initialized.
- return &Int{abs: nat{1}}
- }
- 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
- fallthrough
- case len(z.b.abs) == 0:
- // z is integer; normalize denominator
- z.b.abs = z.b.abs.setWord(1)
- default:
- // z is fraction; normalize numerator and denominator
- 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)
- }
- 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 && len(y) == 0:
- return z.setWord(1)
- case len(x) == 0:
- return z.set(y)
- case len(y) == 0:
- return z.set(x)
- }
- return z.mul(x, y)
- }
- // scaleDenom sets z to the product x*f.
- // If f == 0 (zero value of denominator), z is set to (a copy of) x.
- func (z *Int) scaleDenom(x *Int, f nat) {
- if len(f) == 0 {
- z.Set(x)
- return
- }
- z.abs = z.abs.mul(x.abs, f)
- z.neg = x.neg
- }
- // 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 {
- var a, b Int
- a.scaleDenom(&x.a, y.b.abs)
- b.scaleDenom(&y.a, x.b.abs)
- return a.Cmp(&b)
- }
- // Add sets z to the sum x+y and returns z.
- func (z *Rat) Add(x, y *Rat) *Rat {
- var a1, a2 Int
- 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 {
- var a1, a2 Int
- 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 (no need to call norm())
- z.a.neg = false
- z.a.abs = z.a.abs.sqr(x.a.abs)
- if len(x.b.abs) == 0 {
- z.b.abs = z.b.abs.setWord(1)
- } else {
- 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, Quo panics.
- func (z *Rat) Quo(x, y *Rat) *Rat {
- if len(y.a.abs) == 0 {
- panic("division by zero")
- }
- var a, b Int
- 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()
- }
|