1// Copyright 2010 The Go Authors. All rights reserved.
2// Use of this source code is governed by a BSD-style
3// license that can be found in the LICENSE file.
4
5// This file implements multi-precision rational numbers.
6
7package big
8
9import (
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.
23type 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.
33func 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.
39func (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 {
45	case expMask: // non-finite
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.
77func 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.
175func 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.
273func (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.
289func (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.
303func (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.
319func (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.
333func (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.
340func (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.
347func (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.
354func (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.
366func (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.
373func (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.
381func (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//   - -1 if x < 0;
392//   - 0 if x == 0;
393//   - +1 if x > 0.
394func (x *Rat) Sign() int {
395	return x.a.Sign()
396}
397
398// IsInt reports whether the denominator of x is 1.
399func (x *Rat) IsInt() bool {
400	return len(x.b.abs) == 0 || x.b.abs.cmp(natOne) == 0
401}
402
403// Num returns the numerator of x; it may be <= 0.
404// The result is a reference to x's numerator; it
405// may change if a new value is assigned to x, and vice versa.
406// The sign of the numerator corresponds to the sign of x.
407func (x *Rat) Num() *Int {
408	return &x.a
409}
410
411// Denom returns the denominator of x; it is always > 0.
412// The result is a reference to x's denominator, unless
413// x is an uninitialized (zero value) [Rat], in which case
414// the result is a new [Int] of value 1. (To initialize x,
415// any operation that sets x will do, including x.Set(x).)
416// If the result is a reference to x's denominator it
417// may change if a new value is assigned to x, and vice versa.
418func (x *Rat) Denom() *Int {
419	// Note that x.b.neg is guaranteed false.
420	if len(x.b.abs) == 0 {
421		// Note: If this proves problematic, we could
422		//       panic instead and require the Rat to
423		//       be explicitly initialized.
424		return &Int{abs: nat{1}}
425	}
426	return &x.b
427}
428
429func (z *Rat) norm() *Rat {
430	switch {
431	case len(z.a.abs) == 0:
432		// z == 0; normalize sign and denominator
433		z.a.neg = false
434		fallthrough
435	case len(z.b.abs) == 0:
436		// z is integer; normalize denominator
437		z.b.abs = z.b.abs.setWord(1)
438	default:
439		// z is fraction; normalize numerator and denominator
440		neg := z.a.neg
441		z.a.neg = false
442		z.b.neg = false
443		if f := NewInt(0).lehmerGCD(nil, nil, &z.a, &z.b); f.Cmp(intOne) != 0 {
444			z.a.abs, _ = z.a.abs.div(nil, z.a.abs, f.abs)
445			z.b.abs, _ = z.b.abs.div(nil, z.b.abs, f.abs)
446		}
447		z.a.neg = neg
448	}
449	return z
450}
451
452// mulDenom sets z to the denominator product x*y (by taking into
453// account that 0 values for x or y must be interpreted as 1) and
454// returns z.
455func mulDenom(z, x, y nat) nat {
456	switch {
457	case len(x) == 0 && len(y) == 0:
458		return z.setWord(1)
459	case len(x) == 0:
460		return z.set(y)
461	case len(y) == 0:
462		return z.set(x)
463	}
464	return z.mul(x, y)
465}
466
467// scaleDenom sets z to the product x*f.
468// If f == 0 (zero value of denominator), z is set to (a copy of) x.
469func (z *Int) scaleDenom(x *Int, f nat) {
470	if len(f) == 0 {
471		z.Set(x)
472		return
473	}
474	z.abs = z.abs.mul(x.abs, f)
475	z.neg = x.neg
476}
477
478// Cmp compares x and y and returns:
479//   - -1 if x < y;
480//   - 0 if x == y;
481//   - +1 if x > y.
482func (x *Rat) Cmp(y *Rat) int {
483	var a, b Int
484	a.scaleDenom(&x.a, y.b.abs)
485	b.scaleDenom(&y.a, x.b.abs)
486	return a.Cmp(&b)
487}
488
489// Add sets z to the sum x+y and returns z.
490func (z *Rat) Add(x, y *Rat) *Rat {
491	var a1, a2 Int
492	a1.scaleDenom(&x.a, y.b.abs)
493	a2.scaleDenom(&y.a, x.b.abs)
494	z.a.Add(&a1, &a2)
495	z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
496	return z.norm()
497}
498
499// Sub sets z to the difference x-y and returns z.
500func (z *Rat) Sub(x, y *Rat) *Rat {
501	var a1, a2 Int
502	a1.scaleDenom(&x.a, y.b.abs)
503	a2.scaleDenom(&y.a, x.b.abs)
504	z.a.Sub(&a1, &a2)
505	z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
506	return z.norm()
507}
508
509// Mul sets z to the product x*y and returns z.
510func (z *Rat) Mul(x, y *Rat) *Rat {
511	if x == y {
512		// a squared Rat is positive and can't be reduced (no need to call norm())
513		z.a.neg = false
514		z.a.abs = z.a.abs.sqr(x.a.abs)
515		if len(x.b.abs) == 0 {
516			z.b.abs = z.b.abs.setWord(1)
517		} else {
518			z.b.abs = z.b.abs.sqr(x.b.abs)
519		}
520		return z
521	}
522	z.a.Mul(&x.a, &y.a)
523	z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
524	return z.norm()
525}
526
527// Quo sets z to the quotient x/y and returns z.
528// If y == 0, Quo panics.
529func (z *Rat) Quo(x, y *Rat) *Rat {
530	if len(y.a.abs) == 0 {
531		panic("division by zero")
532	}
533	var a, b Int
534	a.scaleDenom(&x.a, y.b.abs)
535	b.scaleDenom(&y.a, x.b.abs)
536	z.a.abs = a.abs
537	z.b.abs = b.abs
538	z.a.neg = a.neg != b.neg
539	return z.norm()
540}
541