libgo/go/math/big/rat.go

   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         "encoding/binary"
  11         "errors"
  12         "fmt"
  13         "math"
  14         "strings"
  15 )
  16
  17 // A Rat represents a quotient a/b of arbitrary precision.
  18 // The zero value for a Rat represents the value 0.
  19 type Rat struct {
  20         // To make zero values for Rat work w/o initialization,
  21         // a zero value of b (len(b) == 0) acts like b == 1.
  22         // a.neg determines the sign of the Rat, b.neg is ignored.
  23         a, b Int
  24 }
  25
  26 // NewRat creates a new Rat with numerator a and denominator b.
  27 func NewRat(a, b int64) *Rat {
  28         return new(Rat).SetFrac64(a, b)
  29 }
  30
  31 // SetFloat64 sets z to exactly f and returns z.
  32 // If f is not finite, SetFloat returns nil.
  33 func (z *Rat) SetFloat64(f float64) *Rat {
  34         const expMask = 1<<11 - 1
  35         bits := math.Float64bits(f)
  36         mantissa := bits & (1<<52 - 1)
  37         exp := int((bits >> 52) & expMask)
  38         switch exp {
  39         case expMask: // non-finite
  40                 return nil
  41         case 0: // denormal
  42                 exp -= 1022
  43         default: // normal
  44                 mantissa |= 1 << 52
  45                 exp -= 1023
  46         }
  47
  48         shift := 52 - exp
  49
  50         // Optimisation (?): partially pre-normalise.
  51         for mantissa&1 == 0 && shift > 0 {
  52                 mantissa >>= 1
  53                 shift--
  54         }
  55
  56         z.a.SetUint64(mantissa)
  57         z.a.neg = f < 0
  58         z.b.Set(intOne)
  59         if shift > 0 {
  60                 z.b.Lsh(&z.b, uint(shift))
  61         } else {
  62                 z.a.Lsh(&z.a, uint(-shift))
  63         }
  64         return z.norm()
  65 }
  66
  67 // isFinite reports whether f represents a finite rational value.
  68 // It is equivalent to !math.IsNan(f) && !math.IsInf(f, 0).
  69 func isFinite(f float64) bool {
  70         return math.Abs(f) <= math.MaxFloat64
  71 }
  72
  73 // low64 returns the least significant 64 bits of natural number z.
  74 func low64(z nat) uint64 {
  75         if len(z) == 0 {
  76                 return 0
  77         }
  78         if _W == 32 && len(z) > 1 {
  79                 return uint64(z[1])<<32 | uint64(z[0])
  80         }
  81         return uint64(z[0])
  82 }
  83
  84 // quotToFloat returns the non-negative IEEE 754 double-precision
  85 // value nearest to the quotient a/b, using round-to-even in halfway
  86 // cases.  It does not mutate its arguments.
  87 // Preconditions: b is non-zero; a and b have no common factors.
  88 func quotToFloat(a, b nat) (f float64, exact bool) {
  89         // TODO(adonovan): specialize common degenerate cases: 1.0, integers.
  90         alen := a.bitLen()
  91         if alen == 0 {
  92                 return 0, true
  93         }
  94         blen := b.bitLen()
  95         if blen == 0 {
  96                 panic("division by zero")
  97         }
  98
  99         // 1. Left-shift A or B such that quotient A/B is in [1<<53, 1<<55).
 100         // (54 bits if A<B when they are left-aligned, 55 bits if A>=B.)
 101         // This is 2 or 3 more than the float64 mantissa field width of 52:
 102         // - the optional extra bit is shifted away in step 3 below.
 103         // - the high-order 1 is omitted in float64 "normal" representation;
 104         // - the low-order 1 will be used during rounding then discarded.
 105         exp := alen - blen
 106         var a2, b2 nat
 107         a2 = a2.set(a)
 108         b2 = b2.set(b)
 109         if shift := 54 - exp; shift > 0 {
 110                 a2 = a2.shl(a2, uint(shift))
 111         } else if shift < 0 {
 112                 b2 = b2.shl(b2, uint(-shift))
 113         }
 114
 115         // 2. Compute quotient and remainder (q, r).  NB: due to the
 116         // extra shift, the low-order bit of q is logically the
 117         // high-order bit of r.
 118         var q nat
 119         q, r := q.div(a2, a2, b2) // (recycle a2)
 120         mantissa := low64(q)
 121         haveRem := len(r) > 0 // mantissa&1 && !haveRem => remainder is exactly half
 122
 123         // 3. If quotient didn't fit in 54 bits, re-do division by b2<<1
 124         // (in effect---we accomplish this incrementally).
 125         if mantissa>>54 == 1 {
 126                 if mantissa&1 == 1 {
 127                         haveRem = true
 128                 }
 129                 mantissa >>= 1
 130                 exp++
 131         }
 132         if mantissa>>53 != 1 {
 133                 panic("expected exactly 54 bits of result")
 134         }
 135
 136         // 4. Rounding.
 137         if -1022-52 <= exp && exp <= -1022 {
 138                 // Denormal case; lose 'shift' bits of precision.
 139                 shift := uint64(-1022 - (exp - 1)) // [1..53)
 140                 lostbits := mantissa & (1<<shift - 1)
 141                 haveRem = haveRem || lostbits != 0
 142                 mantissa >>= shift
 143                 exp = -1023 + 2
 144         }
 145         // Round q using round-half-to-even.
 146         exact = !haveRem
 147         if mantissa&1 != 0 {
 148                 exact = false
 149                 if haveRem || mantissa&2 != 0 {
 150                         if mantissa++; mantissa >= 1<<54 {
 151                                 // Complete rollover 11...1 => 100...0, so shift is safe
 152                                 mantissa >>= 1
 153                                 exp++
 154                         }
 155                 }
 156         }
 157         mantissa >>= 1 // discard rounding bit.  Mantissa now scaled by 2^53.
 158
 159         f = math.Ldexp(float64(mantissa), exp-53)
 160         if math.IsInf(f, 0) {
 161                 exact = false
 162         }
 163         return
 164 }
 165
 166 // Float64 returns the nearest float64 value for x and a bool indicating
 167 // whether f represents x exactly. If the magnitude of x is too large to
 168 // be represented by a float64, f is an infinity and exact is false.
 169 // The sign of f always matches the sign of x, even if f == 0.
 170 func (x *Rat) Float64() (f float64, exact bool) {
 171         b := x.b.abs
 172         if len(b) == 0 {
 173                 b = b.set(natOne) // materialize denominator
 174         }
 175         f, exact = quotToFloat(x.a.abs, b)
 176         if x.a.neg {
 177                 f = -f
 178         }
 179         return
 180 }
 181
 182 // SetFrac sets z to a/b and returns z.
 183 func (z *Rat) SetFrac(a, b *Int) *Rat {
 184         z.a.neg = a.neg != b.neg
 185         babs := b.abs
 186         if len(babs) == 0 {
 187                 panic("division by zero")
 188         }
 189         if &z.a == b || alias(z.a.abs, babs) {
 190                 babs = nat(nil).set(babs) // make a copy
 191         }
 192         z.a.abs = z.a.abs.set(a.abs)
 193         z.b.abs = z.b.abs.set(babs)
 194         return z.norm()
 195 }
 196
 197 // SetFrac64 sets z to a/b and returns z.
 198 func (z *Rat) SetFrac64(a, b int64) *Rat {
 199         z.a.SetInt64(a)
 200         if b == 0 {
 201                 panic("division by zero")
 202         }
 203         if b < 0 {
 204                 b = -b
 205                 z.a.neg = !z.a.neg
 206         }
 207         z.b.abs = z.b.abs.setUint64(uint64(b))
 208         return z.norm()
 209 }
 210
 211 // SetInt sets z to x (by making a copy of x) and returns z.
 212 func (z *Rat) SetInt(x *Int) *Rat {
 213         z.a.Set(x)
 214         z.b.abs = z.b.abs.make(0)
 215         return z
 216 }
 217
 218 // SetInt64 sets z to x and returns z.
 219 func (z *Rat) SetInt64(x int64) *Rat {
 220         z.a.SetInt64(x)
 221         z.b.abs = z.b.abs.make(0)
 222         return z
 223 }
 224
 225 // Set sets z to x (by making a copy of x) and returns z.
 226 func (z *Rat) Set(x *Rat) *Rat {
 227         if z != x {
 228                 z.a.Set(&x.a)
 229                 z.b.Set(&x.b)
 230         }
 231         return z
 232 }
 233
 234 // Abs sets z to |x| (the absolute value of x) and returns z.
 235 func (z *Rat) Abs(x *Rat) *Rat {
 236         z.Set(x)
 237         z.a.neg = false
 238         return z
 239 }
 240
 241 // Neg sets z to -x and returns z.
 242 func (z *Rat) Neg(x *Rat) *Rat {
 243         z.Set(x)
 244         z.a.neg = len(z.a.abs) > 0 && !z.a.neg // 0 has no sign
 245         return z
 246 }
 247
 248 // Inv sets z to 1/x and returns z.
 249 func (z *Rat) Inv(x *Rat) *Rat {
 250         if len(x.a.abs) == 0 {
 251                 panic("division by zero")
 252         }
 253         z.Set(x)
 254         a := z.b.abs
 255         if len(a) == 0 {
 256                 a = a.set(natOne) // materialize numerator
 257         }
 258         b := z.a.abs
 259         if b.cmp(natOne) == 0 {
 260                 b = b.make(0) // normalize denominator
 261         }
 262         z.a.abs, z.b.abs = a, b // sign doesn't change
 263         return z
 264 }
 265
 266 // Sign returns:
 267 //
 268 //      -1 if x <  0
 269 //       0 if x == 0
 270 //      +1 if x >  0
 271 //
 272 func (x *Rat) Sign() int {
 273         return x.a.Sign()
 274 }
 275
 276 // IsInt returns true if the denominator of x is 1.
 277 func (x *Rat) IsInt() bool {
 278         return len(x.b.abs) == 0 || x.b.abs.cmp(natOne) == 0
 279 }
 280
 281 // Num returns the numerator of x; it may be <= 0.
 282 // The result is a reference to x's numerator; it
 283 // may change if a new value is assigned to x, and vice versa.
 284 // The sign of the numerator corresponds to the sign of x.
 285 func (x *Rat) Num() *Int {
 286         return &x.a
 287 }
 288
 289 // Denom returns the denominator of x; it is always > 0.
 290 // The result is a reference to x's denominator; it
 291 // may change if a new value is assigned to x, and vice versa.
 292 func (x *Rat) Denom() *Int {
 293         x.b.neg = false // the result is always >= 0
 294         if len(x.b.abs) == 0 {
 295                 x.b.abs = x.b.abs.set(natOne) // materialize denominator
 296         }
 297         return &x.b
 298 }
 299
 300 func (z *Rat) norm() *Rat {
 301         switch {
 302         case len(z.a.abs) == 0:
 303                 // z == 0 - normalize sign and denominator
 304                 z.a.neg = false
 305                 z.b.abs = z.b.abs.make(0)
 306         case len(z.b.abs) == 0:
 307                 // z is normalized int - nothing to do
 308         case z.b.abs.cmp(natOne) == 0:
 309                 // z is int - normalize denominator
 310                 z.b.abs = z.b.abs.make(0)
 311         default:
 312                 neg := z.a.neg
 313                 z.a.neg = false
 314                 z.b.neg = false
 315                 if f := NewInt(0).binaryGCD(&z.a, &z.b); f.Cmp(intOne) != 0 {
 316                         z.a.abs, _ = z.a.abs.div(nil, z.a.abs, f.abs)
 317                         z.b.abs, _ = z.b.abs.div(nil, z.b.abs, f.abs)
 318                         if z.b.abs.cmp(natOne) == 0 {
 319                                 // z is int - normalize denominator
 320                                 z.b.abs = z.b.abs.make(0)
 321                         }
 322                 }
 323                 z.a.neg = neg
 324         }
 325         return z
 326 }
 327
 328 // mulDenom sets z to the denominator product x*y (by taking into
 329 // account that 0 values for x or y must be interpreted as 1) and
 330 // returns z.
 331 func mulDenom(z, x, y nat) nat {
 332         switch {
 333         case len(x) == 0:
 334                 return z.set(y)
 335         case len(y) == 0:
 336                 return z.set(x)
 337         }
 338         return z.mul(x, y)
 339 }
 340
 341 // scaleDenom computes x*f.
 342 // If f == 0 (zero value of denominator), the result is (a copy of) x.
 343 func scaleDenom(x *Int, f nat) *Int {
 344         var z Int
 345         if len(f) == 0 {
 346                 return z.Set(x)
 347         }
 348         z.abs = z.abs.mul(x.abs, f)
 349         z.neg = x.neg
 350         return &z
 351 }
 352
 353 // Cmp compares x and y and returns:
 354 //
 355 //   -1 if x <  y
 356 //    0 if x == y
 357 //   +1 if x >  y
 358 //
 359 func (x *Rat) Cmp(y *Rat) int {
 360         return scaleDenom(&x.a, y.b.abs).Cmp(scaleDenom(&y.a, x.b.abs))
 361 }
 362
 363 // Add sets z to the sum x+y and returns z.
 364 func (z *Rat) Add(x, y *Rat) *Rat {
 365         a1 := scaleDenom(&x.a, y.b.abs)
 366         a2 := scaleDenom(&y.a, x.b.abs)
 367         z.a.Add(a1, a2)
 368         z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
 369         return z.norm()
 370 }
 371
 372 // Sub sets z to the difference x-y and returns z.
 373 func (z *Rat) Sub(x, y *Rat) *Rat {
 374         a1 := scaleDenom(&x.a, y.b.abs)
 375         a2 := scaleDenom(&y.a, x.b.abs)
 376         z.a.Sub(a1, a2)
 377         z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
 378         return z.norm()
 379 }
 380
 381 // Mul sets z to the product x*y and returns z.
 382 func (z *Rat) Mul(x, y *Rat) *Rat {
 383         z.a.Mul(&x.a, &y.a)
 384         z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
 385         return z.norm()
 386 }
 387
 388 // Quo sets z to the quotient x/y and returns z.
 389 // If y == 0, a division-by-zero run-time panic occurs.
 390 func (z *Rat) Quo(x, y *Rat) *Rat {
 391         if len(y.a.abs) == 0 {
 392                 panic("division by zero")
 393         }
 394         a := scaleDenom(&x.a, y.b.abs)
 395         b := scaleDenom(&y.a, x.b.abs)
 396         z.a.abs = a.abs
 397         z.b.abs = b.abs
 398         z.a.neg = a.neg != b.neg
 399         return z.norm()
 400 }
 401
 402 func ratTok(ch rune) bool {
 403         return strings.IndexRune("+-/0123456789.eE", ch) >= 0
 404 }
 405
 406 // Scan is a support routine for fmt.Scanner. It accepts the formats
 407 // 'e', 'E', 'f', 'F', 'g', 'G', and 'v'. All formats are equivalent.
 408 func (z *Rat) Scan(s fmt.ScanState, ch rune) error {
 409         tok, err := s.Token(true, ratTok)
 410         if err != nil {
 411                 return err
 412         }
 413         if strings.IndexRune("efgEFGv", ch) < 0 {
 414                 return errors.New("Rat.Scan: invalid verb")
 415         }
 416         if _, ok := z.SetString(string(tok)); !ok {
 417                 return errors.New("Rat.Scan: invalid syntax")
 418         }
 419         return nil
 420 }
 421
 422 // SetString sets z to the value of s and returns z and a boolean indicating
 423 // success. s can be given as a fraction "a/b" or as a floating-point number
 424 // optionally followed by an exponent. If the operation failed, the value of
 425 // z is undefined but the returned value is nil.
 426 func (z *Rat) SetString(s string) (*Rat, bool) {
 427         if len(s) == 0 {
 428                 return nil, false
 429         }
 430
 431         // check for a quotient
 432         sep := strings.Index(s, "/")
 433         if sep >= 0 {
 434                 if _, ok := z.a.SetString(s[0:sep], 10); !ok {
 435                         return nil, false
 436                 }
 437                 s = s[sep+1:]
 438                 var err error
 439                 if z.b.abs, _, err = z.b.abs.scan(strings.NewReader(s), 10); err != nil {
 440                         return nil, false
 441                 }
 442                 return z.norm(), true
 443         }
 444
 445         // check for a decimal point
 446         sep = strings.Index(s, ".")
 447         // check for an exponent
 448         e := strings.IndexAny(s, "eE")
 449         var exp Int
 450         if e >= 0 {
 451                 if e < sep {
 452                         // The E must come after the decimal point.
 453                         return nil, false
 454                 }
 455                 if _, ok := exp.SetString(s[e+1:], 10); !ok {
 456                         return nil, false
 457                 }
 458                 s = s[0:e]
 459         }
 460         if sep >= 0 {
 461                 s = s[0:sep] + s[sep+1:]
 462                 exp.Sub(&exp, NewInt(int64(len(s)-sep)))
 463         }
 464
 465         if _, ok := z.a.SetString(s, 10); !ok {
 466                 return nil, false
 467         }
 468         powTen := nat(nil).expNN(natTen, exp.abs, nil)
 469         if exp.neg {
 470                 z.b.abs = powTen
 471                 z.norm()
 472         } else {
 473                 z.a.abs = z.a.abs.mul(z.a.abs, powTen)
 474                 z.b.abs = z.b.abs.make(0)
 475         }
 476
 477         return z, true
 478 }
 479
 480 // String returns a string representation of z in the form "a/b" (even if b == 1).
 481 func (x *Rat) String() string {
 482         s := "/1"
 483         if len(x.b.abs) != 0 {
 484                 s = "/" + x.b.abs.decimalString()
 485         }
 486         return x.a.String() + s
 487 }
 488
 489 // RatString returns a string representation of z in the form "a/b" if b != 1,
 490 // and in the form "a" if b == 1.
 491 func (x *Rat) RatString() string {
 492         if x.IsInt() {
 493                 return x.a.String()
 494         }
 495         return x.String()
 496 }
 497
 498 // FloatString returns a string representation of z in decimal form with prec
 499 // digits of precision after the decimal point and the last digit rounded.
 500 func (x *Rat) FloatString(prec int) string {
 501         if x.IsInt() {
 502                 s := x.a.String()
 503                 if prec > 0 {
 504                         s += "." + strings.Repeat("0", prec)
 505                 }
 506                 return s
 507         }
 508         // x.b.abs != 0
 509
 510         q, r := nat(nil).div(nat(nil), x.a.abs, x.b.abs)
 511
 512         p := natOne
 513         if prec > 0 {
 514                 p = nat(nil).expNN(natTen, nat(nil).setUint64(uint64(prec)), nil)
 515         }
 516
 517         r = r.mul(r, p)
 518         r, r2 := r.div(nat(nil), r, x.b.abs)
 519
 520         // see if we need to round up
 521         r2 = r2.add(r2, r2)
 522         if x.b.abs.cmp(r2) <= 0 {
 523                 r = r.add(r, natOne)
 524                 if r.cmp(p) >= 0 {
 525                         q = nat(nil).add(q, natOne)
 526                         r = nat(nil).sub(r, p)
 527                 }
 528         }
 529
 530         s := q.decimalString()
 531         if x.a.neg {
 532                 s = "-" + s
 533         }
 534
 535         if prec > 0 {
 536                 rs := r.decimalString()
 537                 leadingZeros := prec - len(rs)
 538                 s += "." + strings.Repeat("0", leadingZeros) + rs
 539         }
 540
 541         return s
 542 }
 543
 544 // Gob codec version. Permits backward-compatible changes to the encoding.
 545 const ratGobVersion byte = 1
 546
 547 // GobEncode implements the gob.GobEncoder interface.
 548 func (x *Rat) GobEncode() ([]byte, error) {
 549         if x == nil {
 550                 return nil, nil
 551         }
 552         buf := make([]byte, 1+4+(len(x.a.abs)+len(x.b.abs))*_S) // extra bytes for version and sign bit (1), and numerator length (4)
 553         i := x.b.abs.bytes(buf)
 554         j := x.a.abs.bytes(buf[0:i])
 555         n := i - j
 556         if int(uint32(n)) != n {
 557                 // this should never happen
 558                 return nil, errors.New("Rat.GobEncode: numerator too large")
 559         }
 560         binary.BigEndian.PutUint32(buf[j-4:j], uint32(n))
 561         j -= 1 + 4
 562         b := ratGobVersion << 1 // make space for sign bit
 563         if x.a.neg {
 564                 b |= 1
 565         }
 566         buf[j] = b
 567         return buf[j:], nil
 568 }
 569
 570 // GobDecode implements the gob.GobDecoder interface.
 571 func (z *Rat) GobDecode(buf []byte) error {
 572         if len(buf) == 0 {
 573                 // Other side sent a nil or default value.
 574                 *z = Rat{}
 575                 return nil
 576         }
 577         b := buf[0]
 578         if b>>1 != ratGobVersion {
 579                 return errors.New(fmt.Sprintf("Rat.GobDecode: encoding version %d not supported", b>>1))
 580         }
 581         const j = 1 + 4
 582         i := j + binary.BigEndian.Uint32(buf[j-4:j])
 583         z.a.neg = b&1 != 0
 584         z.a.abs = z.a.abs.setBytes(buf[j:i])
 585         z.b.abs = z.b.abs.setBytes(buf[i:])
 586         return nil
 587 }