...

# 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
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 {
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)
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
522		if x.b.abs.cmp(r2) <= 0 {
524			if r.cmp(p) >= 0 {
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