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