...

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

View as plain text