...

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

View as plain text