Black Lives Matter. Support the Equal Justice Initiative.

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

## Documentation: math/big

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

View as plain text