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. 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 {
		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.
		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)
	 498  	z.a.Add(&a1, &a2)
	 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