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 // SetUint64 sets z to x and returns z. 343 func (z *Rat) SetUint64(x uint64) *Rat { 344 z.a.SetUint64(x) 345 z.b.abs = z.b.abs[:0] 346 return z 347 } 348 349 // Set sets z to x (by making a copy of x) and returns z. 350 func (z *Rat) Set(x *Rat) *Rat { 351 if z != x { 352 z.a.Set(&x.a) 353 z.b.Set(&x.b) 354 } 355 return z 356 } 357 358 // Abs sets z to |x| (the absolute value of x) and returns z. 359 func (z *Rat) Abs(x *Rat) *Rat { 360 z.Set(x) 361 z.a.neg = false 362 return z 363 } 364 365 // Neg sets z to -x and returns z. 366 func (z *Rat) Neg(x *Rat) *Rat { 367 z.Set(x) 368 z.a.neg = len(z.a.abs) > 0 && !z.a.neg // 0 has no sign 369 return z 370 } 371 372 // Inv sets z to 1/x and returns z. 373 func (z *Rat) Inv(x *Rat) *Rat { 374 if len(x.a.abs) == 0 { 375 panic("division by zero") 376 } 377 z.Set(x) 378 a := z.b.abs 379 if len(a) == 0 { 380 a = a.set(natOne) // materialize numerator 381 } 382 b := z.a.abs 383 if b.cmp(natOne) == 0 { 384 b = b[:0] // normalize denominator 385 } 386 z.a.abs, z.b.abs = a, b // sign doesn't change 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; it 415 // may change if a new value is assigned to x, and vice versa. 416 func (x *Rat) Denom() *Int { 417 x.b.neg = false // the result is always >= 0 418 if len(x.b.abs) == 0 { 419 x.b.abs = x.b.abs.set(natOne) // materialize denominator 420 } 421 return &x.b 422 } 423 424 func (z *Rat) norm() *Rat { 425 switch { 426 case len(z.a.abs) == 0: 427 // z == 0 - normalize sign and denominator 428 z.a.neg = false 429 z.b.abs = z.b.abs[:0] 430 case len(z.b.abs) == 0: 431 // z is normalized int - nothing to do 432 case z.b.abs.cmp(natOne) == 0: 433 // z is int - normalize denominator 434 z.b.abs = z.b.abs[:0] 435 default: 436 neg := z.a.neg 437 z.a.neg = false 438 z.b.neg = false 439 if f := NewInt(0).lehmerGCD(nil, nil, &z.a, &z.b); f.Cmp(intOne) != 0 { 440 z.a.abs, _ = z.a.abs.div(nil, z.a.abs, f.abs) 441 z.b.abs, _ = z.b.abs.div(nil, z.b.abs, f.abs) 442 if z.b.abs.cmp(natOne) == 0 { 443 // z is int - normalize denominator 444 z.b.abs = z.b.abs[:0] 445 } 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. 455 func mulDenom(z, x, y nat) nat { 456 switch { 457 case len(x) == 0: 458 return z.set(y) 459 case len(y) == 0: 460 return z.set(x) 461 } 462 return z.mul(x, y) 463 } 464 465 // scaleDenom sets z to the product x*f. 466 // If f == 0 (zero value of denominator), z is set to (a copy of) x. 467 func (z *Int) scaleDenom(x *Int, f nat) { 468 if len(f) == 0 { 469 z.Set(x) 470 return 471 } 472 z.abs = z.abs.mul(x.abs, f) 473 z.neg = x.neg 474 } 475 476 // Cmp compares x and y and returns: 477 // 478 // -1 if x < y 479 // 0 if x == y 480 // +1 if x > y 481 // 482 func (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. 490 func (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. 500 func (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. 510 func (z *Rat) Mul(x, y *Rat) *Rat { 511 if x == y { 512 // a squared Rat is positive and can't be reduced 513 z.a.neg = false 514 z.a.abs = z.a.abs.sqr(x.a.abs) 515 z.b.abs = z.b.abs.sqr(x.b.abs) 516 return z 517 } 518 z.a.Mul(&x.a, &y.a) 519 z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs) 520 return z.norm() 521 } 522 523 // Quo sets z to the quotient x/y and returns z. 524 // If y == 0, a division-by-zero run-time panic occurs. 525 func (z *Rat) Quo(x, y *Rat) *Rat { 526 if len(y.a.abs) == 0 { 527 panic("division by zero") 528 } 529 var a, b Int 530 a.scaleDenom(&x.a, y.b.abs) 531 b.scaleDenom(&y.a, x.b.abs) 532 z.a.abs = a.abs 533 z.b.abs = b.abs 534 z.a.neg = a.neg != b.neg 535 return z.norm() 536 } 537
View as plain text