Source file src/pkg/math/big/rat.go

     1	// Copyright 2010 The Go Authors. All rights reserved.
     2	// Use of this source code is governed by a BSD-style
     3	// license that can be found in the LICENSE file.
     4	
     5	// This file implements multi-precision rational numbers.
     6	
     7	package big
     8	
     9	import (
    10		"fmt"
    11		"math"
    12	)
    13	
    14	// A Rat represents a quotient a/b of arbitrary precision.
    15	// The zero value for a Rat represents the value 0.
    16	//
    17	// Operations always take pointer arguments (*Rat) rather
    18	// than Rat values, and each unique Rat value requires
    19	// its own unique *Rat pointer. To "copy" a Rat value,
    20	// an existing (or newly allocated) Rat must be set to
    21	// a new value using the Rat.Set method; shallow copies
    22	// of Rats are not supported and may lead to errors.
    23	type Rat struct {
    24		// To make zero values for Rat work w/o initialization,
    25		// a zero value of b (len(b) == 0) acts like b == 1.
    26		// a.neg determines the sign of the Rat, b.neg is ignored.
    27		a, b Int
    28	}
    29	
    30	// NewRat creates a new Rat with numerator a and denominator b.
    31	func NewRat(a, b int64) *Rat {
    32		return new(Rat).SetFrac64(a, b)
    33	}
    34	
    35	// SetFloat64 sets z to exactly f and returns z.
    36	// If f is not finite, SetFloat returns nil.
    37	func (z *Rat) SetFloat64(f float64) *Rat {
    38		const expMask = 1<<11 - 1
    39		bits := math.Float64bits(f)
    40		mantissa := bits & (1<<52 - 1)
    41		exp := int((bits >> 52) & expMask)
    42		switch exp {
    43		case expMask: // non-finite
    44			return nil
    45		case 0: // denormal
    46			exp -= 1022
    47		default: // normal
    48			mantissa |= 1 << 52
    49			exp -= 1023
    50		}
    51	
    52		shift := 52 - exp
    53	
    54		// Optimization (?): partially pre-normalise.
    55		for mantissa&1 == 0 && shift > 0 {
    56			mantissa >>= 1
    57			shift--
    58		}
    59	
    60		z.a.SetUint64(mantissa)
    61		z.a.neg = f < 0
    62		z.b.Set(intOne)
    63		if shift > 0 {
    64			z.b.Lsh(&z.b, uint(shift))
    65		} else {
    66			z.a.Lsh(&z.a, uint(-shift))
    67		}
    68		return z.norm()
    69	}
    70	
    71	// quotToFloat32 returns the non-negative float32 value
    72	// nearest to the quotient a/b, using round-to-even in
    73	// halfway cases. It does not mutate its arguments.
    74	// Preconditions: b is non-zero; a and b have no common factors.
    75	func quotToFloat32(a, b nat) (f float32, exact bool) {
    76		const (
    77			// float size in bits
    78			Fsize = 32
    79	
    80			// mantissa
    81			Msize  = 23
    82			Msize1 = Msize + 1 // incl. implicit 1
    83			Msize2 = Msize1 + 1
    84	
    85			// exponent
    86			Esize = Fsize - Msize1
    87			Ebias = 1<<(Esize-1) - 1
    88			Emin  = 1 - Ebias
    89			Emax  = Ebias
    90		)
    91	
    92		// TODO(adonovan): specialize common degenerate cases: 1.0, integers.
    93		alen := a.bitLen()
    94		if alen == 0 {
    95			return 0, true
    96		}
    97		blen := b.bitLen()
    98		if blen == 0 {
    99			panic("division by zero")
   100		}
   101	
   102		// 1. Left-shift A or B such that quotient A/B is in [1<<Msize1, 1<<(Msize2+1)
   103		// (Msize2 bits if A < B when they are left-aligned, Msize2+1 bits if A >= B).
   104		// This is 2 or 3 more than the float32 mantissa field width of Msize:
   105		// - the optional extra bit is shifted away in step 3 below.
   106		// - the high-order 1 is omitted in "normal" representation;
   107		// - the low-order 1 will be used during rounding then discarded.
   108		exp := alen - blen
   109		var a2, b2 nat
   110		a2 = a2.set(a)
   111		b2 = b2.set(b)
   112		if shift := Msize2 - exp; shift > 0 {
   113			a2 = a2.shl(a2, uint(shift))
   114		} else if shift < 0 {
   115			b2 = b2.shl(b2, uint(-shift))
   116		}
   117	
   118		// 2. Compute quotient and remainder (q, r).  NB: due to the
   119		// extra shift, the low-order bit of q is logically the
   120		// high-order bit of r.
   121		var q nat
   122		q, r := q.div(a2, a2, b2) // (recycle a2)
   123		mantissa := low32(q)
   124		haveRem := len(r) > 0 // mantissa&1 && !haveRem => remainder is exactly half
   125	
   126		// 3. If quotient didn't fit in Msize2 bits, redo division by b2<<1
   127		// (in effect---we accomplish this incrementally).
   128		if mantissa>>Msize2 == 1 {
   129			if mantissa&1 == 1 {
   130				haveRem = true
   131			}
   132			mantissa >>= 1
   133			exp++
   134		}
   135		if mantissa>>Msize1 != 1 {
   136			panic(fmt.Sprintf("expected exactly %d bits of result", Msize2))
   137		}
   138	
   139		// 4. Rounding.
   140		if Emin-Msize <= exp && exp <= Emin {
   141			// Denormal case; lose 'shift' bits of precision.
   142			shift := uint(Emin - (exp - 1)) // [1..Esize1)
   143			lostbits := mantissa & (1<<shift - 1)
   144			haveRem = haveRem || lostbits != 0
   145			mantissa >>= shift
   146			exp = 2 - Ebias // == exp + shift
   147		}
   148		// Round q using round-half-to-even.
   149		exact = !haveRem
   150		if mantissa&1 != 0 {
   151			exact = false
   152			if haveRem || mantissa&2 != 0 {
   153				if mantissa++; mantissa >= 1<<Msize2 {
   154					// Complete rollover 11...1 => 100...0, so shift is safe
   155					mantissa >>= 1
   156					exp++
   157				}
   158			}
   159		}
   160		mantissa >>= 1 // discard rounding bit.  Mantissa now scaled by 1<<Msize1.
   161	
   162		f = float32(math.Ldexp(float64(mantissa), exp-Msize1))
   163		if math.IsInf(float64(f), 0) {
   164			exact = false
   165		}
   166		return
   167	}
   168	
   169	// quotToFloat64 returns the non-negative float64 value
   170	// nearest to the quotient a/b, using round-to-even in
   171	// halfway cases. It does not mutate its arguments.
   172	// Preconditions: b is non-zero; a and b have no common factors.
   173	func quotToFloat64(a, b nat) (f float64, exact bool) {
   174		const (
   175			// float size in bits
   176			Fsize = 64
   177	
   178			// mantissa
   179			Msize  = 52
   180			Msize1 = Msize + 1 // incl. implicit 1
   181			Msize2 = Msize1 + 1
   182	
   183			// exponent
   184			Esize = Fsize - Msize1
   185			Ebias = 1<<(Esize-1) - 1
   186			Emin  = 1 - Ebias
   187			Emax  = Ebias
   188		)
   189	
   190		// TODO(adonovan): specialize common degenerate cases: 1.0, integers.
   191		alen := a.bitLen()
   192		if alen == 0 {
   193			return 0, true
   194		}
   195		blen := b.bitLen()
   196		if blen == 0 {
   197			panic("division by zero")
   198		}
   199	
   200		// 1. Left-shift A or B such that quotient A/B is in [1<<Msize1, 1<<(Msize2+1)
   201		// (Msize2 bits if A < B when they are left-aligned, Msize2+1 bits if A >= B).
   202		// This is 2 or 3 more than the float64 mantissa field width of Msize:
   203		// - the optional extra bit is shifted away in step 3 below.
   204		// - the high-order 1 is omitted in "normal" representation;
   205		// - the low-order 1 will be used during rounding then discarded.
   206		exp := alen - blen
   207		var a2, b2 nat
   208		a2 = a2.set(a)
   209		b2 = b2.set(b)
   210		if shift := Msize2 - exp; shift > 0 {
   211			a2 = a2.shl(a2, uint(shift))
   212		} else if shift < 0 {
   213			b2 = b2.shl(b2, uint(-shift))
   214		}
   215	
   216		// 2. Compute quotient and remainder (q, r).  NB: due to the
   217		// extra shift, the low-order bit of q is logically the
   218		// high-order bit of r.
   219		var q nat
   220		q, r := q.div(a2, a2, b2) // (recycle a2)
   221		mantissa := low64(q)
   222		haveRem := len(r) > 0 // mantissa&1 && !haveRem => remainder is exactly half
   223	
   224		// 3. If quotient didn't fit in Msize2 bits, redo division by b2<<1
   225		// (in effect---we accomplish this incrementally).
   226		if mantissa>>Msize2 == 1 {
   227			if mantissa&1 == 1 {
   228				haveRem = true
   229			}
   230			mantissa >>= 1
   231			exp++
   232		}
   233		if mantissa>>Msize1 != 1 {
   234			panic(fmt.Sprintf("expected exactly %d bits of result", Msize2))
   235		}
   236	
   237		// 4. Rounding.
   238		if Emin-Msize <= exp && exp <= Emin {
   239			// Denormal case; lose 'shift' bits of precision.
   240			shift := uint(Emin - (exp - 1)) // [1..Esize1)
   241			lostbits := mantissa & (1<<shift - 1)
   242			haveRem = haveRem || lostbits != 0
   243			mantissa >>= shift
   244			exp = 2 - Ebias // == exp + shift
   245		}
   246		// Round q using round-half-to-even.
   247		exact = !haveRem
   248		if mantissa&1 != 0 {
   249			exact = false
   250			if haveRem || mantissa&2 != 0 {
   251				if mantissa++; mantissa >= 1<<Msize2 {
   252					// Complete rollover 11...1 => 100...0, so shift is safe
   253					mantissa >>= 1
   254					exp++
   255				}
   256			}
   257		}
   258		mantissa >>= 1 // discard rounding bit.  Mantissa now scaled by 1<<Msize1.
   259	
   260		f = math.Ldexp(float64(mantissa), exp-Msize1)
   261		if math.IsInf(f, 0) {
   262			exact = false
   263		}
   264		return
   265	}
   266	
   267	// Float32 returns the nearest float32 value for x and a bool indicating
   268	// whether f represents x exactly. If the magnitude of x is too large to
   269	// be represented by a float32, f is an infinity and exact is false.
   270	// The sign of f always matches the sign of x, even if f == 0.
   271	func (x *Rat) Float32() (f float32, exact bool) {
   272		b := x.b.abs
   273		if len(b) == 0 {
   274			b = b.set(natOne) // materialize denominator
   275		}
   276		f, exact = quotToFloat32(x.a.abs, b)
   277		if x.a.neg {
   278			f = -f
   279		}
   280		return
   281	}
   282	
   283	// Float64 returns the nearest float64 value for x and a bool indicating
   284	// whether f represents x exactly. If the magnitude of x is too large to
   285	// be represented by a float64, f is an infinity and exact is false.
   286	// The sign of f always matches the sign of x, even if f == 0.
   287	func (x *Rat) Float64() (f float64, exact bool) {
   288		b := x.b.abs
   289		if len(b) == 0 {
   290			b = b.set(natOne) // materialize denominator
   291		}
   292		f, exact = quotToFloat64(x.a.abs, b)
   293		if x.a.neg {
   294			f = -f
   295		}
   296		return
   297	}
   298	
   299	// SetFrac sets z to a/b and returns z.
   300	func (z *Rat) SetFrac(a, b *Int) *Rat {
   301		z.a.neg = a.neg != b.neg
   302		babs := b.abs
   303		if len(babs) == 0 {
   304			panic("division by zero")
   305		}
   306		if &z.a == b || alias(z.a.abs, babs) {
   307			babs = nat(nil).set(babs) // make a copy
   308		}
   309		z.a.abs = z.a.abs.set(a.abs)
   310		z.b.abs = z.b.abs.set(babs)
   311		return z.norm()
   312	}
   313	
   314	// SetFrac64 sets z to a/b and returns z.
   315	func (z *Rat) SetFrac64(a, b int64) *Rat {
   316		z.a.SetInt64(a)
   317		if b == 0 {
   318			panic("division by zero")
   319		}
   320		if b < 0 {
   321			b = -b
   322			z.a.neg = !z.a.neg
   323		}
   324		z.b.abs = z.b.abs.setUint64(uint64(b))
   325		return z.norm()
   326	}
   327	
   328	// SetInt sets z to x (by making a copy of x) and returns z.
   329	func (z *Rat) SetInt(x *Int) *Rat {
   330		z.a.Set(x)
   331		z.b.abs = z.b.abs[:0]
   332		return z
   333	}
   334	
   335	// SetInt64 sets z to x and returns z.
   336	func (z *Rat) SetInt64(x int64) *Rat {
   337		z.a.SetInt64(x)
   338		z.b.abs = z.b.abs[:0]
   339		return z
   340	}
   341	
   342	// SetUint64 sets z to x and returns z.
   343	func (z *Rat) SetUint64(x uint64) *Rat {
   344		z.a.SetUint64(x)
   345		z.b.abs = z.b.abs[:0]
   346		return z
   347	}
   348	
   349	// Set sets z to x (by making a copy of x) and returns z.
   350	func (z *Rat) Set(x *Rat) *Rat {
   351		if z != x {
   352			z.a.Set(&x.a)
   353			z.b.Set(&x.b)
   354		}
   355		return z
   356	}
   357	
   358	// Abs sets z to |x| (the absolute value of x) and returns z.
   359	func (z *Rat) Abs(x *Rat) *Rat {
   360		z.Set(x)
   361		z.a.neg = false
   362		return z
   363	}
   364	
   365	// Neg sets z to -x and returns z.
   366	func (z *Rat) Neg(x *Rat) *Rat {
   367		z.Set(x)
   368		z.a.neg = len(z.a.abs) > 0 && !z.a.neg // 0 has no sign
   369		return z
   370	}
   371	
   372	// Inv sets z to 1/x and returns z.
   373	func (z *Rat) Inv(x *Rat) *Rat {
   374		if len(x.a.abs) == 0 {
   375			panic("division by zero")
   376		}
   377		z.Set(x)
   378		a := z.b.abs
   379		if len(a) == 0 {
   380			a = a.set(natOne) // materialize numerator
   381		}
   382		b := z.a.abs
   383		if b.cmp(natOne) == 0 {
   384			b = b[:0] // normalize denominator
   385		}
   386		z.a.abs, z.b.abs = a, b // sign doesn't change
   387		return z
   388	}
   389	
   390	// Sign returns:
   391	//
   392	//	-1 if x <  0
   393	//	 0 if x == 0
   394	//	+1 if x >  0
   395	//
   396	func (x *Rat) Sign() int {
   397		return x.a.Sign()
   398	}
   399	
   400	// IsInt reports whether the denominator of x is 1.
   401	func (x *Rat) IsInt() bool {
   402		return len(x.b.abs) == 0 || x.b.abs.cmp(natOne) == 0
   403	}
   404	
   405	// Num returns the numerator of x; it may be <= 0.
   406	// The result is a reference to x's numerator; it
   407	// may change if a new value is assigned to x, and vice versa.
   408	// The sign of the numerator corresponds to the sign of x.
   409	func (x *Rat) Num() *Int {
   410		return &x.a
   411	}
   412	
   413	// Denom returns the denominator of x; it is always > 0.
   414	// The result is a reference to x's denominator; it
   415	// may change if a new value is assigned to x, and vice versa.
   416	func (x *Rat) Denom() *Int {
   417		x.b.neg = false // the result is always >= 0
   418		if len(x.b.abs) == 0 {
   419			x.b.abs = x.b.abs.set(natOne) // materialize denominator
   420		}
   421		return &x.b
   422	}
   423	
   424	func (z *Rat) norm() *Rat {
   425		switch {
   426		case len(z.a.abs) == 0:
   427			// z == 0 - normalize sign and denominator
   428			z.a.neg = false
   429			z.b.abs = z.b.abs[:0]
   430		case len(z.b.abs) == 0:
   431			// z is normalized int - nothing to do
   432		case z.b.abs.cmp(natOne) == 0:
   433			// z is int - normalize denominator
   434			z.b.abs = z.b.abs[:0]
   435		default:
   436			neg := z.a.neg
   437			z.a.neg = false
   438			z.b.neg = false
   439			if f := NewInt(0).lehmerGCD(nil, nil, &z.a, &z.b); f.Cmp(intOne) != 0 {
   440				z.a.abs, _ = z.a.abs.div(nil, z.a.abs, f.abs)
   441				z.b.abs, _ = z.b.abs.div(nil, z.b.abs, f.abs)
   442				if z.b.abs.cmp(natOne) == 0 {
   443					// z is int - normalize denominator
   444					z.b.abs = z.b.abs[:0]
   445				}
   446			}
   447			z.a.neg = neg
   448		}
   449		return z
   450	}
   451	
   452	// mulDenom sets z to the denominator product x*y (by taking into
   453	// account that 0 values for x or y must be interpreted as 1) and
   454	// returns z.
   455	func mulDenom(z, x, y nat) nat {
   456		switch {
   457		case len(x) == 0:
   458			return z.set(y)
   459		case len(y) == 0:
   460			return z.set(x)
   461		}
   462		return z.mul(x, y)
   463	}
   464	
   465	// scaleDenom sets z to the product x*f.
   466	// If f == 0 (zero value of denominator), z is set to (a copy of) x.
   467	func (z *Int) scaleDenom(x *Int, f nat) {
   468		if len(f) == 0 {
   469			z.Set(x)
   470			return
   471		}
   472		z.abs = z.abs.mul(x.abs, f)
   473		z.neg = x.neg
   474	}
   475	
   476	// Cmp compares x and y and returns:
   477	//
   478	//   -1 if x <  y
   479	//    0 if x == y
   480	//   +1 if x >  y
   481	//
   482	func (x *Rat) Cmp(y *Rat) int {
   483		var a, b Int
   484		a.scaleDenom(&x.a, y.b.abs)
   485		b.scaleDenom(&y.a, x.b.abs)
   486		return a.Cmp(&b)
   487	}
   488	
   489	// Add sets z to the sum x+y and returns z.
   490	func (z *Rat) Add(x, y *Rat) *Rat {
   491		var a1, a2 Int
   492		a1.scaleDenom(&x.a, y.b.abs)
   493		a2.scaleDenom(&y.a, x.b.abs)
   494		z.a.Add(&a1, &a2)
   495		z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
   496		return z.norm()
   497	}
   498	
   499	// Sub sets z to the difference x-y and returns z.
   500	func (z *Rat) Sub(x, y *Rat) *Rat {
   501		var a1, a2 Int
   502		a1.scaleDenom(&x.a, y.b.abs)
   503		a2.scaleDenom(&y.a, x.b.abs)
   504		z.a.Sub(&a1, &a2)
   505		z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
   506		return z.norm()
   507	}
   508	
   509	// Mul sets z to the product x*y and returns z.
   510	func (z *Rat) Mul(x, y *Rat) *Rat {
   511		if x == y {
   512			// a squared Rat is positive and can't be reduced
   513			z.a.neg = false
   514			z.a.abs = z.a.abs.sqr(x.a.abs)
   515			z.b.abs = z.b.abs.sqr(x.b.abs)
   516			return z
   517		}
   518		z.a.Mul(&x.a, &y.a)
   519		z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
   520		return z.norm()
   521	}
   522	
   523	// Quo sets z to the quotient x/y and returns z.
   524	// If y == 0, a division-by-zero run-time panic occurs.
   525	func (z *Rat) Quo(x, y *Rat) *Rat {
   526		if len(y.a.abs) == 0 {
   527			panic("division by zero")
   528		}
   529		var a, b Int
   530		a.scaleDenom(&x.a, y.b.abs)
   531		b.scaleDenom(&y.a, x.b.abs)
   532		z.a.abs = a.abs
   533		z.b.abs = b.abs
   534		z.a.neg = a.neg != b.neg
   535		return z.norm()
   536	}
   537	

View as plain text