...

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

View as plain text