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  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 {
    36  	case expMask: // non-finite
    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)
   477  	z.a.Add(a1, a2)
   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