Source: rat.go in package math/big

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.

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) {

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.

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 {

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) {

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.

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 {

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 {
	x.b.neg = false // the result is always >= 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 {

z == 0; normalize sign and denominator

		z.a.neg = false
		fallthrough

z is integer; normalize denominator

		z.b.abs = z.b.abs.setWord(1)

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 {

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()