...

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

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

View as plain text