...
Run Format

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		"encoding/binary"
    11		"errors"
    12		"fmt"
    13		"math"
    14		"strings"
    15	)
    16	
    17	// A Rat represents a quotient a/b of arbitrary precision.
    18	// The zero value for a Rat represents the value 0.
    19	type Rat struct {
    20		// To make zero values for Rat work w/o initialization,
    21		// a zero value of b (len(b) == 0) acts like b == 1.
    22		// a.neg determines the sign of the Rat, b.neg is ignored.
    23		a, b Int
    24	}
    25	
    26	// NewRat creates a new Rat with numerator a and denominator b.
    27	func NewRat(a, b int64) *Rat {
    28		return new(Rat).SetFrac64(a, b)
    29	}
    30	
    31	// SetFloat64 sets z to exactly f and returns z.
    32	// If f is not finite, SetFloat returns nil.
    33	func (z *Rat) SetFloat64(f float64) *Rat {
    34		const expMask = 1<<11 - 1
    35		bits := math.Float64bits(f)
    36		mantissa := bits & (1<<52 - 1)
    37		exp := int((bits >> 52) & expMask)
    38		switch exp {
    39		case expMask: // non-finite
    40			return nil
    41		case 0: // denormal
    42			exp -= 1022
    43		default: // normal
    44			mantissa |= 1 << 52
    45			exp -= 1023
    46		}
    47	
    48		shift := 52 - exp
    49	
    50		// Optimization (?): partially pre-normalise.
    51		for mantissa&1 == 0 && shift > 0 {
    52			mantissa >>= 1
    53			shift--
    54		}
    55	
    56		z.a.SetUint64(mantissa)
    57		z.a.neg = f < 0
    58		z.b.Set(intOne)
    59		if shift > 0 {
    60			z.b.Lsh(&z.b, uint(shift))
    61		} else {
    62			z.a.Lsh(&z.a, uint(-shift))
    63		}
    64		return z.norm()
    65	}
    66	
    67	// isFinite reports whether f represents a finite rational value.
    68	// It is equivalent to !math.IsNan(f) && !math.IsInf(f, 0).
    69	func isFinite(f float64) bool {
    70		return math.Abs(f) <= math.MaxFloat64
    71	}
    72	
    73	// low64 returns the least significant 64 bits of natural number z.
    74	func low64(z nat) uint64 {
    75		if len(z) == 0 {
    76			return 0
    77		}
    78		if _W == 32 && len(z) > 1 {
    79			return uint64(z[1])<<32 | uint64(z[0])
    80		}
    81		return uint64(z[0])
    82	}
    83	
    84	// quotToFloat returns the non-negative IEEE 754 double-precision
    85	// value nearest to the quotient a/b, using round-to-even in halfway
    86	// cases.  It does not mutate its arguments.
    87	// Preconditions: b is non-zero; a and b have no common factors.
    88	func quotToFloat(a, b nat) (f float64, exact bool) {
    89		// TODO(adonovan): specialize common degenerate cases: 1.0, integers.
    90		alen := a.bitLen()
    91		if alen == 0 {
    92			return 0, true
    93		}
    94		blen := b.bitLen()
    95		if blen == 0 {
    96			panic("division by zero")
    97		}
    98	
    99		// 1. Left-shift A or B such that quotient A/B is in [1<<53, 1<<55).
   100		// (54 bits if A<B when they are left-aligned, 55 bits if A>=B.)
   101		// This is 2 or 3 more than the float64 mantissa field width of 52:
   102		// - the optional extra bit is shifted away in step 3 below.
   103		// - the high-order 1 is omitted in float64 "normal" representation;
   104		// - the low-order 1 will be used during rounding then discarded.
   105		exp := alen - blen
   106		var a2, b2 nat
   107		a2 = a2.set(a)
   108		b2 = b2.set(b)
   109		if shift := 54 - exp; shift > 0 {
   110			a2 = a2.shl(a2, uint(shift))
   111		} else if shift < 0 {
   112			b2 = b2.shl(b2, uint(-shift))
   113		}
   114	
   115		// 2. Compute quotient and remainder (q, r).  NB: due to the
   116		// extra shift, the low-order bit of q is logically the
   117		// high-order bit of r.
   118		var q nat
   119		q, r := q.div(a2, a2, b2) // (recycle a2)
   120		mantissa := low64(q)
   121		haveRem := len(r) > 0 // mantissa&1 && !haveRem => remainder is exactly half
   122	
   123		// 3. If quotient didn't fit in 54 bits, re-do division by b2<<1
   124		// (in effect---we accomplish this incrementally).
   125		if mantissa>>54 == 1 {
   126			if mantissa&1 == 1 {
   127				haveRem = true
   128			}
   129			mantissa >>= 1
   130			exp++
   131		}
   132		if mantissa>>53 != 1 {
   133			panic("expected exactly 54 bits of result")
   134		}
   135	
   136		// 4. Rounding.
   137		if -1022-52 <= exp && exp <= -1022 {
   138			// Denormal case; lose 'shift' bits of precision.
   139			shift := uint64(-1022 - (exp - 1)) // [1..53)
   140			lostbits := mantissa & (1<<shift - 1)
   141			haveRem = haveRem || lostbits != 0
   142			mantissa >>= shift
   143			exp = -1023 + 2
   144		}
   145		// Round q using round-half-to-even.
   146		exact = !haveRem
   147		if mantissa&1 != 0 {
   148			exact = false
   149			if haveRem || mantissa&2 != 0 {
   150				if mantissa++; mantissa >= 1<<54 {
   151					// Complete rollover 11...1 => 100...0, so shift is safe
   152					mantissa >>= 1
   153					exp++
   154				}
   155			}
   156		}
   157		mantissa >>= 1 // discard rounding bit.  Mantissa now scaled by 2^53.
   158	
   159		f = math.Ldexp(float64(mantissa), exp-53)
   160		if math.IsInf(f, 0) {
   161			exact = false
   162		}
   163		return
   164	}
   165	
   166	// Float64 returns the nearest float64 value for x and a bool indicating
   167	// whether f represents x exactly. If the magnitude of x is too large to
   168	// be represented by a float64, f is an infinity and exact is false.
   169	// The sign of f always matches the sign of x, even if f == 0.
   170	func (x *Rat) Float64() (f float64, exact bool) {
   171		b := x.b.abs
   172		if len(b) == 0 {
   173			b = b.set(natOne) // materialize denominator
   174		}
   175		f, exact = quotToFloat(x.a.abs, b)
   176		if x.a.neg {
   177			f = -f
   178		}
   179		return
   180	}
   181	
   182	// SetFrac sets z to a/b and returns z.
   183	func (z *Rat) SetFrac(a, b *Int) *Rat {
   184		z.a.neg = a.neg != b.neg
   185		babs := b.abs
   186		if len(babs) == 0 {
   187			panic("division by zero")
   188		}
   189		if &z.a == b || alias(z.a.abs, babs) {
   190			babs = nat(nil).set(babs) // make a copy
   191		}
   192		z.a.abs = z.a.abs.set(a.abs)
   193		z.b.abs = z.b.abs.set(babs)
   194		return z.norm()
   195	}
   196	
   197	// SetFrac64 sets z to a/b and returns z.
   198	func (z *Rat) SetFrac64(a, b int64) *Rat {
   199		z.a.SetInt64(a)
   200		if b == 0 {
   201			panic("division by zero")
   202		}
   203		if b < 0 {
   204			b = -b
   205			z.a.neg = !z.a.neg
   206		}
   207		z.b.abs = z.b.abs.setUint64(uint64(b))
   208		return z.norm()
   209	}
   210	
   211	// SetInt sets z to x (by making a copy of x) and returns z.
   212	func (z *Rat) SetInt(x *Int) *Rat {
   213		z.a.Set(x)
   214		z.b.abs = z.b.abs.make(0)
   215		return z
   216	}
   217	
   218	// SetInt64 sets z to x and returns z.
   219	func (z *Rat) SetInt64(x int64) *Rat {
   220		z.a.SetInt64(x)
   221		z.b.abs = z.b.abs.make(0)
   222		return z
   223	}
   224	
   225	// Set sets z to x (by making a copy of x) and returns z.
   226	func (z *Rat) Set(x *Rat) *Rat {
   227		if z != x {
   228			z.a.Set(&x.a)
   229			z.b.Set(&x.b)
   230		}
   231		return z
   232	}
   233	
   234	// Abs sets z to |x| (the absolute value of x) and returns z.
   235	func (z *Rat) Abs(x *Rat) *Rat {
   236		z.Set(x)
   237		z.a.neg = false
   238		return z
   239	}
   240	
   241	// Neg sets z to -x and returns z.
   242	func (z *Rat) Neg(x *Rat) *Rat {
   243		z.Set(x)
   244		z.a.neg = len(z.a.abs) > 0 && !z.a.neg // 0 has no sign
   245		return z
   246	}
   247	
   248	// Inv sets z to 1/x and returns z.
   249	func (z *Rat) Inv(x *Rat) *Rat {
   250		if len(x.a.abs) == 0 {
   251			panic("division by zero")
   252		}
   253		z.Set(x)
   254		a := z.b.abs
   255		if len(a) == 0 {
   256			a = a.set(natOne) // materialize numerator
   257		}
   258		b := z.a.abs
   259		if b.cmp(natOne) == 0 {
   260			b = b.make(0) // normalize denominator
   261		}
   262		z.a.abs, z.b.abs = a, b // sign doesn't change
   263		return z
   264	}
   265	
   266	// Sign returns:
   267	//
   268	//	-1 if x <  0
   269	//	 0 if x == 0
   270	//	+1 if x >  0
   271	//
   272	func (x *Rat) Sign() int {
   273		return x.a.Sign()
   274	}
   275	
   276	// IsInt returns true if the denominator of x is 1.
   277	func (x *Rat) IsInt() bool {
   278		return len(x.b.abs) == 0 || x.b.abs.cmp(natOne) == 0
   279	}
   280	
   281	// Num returns the numerator of x; it may be <= 0.
   282	// The result is a reference to x's numerator; it
   283	// may change if a new value is assigned to x, and vice versa.
   284	// The sign of the numerator corresponds to the sign of x.
   285	func (x *Rat) Num() *Int {
   286		return &x.a
   287	}
   288	
   289	// Denom returns the denominator of x; it is always > 0.
   290	// The result is a reference to x's denominator; it
   291	// may change if a new value is assigned to x, and vice versa.
   292	func (x *Rat) Denom() *Int {
   293		x.b.neg = false // the result is always >= 0
   294		if len(x.b.abs) == 0 {
   295			x.b.abs = x.b.abs.set(natOne) // materialize denominator
   296		}
   297		return &x.b
   298	}
   299	
   300	func (z *Rat) norm() *Rat {
   301		switch {
   302		case len(z.a.abs) == 0:
   303			// z == 0 - normalize sign and denominator
   304			z.a.neg = false
   305			z.b.abs = z.b.abs.make(0)
   306		case len(z.b.abs) == 0:
   307			// z is normalized int - nothing to do
   308		case z.b.abs.cmp(natOne) == 0:
   309			// z is int - normalize denominator
   310			z.b.abs = z.b.abs.make(0)
   311		default:
   312			neg := z.a.neg
   313			z.a.neg = false
   314			z.b.neg = false
   315			if f := NewInt(0).binaryGCD(&z.a, &z.b); f.Cmp(intOne) != 0 {
   316				z.a.abs, _ = z.a.abs.div(nil, z.a.abs, f.abs)
   317				z.b.abs, _ = z.b.abs.div(nil, z.b.abs, f.abs)
   318				if z.b.abs.cmp(natOne) == 0 {
   319					// z is int - normalize denominator
   320					z.b.abs = z.b.abs.make(0)
   321				}
   322			}
   323			z.a.neg = neg
   324		}
   325		return z
   326	}
   327	
   328	// mulDenom sets z to the denominator product x*y (by taking into
   329	// account that 0 values for x or y must be interpreted as 1) and
   330	// returns z.
   331	func mulDenom(z, x, y nat) nat {
   332		switch {
   333		case len(x) == 0:
   334			return z.set(y)
   335		case len(y) == 0:
   336			return z.set(x)
   337		}
   338		return z.mul(x, y)
   339	}
   340	
   341	// scaleDenom computes x*f.
   342	// If f == 0 (zero value of denominator), the result is (a copy of) x.
   343	func scaleDenom(x *Int, f nat) *Int {
   344		var z Int
   345		if len(f) == 0 {
   346			return z.Set(x)
   347		}
   348		z.abs = z.abs.mul(x.abs, f)
   349		z.neg = x.neg
   350		return &z
   351	}
   352	
   353	// Cmp compares x and y and returns:
   354	//
   355	//   -1 if x <  y
   356	//    0 if x == y
   357	//   +1 if x >  y
   358	//
   359	func (x *Rat) Cmp(y *Rat) int {
   360		return scaleDenom(&x.a, y.b.abs).Cmp(scaleDenom(&y.a, x.b.abs))
   361	}
   362	
   363	// Add sets z to the sum x+y and returns z.
   364	func (z *Rat) Add(x, y *Rat) *Rat {
   365		a1 := scaleDenom(&x.a, y.b.abs)
   366		a2 := scaleDenom(&y.a, x.b.abs)
   367		z.a.Add(a1, a2)
   368		z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
   369		return z.norm()
   370	}
   371	
   372	// Sub sets z to the difference x-y and returns z.
   373	func (z *Rat) Sub(x, y *Rat) *Rat {
   374		a1 := scaleDenom(&x.a, y.b.abs)
   375		a2 := scaleDenom(&y.a, x.b.abs)
   376		z.a.Sub(a1, a2)
   377		z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
   378		return z.norm()
   379	}
   380	
   381	// Mul sets z to the product x*y and returns z.
   382	func (z *Rat) Mul(x, y *Rat) *Rat {
   383		z.a.Mul(&x.a, &y.a)
   384		z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
   385		return z.norm()
   386	}
   387	
   388	// Quo sets z to the quotient x/y and returns z.
   389	// If y == 0, a division-by-zero run-time panic occurs.
   390	func (z *Rat) Quo(x, y *Rat) *Rat {
   391		if len(y.a.abs) == 0 {
   392			panic("division by zero")
   393		}
   394		a := scaleDenom(&x.a, y.b.abs)
   395		b := scaleDenom(&y.a, x.b.abs)
   396		z.a.abs = a.abs
   397		z.b.abs = b.abs
   398		z.a.neg = a.neg != b.neg
   399		return z.norm()
   400	}
   401	
   402	func ratTok(ch rune) bool {
   403		return strings.IndexRune("+-/0123456789.eE", ch) >= 0
   404	}
   405	
   406	// Scan is a support routine for fmt.Scanner. It accepts the formats
   407	// 'e', 'E', 'f', 'F', 'g', 'G', and 'v'. All formats are equivalent.
   408	func (z *Rat) Scan(s fmt.ScanState, ch rune) error {
   409		tok, err := s.Token(true, ratTok)
   410		if err != nil {
   411			return err
   412		}
   413		if strings.IndexRune("efgEFGv", ch) < 0 {
   414			return errors.New("Rat.Scan: invalid verb")
   415		}
   416		if _, ok := z.SetString(string(tok)); !ok {
   417			return errors.New("Rat.Scan: invalid syntax")
   418		}
   419		return nil
   420	}
   421	
   422	// SetString sets z to the value of s and returns z and a boolean indicating
   423	// success. s can be given as a fraction "a/b" or as a floating-point number
   424	// optionally followed by an exponent. If the operation failed, the value of
   425	// z is undefined but the returned value is nil.
   426	func (z *Rat) SetString(s string) (*Rat, bool) {
   427		if len(s) == 0 {
   428			return nil, false
   429		}
   430	
   431		// check for a quotient
   432		sep := strings.Index(s, "/")
   433		if sep >= 0 {
   434			if _, ok := z.a.SetString(s[0:sep], 10); !ok {
   435				return nil, false
   436			}
   437			s = s[sep+1:]
   438			var err error
   439			if z.b.abs, _, err = z.b.abs.scan(strings.NewReader(s), 10); err != nil {
   440				return nil, false
   441			}
   442			return z.norm(), true
   443		}
   444	
   445		// check for a decimal point
   446		sep = strings.Index(s, ".")
   447		// check for an exponent
   448		e := strings.IndexAny(s, "eE")
   449		var exp Int
   450		if e >= 0 {
   451			if e < sep {
   452				// The E must come after the decimal point.
   453				return nil, false
   454			}
   455			if _, ok := exp.SetString(s[e+1:], 10); !ok {
   456				return nil, false
   457			}
   458			s = s[0:e]
   459		}
   460		if sep >= 0 {
   461			s = s[0:sep] + s[sep+1:]
   462			exp.Sub(&exp, NewInt(int64(len(s)-sep)))
   463		}
   464	
   465		if _, ok := z.a.SetString(s, 10); !ok {
   466			return nil, false
   467		}
   468		powTen := nat(nil).expNN(natTen, exp.abs, nil)
   469		if exp.neg {
   470			z.b.abs = powTen
   471			z.norm()
   472		} else {
   473			z.a.abs = z.a.abs.mul(z.a.abs, powTen)
   474			z.b.abs = z.b.abs.make(0)
   475		}
   476	
   477		return z, true
   478	}
   479	
   480	// String returns a string representation of x in the form "a/b" (even if b == 1).
   481	func (x *Rat) String() string {
   482		s := "/1"
   483		if len(x.b.abs) != 0 {
   484			s = "/" + x.b.abs.decimalString()
   485		}
   486		return x.a.String() + s
   487	}
   488	
   489	// RatString returns a string representation of x in the form "a/b" if b != 1,
   490	// and in the form "a" if b == 1.
   491	func (x *Rat) RatString() string {
   492		if x.IsInt() {
   493			return x.a.String()
   494		}
   495		return x.String()
   496	}
   497	
   498	// FloatString returns a string representation of x in decimal form with prec
   499	// digits of precision after the decimal point and the last digit rounded.
   500	func (x *Rat) FloatString(prec int) string {
   501		if x.IsInt() {
   502			s := x.a.String()
   503			if prec > 0 {
   504				s += "." + strings.Repeat("0", prec)
   505			}
   506			return s
   507		}
   508		// x.b.abs != 0
   509	
   510		q, r := nat(nil).div(nat(nil), x.a.abs, x.b.abs)
   511	
   512		p := natOne
   513		if prec > 0 {
   514			p = nat(nil).expNN(natTen, nat(nil).setUint64(uint64(prec)), nil)
   515		}
   516	
   517		r = r.mul(r, p)
   518		r, r2 := r.div(nat(nil), r, x.b.abs)
   519	
   520		// see if we need to round up
   521		r2 = r2.add(r2, r2)
   522		if x.b.abs.cmp(r2) <= 0 {
   523			r = r.add(r, natOne)
   524			if r.cmp(p) >= 0 {
   525				q = nat(nil).add(q, natOne)
   526				r = nat(nil).sub(r, p)
   527			}
   528		}
   529	
   530		s := q.decimalString()
   531		if x.a.neg {
   532			s = "-" + s
   533		}
   534	
   535		if prec > 0 {
   536			rs := r.decimalString()
   537			leadingZeros := prec - len(rs)
   538			s += "." + strings.Repeat("0", leadingZeros) + rs
   539		}
   540	
   541		return s
   542	}
   543	
   544	// Gob codec version. Permits backward-compatible changes to the encoding.
   545	const ratGobVersion byte = 1
   546	
   547	// GobEncode implements the gob.GobEncoder interface.
   548	func (x *Rat) GobEncode() ([]byte, error) {
   549		if x == nil {
   550			return nil, nil
   551		}
   552		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)
   553		i := x.b.abs.bytes(buf)
   554		j := x.a.abs.bytes(buf[0:i])
   555		n := i - j
   556		if int(uint32(n)) != n {
   557			// this should never happen
   558			return nil, errors.New("Rat.GobEncode: numerator too large")
   559		}
   560		binary.BigEndian.PutUint32(buf[j-4:j], uint32(n))
   561		j -= 1 + 4
   562		b := ratGobVersion << 1 // make space for sign bit
   563		if x.a.neg {
   564			b |= 1
   565		}
   566		buf[j] = b
   567		return buf[j:], nil
   568	}
   569	
   570	// GobDecode implements the gob.GobDecoder interface.
   571	func (z *Rat) GobDecode(buf []byte) error {
   572		if len(buf) == 0 {
   573			// Other side sent a nil or default value.
   574			*z = Rat{}
   575			return nil
   576		}
   577		b := buf[0]
   578		if b>>1 != ratGobVersion {
   579			return errors.New(fmt.Sprintf("Rat.GobDecode: encoding version %d not supported", b>>1))
   580		}
   581		const j = 1 + 4
   582		i := j + binary.BigEndian.Uint32(buf[j-4:j])
   583		z.a.neg = b&1 != 0
   584		z.a.abs = z.a.abs.setBytes(buf[j:i])
   585		z.b.abs = z.b.abs.setBytes(buf[i:])
   586		return nil
   587	}
   588	
   589	// MarshalText implements the encoding.TextMarshaler interface
   590	func (r *Rat) MarshalText() (text []byte, err error) {
   591		return []byte(r.RatString()), nil
   592	}
   593	
   594	// UnmarshalText implements the encoding.TextUnmarshaler interface
   595	func (r *Rat) UnmarshalText(text []byte) error {
   596		if _, ok := r.SetString(string(text)); !ok {
   597			return fmt.Errorf("math/big: cannot unmarshal %q into a *big.Rat", text)
   598		}
   599		return nil
   600	}
   601	

View as plain text