1 // Copyright 2014 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 floating-point numbers. 6 // Like in the GNU MPFR library (https://www.mpfr.org/), operands 7 // can be of mixed precision. Unlike MPFR, the rounding mode is 8 // not specified with each operation, but with each operand. The 9 // rounding mode of the result operand determines the rounding 10 // mode of an operation. This is a from-scratch implementation. 11 12 package big 13 14 import ( 15 "fmt" 16 "math" 17 "math/bits" 18 ) 19 20 const debugFloat = false // enable for debugging 21 22 // A nonzero finite Float represents a multi-precision floating point number 23 // 24 // sign × mantissa × 2**exponent 25 // 26 // with 0.5 <= mantissa < 1.0, and MinExp <= exponent <= MaxExp. 27 // A Float may also be zero (+0, -0) or infinite (+Inf, -Inf). 28 // All Floats are ordered, and the ordering of two Floats x and y 29 // is defined by x.Cmp(y). 30 // 31 // Each Float value also has a precision, rounding mode, and accuracy. 32 // The precision is the maximum number of mantissa bits available to 33 // represent the value. The rounding mode specifies how a result should 34 // be rounded to fit into the mantissa bits, and accuracy describes the 35 // rounding error with respect to the exact result. 36 // 37 // Unless specified otherwise, all operations (including setters) that 38 // specify a *Float variable for the result (usually via the receiver 39 // with the exception of MantExp), round the numeric result according 40 // to the precision and rounding mode of the result variable. 41 // 42 // If the provided result precision is 0 (see below), it is set to the 43 // precision of the argument with the largest precision value before any 44 // rounding takes place, and the rounding mode remains unchanged. Thus, 45 // uninitialized Floats provided as result arguments will have their 46 // precision set to a reasonable value determined by the operands, and 47 // their mode is the zero value for RoundingMode (ToNearestEven). 48 // 49 // By setting the desired precision to 24 or 53 and using matching rounding 50 // mode (typically ToNearestEven), Float operations produce the same results 51 // as the corresponding float32 or float64 IEEE-754 arithmetic for operands 52 // that correspond to normal (i.e., not denormal) float32 or float64 numbers. 53 // Exponent underflow and overflow lead to a 0 or an Infinity for different 54 // values than IEEE-754 because Float exponents have a much larger range. 55 // 56 // The zero (uninitialized) value for a Float is ready to use and represents 57 // the number +0.0 exactly, with precision 0 and rounding mode ToNearestEven. 58 // 59 // Operations always take pointer arguments (*Float) rather 60 // than Float values, and each unique Float value requires 61 // its own unique *Float pointer. To "copy" a Float value, 62 // an existing (or newly allocated) Float must be set to 63 // a new value using the Float.Set method; shallow copies 64 // of Floats are not supported and may lead to errors. 65 type Float struct { 66 prec uint32 67 mode RoundingMode 68 acc Accuracy 69 form form 70 neg bool 71 mant nat 72 exp int32 73 } 74 75 // An ErrNaN panic is raised by a Float operation that would lead to 76 // a NaN under IEEE-754 rules. An ErrNaN implements the error interface. 77 type ErrNaN struct { 78 msg string 79 } 80 81 func (err ErrNaN) Error() string { 82 return err.msg 83 } 84 85 // NewFloat allocates and returns a new Float set to x, 86 // with precision 53 and rounding mode ToNearestEven. 87 // NewFloat panics with ErrNaN if x is a NaN. 88 func NewFloat(x float64) *Float { 89 if math.IsNaN(x) { 90 panic(ErrNaN{"NewFloat(NaN)"}) 91 } 92 return new(Float).SetFloat64(x) 93 } 94 95 // Exponent and precision limits. 96 const ( 97 MaxExp = math.MaxInt32 // largest supported exponent 98 MinExp = math.MinInt32 // smallest supported exponent 99 MaxPrec = math.MaxUint32 // largest (theoretically) supported precision; likely memory-limited 100 ) 101 102 // Internal representation: The mantissa bits x.mant of a nonzero finite 103 // Float x are stored in a nat slice long enough to hold up to x.prec bits; 104 // the slice may (but doesn't have to) be shorter if the mantissa contains 105 // trailing 0 bits. x.mant is normalized if the msb of x.mant == 1 (i.e., 106 // the msb is shifted all the way "to the left"). Thus, if the mantissa has 107 // trailing 0 bits or x.prec is not a multiple of the Word size _W, 108 // x.mant[0] has trailing zero bits. The msb of the mantissa corresponds 109 // to the value 0.5; the exponent x.exp shifts the binary point as needed. 110 // 111 // A zero or non-finite Float x ignores x.mant and x.exp. 112 // 113 // x form neg mant exp 114 // ---------------------------------------------------------- 115 // ±0 zero sign - - 116 // 0 < |x| < +Inf finite sign mantissa exponent 117 // ±Inf inf sign - - 118 119 // A form value describes the internal representation. 120 type form byte 121 122 // The form value order is relevant - do not change! 123 const ( 124 zero form = iota 125 finite 126 inf 127 ) 128 129 // RoundingMode determines how a Float value is rounded to the 130 // desired precision. Rounding may change the Float value; the 131 // rounding error is described by the Float's Accuracy. 132 type RoundingMode byte 133 134 // These constants define supported rounding modes. 135 const ( 136 ToNearestEven RoundingMode = iota // == IEEE 754-2008 roundTiesToEven 137 ToNearestAway // == IEEE 754-2008 roundTiesToAway 138 ToZero // == IEEE 754-2008 roundTowardZero 139 AwayFromZero // no IEEE 754-2008 equivalent 140 ToNegativeInf // == IEEE 754-2008 roundTowardNegative 141 ToPositiveInf // == IEEE 754-2008 roundTowardPositive 142 ) 143 144 //go:generate stringer -type=RoundingMode 145 146 // Accuracy describes the rounding error produced by the most recent 147 // operation that generated a Float value, relative to the exact value. 148 type Accuracy int8 149 150 // Constants describing the Accuracy of a Float. 151 const ( 152 Below Accuracy = -1 153 Exact Accuracy = 0 154 Above Accuracy = +1 155 ) 156 157 //go:generate stringer -type=Accuracy 158 159 // SetPrec sets z's precision to prec and returns the (possibly) rounded 160 // value of z. Rounding occurs according to z's rounding mode if the mantissa 161 // cannot be represented in prec bits without loss of precision. 162 // SetPrec(0) maps all finite values to ±0; infinite values remain unchanged. 163 // If prec > MaxPrec, it is set to MaxPrec. 164 func (z *Float) SetPrec(prec uint) *Float { 165 z.acc = Exact // optimistically assume no rounding is needed 166 167 // special case 168 if prec == 0 { 169 z.prec = 0 170 if z.form == finite { 171 // truncate z to 0 172 z.acc = makeAcc(z.neg) 173 z.form = zero 174 } 175 return z 176 } 177 178 // general case 179 if prec > MaxPrec { 180 prec = MaxPrec 181 } 182 old := z.prec 183 z.prec = uint32(prec) 184 if z.prec < old { 185 z.round(0) 186 } 187 return z 188 } 189 190 func makeAcc(above bool) Accuracy { 191 if above { 192 return Above 193 } 194 return Below 195 } 196 197 // SetMode sets z's rounding mode to mode and returns an exact z. 198 // z remains unchanged otherwise. 199 // z.SetMode(z.Mode()) is a cheap way to set z's accuracy to Exact. 200 func (z *Float) SetMode(mode RoundingMode) *Float { 201 z.mode = mode 202 z.acc = Exact 203 return z 204 } 205 206 // Prec returns the mantissa precision of x in bits. 207 // The result may be 0 for |x| == 0 and |x| == Inf. 208 func (x *Float) Prec() uint { 209 return uint(x.prec) 210 } 211 212 // MinPrec returns the minimum precision required to represent x exactly 213 // (i.e., the smallest prec before x.SetPrec(prec) would start rounding x). 214 // The result is 0 for |x| == 0 and |x| == Inf. 215 func (x *Float) MinPrec() uint { 216 if x.form != finite { 217 return 0 218 } 219 return uint(len(x.mant))*_W - x.mant.trailingZeroBits() 220 } 221 222 // Mode returns the rounding mode of x. 223 func (x *Float) Mode() RoundingMode { 224 return x.mode 225 } 226 227 // Acc returns the accuracy of x produced by the most recent 228 // operation, unless explicitly documented otherwise by that 229 // operation. 230 func (x *Float) Acc() Accuracy { 231 return x.acc 232 } 233 234 // Sign returns: 235 // 236 // -1 if x < 0 237 // 0 if x is ±0 238 // +1 if x > 0 239 func (x *Float) Sign() int { 240 if debugFloat { 241 x.validate() 242 } 243 if x.form == zero { 244 return 0 245 } 246 if x.neg { 247 return -1 248 } 249 return 1 250 } 251 252 // MantExp breaks x into its mantissa and exponent components 253 // and returns the exponent. If a non-nil mant argument is 254 // provided its value is set to the mantissa of x, with the 255 // same precision and rounding mode as x. The components 256 // satisfy x == mant × 2**exp, with 0.5 <= |mant| < 1.0. 257 // Calling MantExp with a nil argument is an efficient way to 258 // get the exponent of the receiver. 259 // 260 // Special cases are: 261 // 262 // ( ±0).MantExp(mant) = 0, with mant set to ±0 263 // (±Inf).MantExp(mant) = 0, with mant set to ±Inf 264 // 265 // x and mant may be the same in which case x is set to its 266 // mantissa value. 267 func (x *Float) MantExp(mant *Float) (exp int) { 268 if debugFloat { 269 x.validate() 270 } 271 if x.form == finite { 272 exp = int(x.exp) 273 } 274 if mant != nil { 275 mant.Copy(x) 276 if mant.form == finite { 277 mant.exp = 0 278 } 279 } 280 return 281 } 282 283 func (z *Float) setExpAndRound(exp int64, sbit uint) { 284 if exp < MinExp { 285 // underflow 286 z.acc = makeAcc(z.neg) 287 z.form = zero 288 return 289 } 290 291 if exp > MaxExp { 292 // overflow 293 z.acc = makeAcc(!z.neg) 294 z.form = inf 295 return 296 } 297 298 z.form = finite 299 z.exp = int32(exp) 300 z.round(sbit) 301 } 302 303 // SetMantExp sets z to mant × 2**exp and returns z. 304 // The result z has the same precision and rounding mode 305 // as mant. SetMantExp is an inverse of MantExp but does 306 // not require 0.5 <= |mant| < 1.0. Specifically, for a 307 // given x of type *Float, SetMantExp relates to MantExp 308 // as follows: 309 // 310 // mant := new(Float) 311 // new(Float).SetMantExp(mant, x.MantExp(mant)).Cmp(x) == 0 312 // 313 // Special cases are: 314 // 315 // z.SetMantExp( ±0, exp) = ±0 316 // z.SetMantExp(±Inf, exp) = ±Inf 317 // 318 // z and mant may be the same in which case z's exponent 319 // is set to exp. 320 func (z *Float) SetMantExp(mant *Float, exp int) *Float { 321 if debugFloat { 322 z.validate() 323 mant.validate() 324 } 325 z.Copy(mant) 326 327 if z.form == finite { 328 // 0 < |mant| < +Inf 329 z.setExpAndRound(int64(z.exp)+int64(exp), 0) 330 } 331 return z 332 } 333 334 // Signbit reports whether x is negative or negative zero. 335 func (x *Float) Signbit() bool { 336 return x.neg 337 } 338 339 // IsInf reports whether x is +Inf or -Inf. 340 func (x *Float) IsInf() bool { 341 return x.form == inf 342 } 343 344 // IsInt reports whether x is an integer. 345 // ±Inf values are not integers. 346 func (x *Float) IsInt() bool { 347 if debugFloat { 348 x.validate() 349 } 350 // special cases 351 if x.form != finite { 352 return x.form == zero 353 } 354 // x.form == finite 355 if x.exp <= 0 { 356 return false 357 } 358 // x.exp > 0 359 return x.prec <= uint32(x.exp) || x.MinPrec() <= uint(x.exp) // not enough bits for fractional mantissa 360 } 361 362 // debugging support 363 func (x *Float) validate() { 364 if !debugFloat { 365 // avoid performance bugs 366 panic("validate called but debugFloat is not set") 367 } 368 if x.form != finite { 369 return 370 } 371 m := len(x.mant) 372 if m == 0 { 373 panic("nonzero finite number with empty mantissa") 374 } 375 const msb = 1 << (_W - 1) 376 if x.mant[m-1]&msb == 0 { 377 panic(fmt.Sprintf("msb not set in last word %#x of %s", x.mant[m-1], x.Text('p', 0))) 378 } 379 if x.prec == 0 { 380 panic("zero precision finite number") 381 } 382 } 383 384 // round rounds z according to z.mode to z.prec bits and sets z.acc accordingly. 385 // sbit must be 0 or 1 and summarizes any "sticky bit" information one might 386 // have before calling round. z's mantissa must be normalized (with the msb set) 387 // or empty. 388 // 389 // CAUTION: The rounding modes ToNegativeInf, ToPositiveInf are affected by the 390 // sign of z. For correct rounding, the sign of z must be set correctly before 391 // calling round. 392 func (z *Float) round(sbit uint) { 393 if debugFloat { 394 z.validate() 395 } 396 397 z.acc = Exact 398 if z.form != finite { 399 // ±0 or ±Inf => nothing left to do 400 return 401 } 402 // z.form == finite && len(z.mant) > 0 403 // m > 0 implies z.prec > 0 (checked by validate) 404 405 m := uint32(len(z.mant)) // present mantissa length in words 406 bits := m * _W // present mantissa bits; bits > 0 407 if bits <= z.prec { 408 // mantissa fits => nothing to do 409 return 410 } 411 // bits > z.prec 412 413 // Rounding is based on two bits: the rounding bit (rbit) and the 414 // sticky bit (sbit). The rbit is the bit immediately before the 415 // z.prec leading mantissa bits (the "0.5"). The sbit is set if any 416 // of the bits before the rbit are set (the "0.25", "0.125", etc.): 417 // 418 // rbit sbit => "fractional part" 419 // 420 // 0 0 == 0 421 // 0 1 > 0 , < 0.5 422 // 1 0 == 0.5 423 // 1 1 > 0.5, < 1.0 424 425 // bits > z.prec: mantissa too large => round 426 r := uint(bits - z.prec - 1) // rounding bit position; r >= 0 427 rbit := z.mant.bit(r) & 1 // rounding bit; be safe and ensure it's a single bit 428 // The sticky bit is only needed for rounding ToNearestEven 429 // or when the rounding bit is zero. Avoid computation otherwise. 430 if sbit == 0 && (rbit == 0 || z.mode == ToNearestEven) { 431 sbit = z.mant.sticky(r) 432 } 433 sbit &= 1 // be safe and ensure it's a single bit 434 435 // cut off extra words 436 n := (z.prec + (_W - 1)) / _W // mantissa length in words for desired precision 437 if m > n { 438 copy(z.mant, z.mant[m-n:]) // move n last words to front 439 z.mant = z.mant[:n] 440 } 441 442 // determine number of trailing zero bits (ntz) and compute lsb mask of mantissa's least-significant word 443 ntz := n*_W - z.prec // 0 <= ntz < _W 444 lsb := Word(1) << ntz 445 446 // round if result is inexact 447 if rbit|sbit != 0 { 448 // Make rounding decision: The result mantissa is truncated ("rounded down") 449 // by default. Decide if we need to increment, or "round up", the (unsigned) 450 // mantissa. 451 inc := false 452 switch z.mode { 453 case ToNegativeInf: 454 inc = z.neg 455 case ToZero: 456 // nothing to do 457 case ToNearestEven: 458 inc = rbit != 0 && (sbit != 0 || z.mant[0]&lsb != 0) 459 case ToNearestAway: 460 inc = rbit != 0 461 case AwayFromZero: 462 inc = true 463 case ToPositiveInf: 464 inc = !z.neg 465 default: 466 panic("unreachable") 467 } 468 469 // A positive result (!z.neg) is Above the exact result if we increment, 470 // and it's Below if we truncate (Exact results require no rounding). 471 // For a negative result (z.neg) it is exactly the opposite. 472 z.acc = makeAcc(inc != z.neg) 473 474 if inc { 475 // add 1 to mantissa 476 if addVW(z.mant, z.mant, lsb) != 0 { 477 // mantissa overflow => adjust exponent 478 if z.exp >= MaxExp { 479 // exponent overflow 480 z.form = inf 481 return 482 } 483 z.exp++ 484 // adjust mantissa: divide by 2 to compensate for exponent adjustment 485 shrVU(z.mant, z.mant, 1) 486 // set msb == carry == 1 from the mantissa overflow above 487 const msb = 1 << (_W - 1) 488 z.mant[n-1] |= msb 489 } 490 } 491 } 492 493 // zero out trailing bits in least-significant word 494 z.mant[0] &^= lsb - 1 495 496 if debugFloat { 497 z.validate() 498 } 499 } 500 501 func (z *Float) setBits64(neg bool, x uint64) *Float { 502 if z.prec == 0 { 503 z.prec = 64 504 } 505 z.acc = Exact 506 z.neg = neg 507 if x == 0 { 508 z.form = zero 509 return z 510 } 511 // x != 0 512 z.form = finite 513 s := bits.LeadingZeros64(x) 514 z.mant = z.mant.setUint64(x << uint(s)) 515 z.exp = int32(64 - s) // always fits 516 if z.prec < 64 { 517 z.round(0) 518 } 519 return z 520 } 521 522 // SetUint64 sets z to the (possibly rounded) value of x and returns z. 523 // If z's precision is 0, it is changed to 64 (and rounding will have 524 // no effect). 525 func (z *Float) SetUint64(x uint64) *Float { 526 return z.setBits64(false, x) 527 } 528 529 // SetInt64 sets z to the (possibly rounded) value of x and returns z. 530 // If z's precision is 0, it is changed to 64 (and rounding will have 531 // no effect). 532 func (z *Float) SetInt64(x int64) *Float { 533 u := x 534 if u < 0 { 535 u = -u 536 } 537 // We cannot simply call z.SetUint64(uint64(u)) and change 538 // the sign afterwards because the sign affects rounding. 539 return z.setBits64(x < 0, uint64(u)) 540 } 541 542 // SetFloat64 sets z to the (possibly rounded) value of x and returns z. 543 // If z's precision is 0, it is changed to 53 (and rounding will have 544 // no effect). SetFloat64 panics with ErrNaN if x is a NaN. 545 func (z *Float) SetFloat64(x float64) *Float { 546 if z.prec == 0 { 547 z.prec = 53 548 } 549 if math.IsNaN(x) { 550 panic(ErrNaN{"Float.SetFloat64(NaN)"}) 551 } 552 z.acc = Exact 553 z.neg = math.Signbit(x) // handle -0, -Inf correctly 554 if x == 0 { 555 z.form = zero 556 return z 557 } 558 if math.IsInf(x, 0) { 559 z.form = inf 560 return z 561 } 562 // normalized x != 0 563 z.form = finite 564 fmant, exp := math.Frexp(x) // get normalized mantissa 565 z.mant = z.mant.setUint64(1<<63 | math.Float64bits(fmant)<<11) 566 z.exp = int32(exp) // always fits 567 if z.prec < 53 { 568 z.round(0) 569 } 570 return z 571 } 572 573 // fnorm normalizes mantissa m by shifting it to the left 574 // such that the msb of the most-significant word (msw) is 1. 575 // It returns the shift amount. It assumes that len(m) != 0. 576 func fnorm(m nat) int64 { 577 if debugFloat && (len(m) == 0 || m[len(m)-1] == 0) { 578 panic("msw of mantissa is 0") 579 } 580 s := nlz(m[len(m)-1]) 581 if s > 0 { 582 c := shlVU(m, m, s) 583 if debugFloat && c != 0 { 584 panic("nlz or shlVU incorrect") 585 } 586 } 587 return int64(s) 588 } 589 590 // SetInt sets z to the (possibly rounded) value of x and returns z. 591 // If z's precision is 0, it is changed to the larger of x.BitLen() 592 // or 64 (and rounding will have no effect). 593 func (z *Float) SetInt(x *Int) *Float { 594 // TODO(gri) can be more efficient if z.prec > 0 595 // but small compared to the size of x, or if there 596 // are many trailing 0's. 597 bits := uint32(x.BitLen()) 598 if z.prec == 0 { 599 z.prec = umax32(bits, 64) 600 } 601 z.acc = Exact 602 z.neg = x.neg 603 if len(x.abs) == 0 { 604 z.form = zero 605 return z 606 } 607 // x != 0 608 z.mant = z.mant.set(x.abs) 609 fnorm(z.mant) 610 z.setExpAndRound(int64(bits), 0) 611 return z 612 } 613 614 // SetRat sets z to the (possibly rounded) value of x and returns z. 615 // If z's precision is 0, it is changed to the largest of a.BitLen(), 616 // b.BitLen(), or 64; with x = a/b. 617 func (z *Float) SetRat(x *Rat) *Float { 618 if x.IsInt() { 619 return z.SetInt(x.Num()) 620 } 621 var a, b Float 622 a.SetInt(x.Num()) 623 b.SetInt(x.Denom()) 624 if z.prec == 0 { 625 z.prec = umax32(a.prec, b.prec) 626 } 627 return z.Quo(&a, &b) 628 } 629 630 // SetInf sets z to the infinite Float -Inf if signbit is 631 // set, or +Inf if signbit is not set, and returns z. The 632 // precision of z is unchanged and the result is always 633 // Exact. 634 func (z *Float) SetInf(signbit bool) *Float { 635 z.acc = Exact 636 z.form = inf 637 z.neg = signbit 638 return z 639 } 640 641 // Set sets z to the (possibly rounded) value of x and returns z. 642 // If z's precision is 0, it is changed to the precision of x 643 // before setting z (and rounding will have no effect). 644 // Rounding is performed according to z's precision and rounding 645 // mode; and z's accuracy reports the result error relative to the 646 // exact (not rounded) result. 647 func (z *Float) Set(x *Float) *Float { 648 if debugFloat { 649 x.validate() 650 } 651 z.acc = Exact 652 if z != x { 653 z.form = x.form 654 z.neg = x.neg 655 if x.form == finite { 656 z.exp = x.exp 657 z.mant = z.mant.set(x.mant) 658 } 659 if z.prec == 0 { 660 z.prec = x.prec 661 } else if z.prec < x.prec { 662 z.round(0) 663 } 664 } 665 return z 666 } 667 668 // Copy sets z to x, with the same precision, rounding mode, and 669 // accuracy as x, and returns z. x is not changed even if z and 670 // x are the same. 671 func (z *Float) Copy(x *Float) *Float { 672 if debugFloat { 673 x.validate() 674 } 675 if z != x { 676 z.prec = x.prec 677 z.mode = x.mode 678 z.acc = x.acc 679 z.form = x.form 680 z.neg = x.neg 681 if z.form == finite { 682 z.mant = z.mant.set(x.mant) 683 z.exp = x.exp 684 } 685 } 686 return z 687 } 688 689 // msb32 returns the 32 most significant bits of x. 690 func msb32(x nat) uint32 { 691 i := len(x) - 1 692 if i < 0 { 693 return 0 694 } 695 if debugFloat && x[i]&(1<<(_W-1)) == 0 { 696 panic("x not normalized") 697 } 698 switch _W { 699 case 32: 700 return uint32(x[i]) 701 case 64: 702 return uint32(x[i] >> 32) 703 } 704 panic("unreachable") 705 } 706 707 // msb64 returns the 64 most significant bits of x. 708 func msb64(x nat) uint64 { 709 i := len(x) - 1 710 if i < 0 { 711 return 0 712 } 713 if debugFloat && x[i]&(1<<(_W-1)) == 0 { 714 panic("x not normalized") 715 } 716 switch _W { 717 case 32: 718 v := uint64(x[i]) << 32 719 if i > 0 { 720 v |= uint64(x[i-1]) 721 } 722 return v 723 case 64: 724 return uint64(x[i]) 725 } 726 panic("unreachable") 727 } 728 729 // Uint64 returns the unsigned integer resulting from truncating x 730 // towards zero. If 0 <= x <= math.MaxUint64, the result is Exact 731 // if x is an integer and Below otherwise. 732 // The result is (0, Above) for x < 0, and (math.MaxUint64, Below) 733 // for x > math.MaxUint64. 734 func (x *Float) Uint64() (uint64, Accuracy) { 735 if debugFloat { 736 x.validate() 737 } 738 739 switch x.form { 740 case finite: 741 if x.neg { 742 return 0, Above 743 } 744 // 0 < x < +Inf 745 if x.exp <= 0 { 746 // 0 < x < 1 747 return 0, Below 748 } 749 // 1 <= x < Inf 750 if x.exp <= 64 { 751 // u = trunc(x) fits into a uint64 752 u := msb64(x.mant) >> (64 - uint32(x.exp)) 753 if x.MinPrec() <= 64 { 754 return u, Exact 755 } 756 return u, Below // x truncated 757 } 758 // x too large 759 return math.MaxUint64, Below 760 761 case zero: 762 return 0, Exact 763 764 case inf: 765 if x.neg { 766 return 0, Above 767 } 768 return math.MaxUint64, Below 769 } 770 771 panic("unreachable") 772 } 773 774 // Int64 returns the integer resulting from truncating x towards zero. 775 // If math.MinInt64 <= x <= math.MaxInt64, the result is Exact if x is 776 // an integer, and Above (x < 0) or Below (x > 0) otherwise. 777 // The result is (math.MinInt64, Above) for x < math.MinInt64, 778 // and (math.MaxInt64, Below) for x > math.MaxInt64. 779 func (x *Float) Int64() (int64, Accuracy) { 780 if debugFloat { 781 x.validate() 782 } 783 784 switch x.form { 785 case finite: 786 // 0 < |x| < +Inf 787 acc := makeAcc(x.neg) 788 if x.exp <= 0 { 789 // 0 < |x| < 1 790 return 0, acc 791 } 792 // x.exp > 0 793 794 // 1 <= |x| < +Inf 795 if x.exp <= 63 { 796 // i = trunc(x) fits into an int64 (excluding math.MinInt64) 797 i := int64(msb64(x.mant) >> (64 - uint32(x.exp))) 798 if x.neg { 799 i = -i 800 } 801 if x.MinPrec() <= uint(x.exp) { 802 return i, Exact 803 } 804 return i, acc // x truncated 805 } 806 if x.neg { 807 // check for special case x == math.MinInt64 (i.e., x == -(0.5 << 64)) 808 if x.exp == 64 && x.MinPrec() == 1 { 809 acc = Exact 810 } 811 return math.MinInt64, acc 812 } 813 // x too large 814 return math.MaxInt64, Below 815 816 case zero: 817 return 0, Exact 818 819 case inf: 820 if x.neg { 821 return math.MinInt64, Above 822 } 823 return math.MaxInt64, Below 824 } 825 826 panic("unreachable") 827 } 828 829 // Float32 returns the float32 value nearest to x. If x is too small to be 830 // represented by a float32 (|x| < math.SmallestNonzeroFloat32), the result 831 // is (0, Below) or (-0, Above), respectively, depending on the sign of x. 832 // If x is too large to be represented by a float32 (|x| > math.MaxFloat32), 833 // the result is (+Inf, Above) or (-Inf, Below), depending on the sign of x. 834 func (x *Float) Float32() (float32, Accuracy) { 835 if debugFloat { 836 x.validate() 837 } 838 839 switch x.form { 840 case finite: 841 // 0 < |x| < +Inf 842 843 const ( 844 fbits = 32 // float size 845 mbits = 23 // mantissa size (excluding implicit msb) 846 ebits = fbits - mbits - 1 // 8 exponent size 847 bias = 1<<(ebits-1) - 1 // 127 exponent bias 848 dmin = 1 - bias - mbits // -149 smallest unbiased exponent (denormal) 849 emin = 1 - bias // -126 smallest unbiased exponent (normal) 850 emax = bias // 127 largest unbiased exponent (normal) 851 ) 852 853 // Float mantissa m is 0.5 <= m < 1.0; compute exponent e for float32 mantissa. 854 e := x.exp - 1 // exponent for normal mantissa m with 1.0 <= m < 2.0 855 856 // Compute precision p for float32 mantissa. 857 // If the exponent is too small, we have a denormal number before 858 // rounding and fewer than p mantissa bits of precision available 859 // (the exponent remains fixed but the mantissa gets shifted right). 860 p := mbits + 1 // precision of normal float 861 if e < emin { 862 // recompute precision 863 p = mbits + 1 - emin + int(e) 864 // If p == 0, the mantissa of x is shifted so much to the right 865 // that its msb falls immediately to the right of the float32 866 // mantissa space. In other words, if the smallest denormal is 867 // considered "1.0", for p == 0, the mantissa value m is >= 0.5. 868 // If m > 0.5, it is rounded up to 1.0; i.e., the smallest denormal. 869 // If m == 0.5, it is rounded down to even, i.e., 0.0. 870 // If p < 0, the mantissa value m is <= "0.25" which is never rounded up. 871 if p < 0 /* m <= 0.25 */ || p == 0 && x.mant.sticky(uint(len(x.mant))*_W-1) == 0 /* m == 0.5 */ { 872 // underflow to ±0 873 if x.neg { 874 var z float32 875 return -z, Above 876 } 877 return 0.0, Below 878 } 879 // otherwise, round up 880 // We handle p == 0 explicitly because it's easy and because 881 // Float.round doesn't support rounding to 0 bits of precision. 882 if p == 0 { 883 if x.neg { 884 return -math.SmallestNonzeroFloat32, Below 885 } 886 return math.SmallestNonzeroFloat32, Above 887 } 888 } 889 // p > 0 890 891 // round 892 var r Float 893 r.prec = uint32(p) 894 r.Set(x) 895 e = r.exp - 1 896 897 // Rounding may have caused r to overflow to ±Inf 898 // (rounding never causes underflows to 0). 899 // If the exponent is too large, also overflow to ±Inf. 900 if r.form == inf || e > emax { 901 // overflow 902 if x.neg { 903 return float32(math.Inf(-1)), Below 904 } 905 return float32(math.Inf(+1)), Above 906 } 907 // e <= emax 908 909 // Determine sign, biased exponent, and mantissa. 910 var sign, bexp, mant uint32 911 if x.neg { 912 sign = 1 << (fbits - 1) 913 } 914 915 // Rounding may have caused a denormal number to 916 // become normal. Check again. 917 if e < emin { 918 // denormal number: recompute precision 919 // Since rounding may have at best increased precision 920 // and we have eliminated p <= 0 early, we know p > 0. 921 // bexp == 0 for denormals 922 p = mbits + 1 - emin + int(e) 923 mant = msb32(r.mant) >> uint(fbits-p) 924 } else { 925 // normal number: emin <= e <= emax 926 bexp = uint32(e+bias) << mbits 927 mant = msb32(r.mant) >> ebits & (1<<mbits - 1) // cut off msb (implicit 1 bit) 928 } 929 930 return math.Float32frombits(sign | bexp | mant), r.acc 931 932 case zero: 933 if x.neg { 934 var z float32 935 return -z, Exact 936 } 937 return 0.0, Exact 938 939 case inf: 940 if x.neg { 941 return float32(math.Inf(-1)), Exact 942 } 943 return float32(math.Inf(+1)), Exact 944 } 945 946 panic("unreachable") 947 } 948 949 // Float64 returns the float64 value nearest to x. If x is too small to be 950 // represented by a float64 (|x| < math.SmallestNonzeroFloat64), the result 951 // is (0, Below) or (-0, Above), respectively, depending on the sign of x. 952 // If x is too large to be represented by a float64 (|x| > math.MaxFloat64), 953 // the result is (+Inf, Above) or (-Inf, Below), depending on the sign of x. 954 func (x *Float) Float64() (float64, Accuracy) { 955 if debugFloat { 956 x.validate() 957 } 958 959 switch x.form { 960 case finite: 961 // 0 < |x| < +Inf 962 963 const ( 964 fbits = 64 // float size 965 mbits = 52 // mantissa size (excluding implicit msb) 966 ebits = fbits - mbits - 1 // 11 exponent size 967 bias = 1<<(ebits-1) - 1 // 1023 exponent bias 968 dmin = 1 - bias - mbits // -1074 smallest unbiased exponent (denormal) 969 emin = 1 - bias // -1022 smallest unbiased exponent (normal) 970 emax = bias // 1023 largest unbiased exponent (normal) 971 ) 972 973 // Float mantissa m is 0.5 <= m < 1.0; compute exponent e for float64 mantissa. 974 e := x.exp - 1 // exponent for normal mantissa m with 1.0 <= m < 2.0 975 976 // Compute precision p for float64 mantissa. 977 // If the exponent is too small, we have a denormal number before 978 // rounding and fewer than p mantissa bits of precision available 979 // (the exponent remains fixed but the mantissa gets shifted right). 980 p := mbits + 1 // precision of normal float 981 if e < emin { 982 // recompute precision 983 p = mbits + 1 - emin + int(e) 984 // If p == 0, the mantissa of x is shifted so much to the right 985 // that its msb falls immediately to the right of the float64 986 // mantissa space. In other words, if the smallest denormal is 987 // considered "1.0", for p == 0, the mantissa value m is >= 0.5. 988 // If m > 0.5, it is rounded up to 1.0; i.e., the smallest denormal. 989 // If m == 0.5, it is rounded down to even, i.e., 0.0. 990 // If p < 0, the mantissa value m is <= "0.25" which is never rounded up. 991 if p < 0 /* m <= 0.25 */ || p == 0 && x.mant.sticky(uint(len(x.mant))*_W-1) == 0 /* m == 0.5 */ { 992 // underflow to ±0 993 if x.neg { 994 var z float64 995 return -z, Above 996 } 997 return 0.0, Below 998 } 999 // otherwise, round up 1000 // We handle p == 0 explicitly because it's easy and because 1001 // Float.round doesn't support rounding to 0 bits of precision. 1002 if p == 0 { 1003 if x.neg { 1004 return -math.SmallestNonzeroFloat64, Below 1005 } 1006 return math.SmallestNonzeroFloat64, Above 1007 } 1008 } 1009 // p > 0 1010 1011 // round 1012 var r Float 1013 r.prec = uint32(p) 1014 r.Set(x) 1015 e = r.exp - 1 1016 1017 // Rounding may have caused r to overflow to ±Inf 1018 // (rounding never causes underflows to 0). 1019 // If the exponent is too large, also overflow to ±Inf. 1020 if r.form == inf || e > emax { 1021 // overflow 1022 if x.neg { 1023 return math.Inf(-1), Below 1024 } 1025 return math.Inf(+1), Above 1026 } 1027 // e <= emax 1028 1029 // Determine sign, biased exponent, and mantissa. 1030 var sign, bexp, mant uint64 1031 if x.neg { 1032 sign = 1 << (fbits - 1) 1033 } 1034 1035 // Rounding may have caused a denormal number to 1036 // become normal. Check again. 1037 if e < emin { 1038 // denormal number: recompute precision 1039 // Since rounding may have at best increased precision 1040 // and we have eliminated p <= 0 early, we know p > 0. 1041 // bexp == 0 for denormals 1042 p = mbits + 1 - emin + int(e) 1043 mant = msb64(r.mant) >> uint(fbits-p) 1044 } else { 1045 // normal number: emin <= e <= emax 1046 bexp = uint64(e+bias) << mbits 1047 mant = msb64(r.mant) >> ebits & (1<<mbits - 1) // cut off msb (implicit 1 bit) 1048 } 1049 1050 return math.Float64frombits(sign | bexp | mant), r.acc 1051 1052 case zero: 1053 if x.neg { 1054 var z float64 1055 return -z, Exact 1056 } 1057 return 0.0, Exact 1058 1059 case inf: 1060 if x.neg { 1061 return math.Inf(-1), Exact 1062 } 1063 return math.Inf(+1), Exact 1064 } 1065 1066 panic("unreachable") 1067 } 1068 1069 // Int returns the result of truncating x towards zero; 1070 // or nil if x is an infinity. 1071 // The result is Exact if x.IsInt(); otherwise it is Below 1072 // for x > 0, and Above for x < 0. 1073 // If a non-nil *Int argument z is provided, Int stores 1074 // the result in z instead of allocating a new Int. 1075 func (x *Float) Int(z *Int) (*Int, Accuracy) { 1076 if debugFloat { 1077 x.validate() 1078 } 1079 1080 if z == nil && x.form <= finite { 1081 z = new(Int) 1082 } 1083 1084 switch x.form { 1085 case finite: 1086 // 0 < |x| < +Inf 1087 acc := makeAcc(x.neg) 1088 if x.exp <= 0 { 1089 // 0 < |x| < 1 1090 return z.SetInt64(0), acc 1091 } 1092 // x.exp > 0 1093 1094 // 1 <= |x| < +Inf 1095 // determine minimum required precision for x 1096 allBits := uint(len(x.mant)) * _W 1097 exp := uint(x.exp) 1098 if x.MinPrec() <= exp { 1099 acc = Exact 1100 } 1101 // shift mantissa as needed 1102 if z == nil { 1103 z = new(Int) 1104 } 1105 z.neg = x.neg 1106 switch { 1107 case exp > allBits: 1108 z.abs = z.abs.shl(x.mant, exp-allBits) 1109 default: 1110 z.abs = z.abs.set(x.mant) 1111 case exp < allBits: 1112 z.abs = z.abs.shr(x.mant, allBits-exp) 1113 } 1114 return z, acc 1115 1116 case zero: 1117 return z.SetInt64(0), Exact 1118 1119 case inf: 1120 return nil, makeAcc(x.neg) 1121 } 1122 1123 panic("unreachable") 1124 } 1125 1126 // Rat returns the rational number corresponding to x; 1127 // or nil if x is an infinity. 1128 // The result is Exact if x is not an Inf. 1129 // If a non-nil *Rat argument z is provided, Rat stores 1130 // the result in z instead of allocating a new Rat. 1131 func (x *Float) Rat(z *Rat) (*Rat, Accuracy) { 1132 if debugFloat { 1133 x.validate() 1134 } 1135 1136 if z == nil && x.form <= finite { 1137 z = new(Rat) 1138 } 1139 1140 switch x.form { 1141 case finite: 1142 // 0 < |x| < +Inf 1143 allBits := int32(len(x.mant)) * _W 1144 // build up numerator and denominator 1145 z.a.neg = x.neg 1146 switch { 1147 case x.exp > allBits: 1148 z.a.abs = z.a.abs.shl(x.mant, uint(x.exp-allBits)) 1149 z.b.abs = z.b.abs[:0] // == 1 (see Rat) 1150 // z already in normal form 1151 default: 1152 z.a.abs = z.a.abs.set(x.mant) 1153 z.b.abs = z.b.abs[:0] // == 1 (see Rat) 1154 // z already in normal form 1155 case x.exp < allBits: 1156 z.a.abs = z.a.abs.set(x.mant) 1157 t := z.b.abs.setUint64(1) 1158 z.b.abs = t.shl(t, uint(allBits-x.exp)) 1159 z.norm() 1160 } 1161 return z, Exact 1162 1163 case zero: 1164 return z.SetInt64(0), Exact 1165 1166 case inf: 1167 return nil, makeAcc(x.neg) 1168 } 1169 1170 panic("unreachable") 1171 } 1172 1173 // Abs sets z to the (possibly rounded) value |x| (the absolute value of x) 1174 // and returns z. 1175 func (z *Float) Abs(x *Float) *Float { 1176 z.Set(x) 1177 z.neg = false 1178 return z 1179 } 1180 1181 // Neg sets z to the (possibly rounded) value of x with its sign negated, 1182 // and returns z. 1183 func (z *Float) Neg(x *Float) *Float { 1184 z.Set(x) 1185 z.neg = !z.neg 1186 return z 1187 } 1188 1189 func validateBinaryOperands(x, y *Float) { 1190 if !debugFloat { 1191 // avoid performance bugs 1192 panic("validateBinaryOperands called but debugFloat is not set") 1193 } 1194 if len(x.mant) == 0 { 1195 panic("empty mantissa for x") 1196 } 1197 if len(y.mant) == 0 { 1198 panic("empty mantissa for y") 1199 } 1200 } 1201 1202 // z = x + y, ignoring signs of x and y for the addition 1203 // but using the sign of z for rounding the result. 1204 // x and y must have a non-empty mantissa and valid exponent. 1205 func (z *Float) uadd(x, y *Float) { 1206 // Note: This implementation requires 2 shifts most of the 1207 // time. It is also inefficient if exponents or precisions 1208 // differ by wide margins. The following article describes 1209 // an efficient (but much more complicated) implementation 1210 // compatible with the internal representation used here: 1211 // 1212 // Vincent Lefèvre: "The Generic Multiple-Precision Floating- 1213 // Point Addition With Exact Rounding (as in the MPFR Library)" 1214 // http://www.vinc17.net/research/papers/rnc6.pdf 1215 1216 if debugFloat { 1217 validateBinaryOperands(x, y) 1218 } 1219 1220 // compute exponents ex, ey for mantissa with "binary point" 1221 // on the right (mantissa.0) - use int64 to avoid overflow 1222 ex := int64(x.exp) - int64(len(x.mant))*_W 1223 ey := int64(y.exp) - int64(len(y.mant))*_W 1224 1225 al := alias(z.mant, x.mant) || alias(z.mant, y.mant) 1226 1227 // TODO(gri) having a combined add-and-shift primitive 1228 // could make this code significantly faster 1229 switch { 1230 case ex < ey: 1231 if al { 1232 t := nat(nil).shl(y.mant, uint(ey-ex)) 1233 z.mant = z.mant.add(x.mant, t) 1234 } else { 1235 z.mant = z.mant.shl(y.mant, uint(ey-ex)) 1236 z.mant = z.mant.add(x.mant, z.mant) 1237 } 1238 default: 1239 // ex == ey, no shift needed 1240 z.mant = z.mant.add(x.mant, y.mant) 1241 case ex > ey: 1242 if al { 1243 t := nat(nil).shl(x.mant, uint(ex-ey)) 1244 z.mant = z.mant.add(t, y.mant) 1245 } else { 1246 z.mant = z.mant.shl(x.mant, uint(ex-ey)) 1247 z.mant = z.mant.add(z.mant, y.mant) 1248 } 1249 ex = ey 1250 } 1251 // len(z.mant) > 0 1252 1253 z.setExpAndRound(ex+int64(len(z.mant))*_W-fnorm(z.mant), 0) 1254 } 1255 1256 // z = x - y for |x| > |y|, ignoring signs of x and y for the subtraction 1257 // but using the sign of z for rounding the result. 1258 // x and y must have a non-empty mantissa and valid exponent. 1259 func (z *Float) usub(x, y *Float) { 1260 // This code is symmetric to uadd. 1261 // We have not factored the common code out because 1262 // eventually uadd (and usub) should be optimized 1263 // by special-casing, and the code will diverge. 1264 1265 if debugFloat { 1266 validateBinaryOperands(x, y) 1267 } 1268 1269 ex := int64(x.exp) - int64(len(x.mant))*_W 1270 ey := int64(y.exp) - int64(len(y.mant))*_W 1271 1272 al := alias(z.mant, x.mant) || alias(z.mant, y.mant) 1273 1274 switch { 1275 case ex < ey: 1276 if al { 1277 t := nat(nil).shl(y.mant, uint(ey-ex)) 1278 z.mant = t.sub(x.mant, t) 1279 } else { 1280 z.mant = z.mant.shl(y.mant, uint(ey-ex)) 1281 z.mant = z.mant.sub(x.mant, z.mant) 1282 } 1283 default: 1284 // ex == ey, no shift needed 1285 z.mant = z.mant.sub(x.mant, y.mant) 1286 case ex > ey: 1287 if al { 1288 t := nat(nil).shl(x.mant, uint(ex-ey)) 1289 z.mant = t.sub(t, y.mant) 1290 } else { 1291 z.mant = z.mant.shl(x.mant, uint(ex-ey)) 1292 z.mant = z.mant.sub(z.mant, y.mant) 1293 } 1294 ex = ey 1295 } 1296 1297 // operands may have canceled each other out 1298 if len(z.mant) == 0 { 1299 z.acc = Exact 1300 z.form = zero 1301 z.neg = false 1302 return 1303 } 1304 // len(z.mant) > 0 1305 1306 z.setExpAndRound(ex+int64(len(z.mant))*_W-fnorm(z.mant), 0) 1307 } 1308 1309 // z = x * y, ignoring signs of x and y for the multiplication 1310 // but using the sign of z for rounding the result. 1311 // x and y must have a non-empty mantissa and valid exponent. 1312 func (z *Float) umul(x, y *Float) { 1313 if debugFloat { 1314 validateBinaryOperands(x, y) 1315 } 1316 1317 // Note: This is doing too much work if the precision 1318 // of z is less than the sum of the precisions of x 1319 // and y which is often the case (e.g., if all floats 1320 // have the same precision). 1321 // TODO(gri) Optimize this for the common case. 1322 1323 e := int64(x.exp) + int64(y.exp) 1324 if x == y { 1325 z.mant = z.mant.sqr(x.mant) 1326 } else { 1327 z.mant = z.mant.mul(x.mant, y.mant) 1328 } 1329 z.setExpAndRound(e-fnorm(z.mant), 0) 1330 } 1331 1332 // z = x / y, ignoring signs of x and y for the division 1333 // but using the sign of z for rounding the result. 1334 // x and y must have a non-empty mantissa and valid exponent. 1335 func (z *Float) uquo(x, y *Float) { 1336 if debugFloat { 1337 validateBinaryOperands(x, y) 1338 } 1339 1340 // mantissa length in words for desired result precision + 1 1341 // (at least one extra bit so we get the rounding bit after 1342 // the division) 1343 n := int(z.prec/_W) + 1 1344 1345 // compute adjusted x.mant such that we get enough result precision 1346 xadj := x.mant 1347 if d := n - len(x.mant) + len(y.mant); d > 0 { 1348 // d extra words needed => add d "0 digits" to x 1349 xadj = make(nat, len(x.mant)+d) 1350 copy(xadj[d:], x.mant) 1351 } 1352 // TODO(gri): If we have too many digits (d < 0), we should be able 1353 // to shorten x for faster division. But we must be extra careful 1354 // with rounding in that case. 1355 1356 // Compute d before division since there may be aliasing of x.mant 1357 // (via xadj) or y.mant with z.mant. 1358 d := len(xadj) - len(y.mant) 1359 1360 // divide 1361 var r nat 1362 z.mant, r = z.mant.div(nil, xadj, y.mant) 1363 e := int64(x.exp) - int64(y.exp) - int64(d-len(z.mant))*_W 1364 1365 // The result is long enough to include (at least) the rounding bit. 1366 // If there's a non-zero remainder, the corresponding fractional part 1367 // (if it were computed), would have a non-zero sticky bit (if it were 1368 // zero, it couldn't have a non-zero remainder). 1369 var sbit uint 1370 if len(r) > 0 { 1371 sbit = 1 1372 } 1373 1374 z.setExpAndRound(e-fnorm(z.mant), sbit) 1375 } 1376 1377 // ucmp returns -1, 0, or +1, depending on whether 1378 // |x| < |y|, |x| == |y|, or |x| > |y|. 1379 // x and y must have a non-empty mantissa and valid exponent. 1380 func (x *Float) ucmp(y *Float) int { 1381 if debugFloat { 1382 validateBinaryOperands(x, y) 1383 } 1384 1385 switch { 1386 case x.exp < y.exp: 1387 return -1 1388 case x.exp > y.exp: 1389 return +1 1390 } 1391 // x.exp == y.exp 1392 1393 // compare mantissas 1394 i := len(x.mant) 1395 j := len(y.mant) 1396 for i > 0 || j > 0 { 1397 var xm, ym Word 1398 if i > 0 { 1399 i-- 1400 xm = x.mant[i] 1401 } 1402 if j > 0 { 1403 j-- 1404 ym = y.mant[j] 1405 } 1406 switch { 1407 case xm < ym: 1408 return -1 1409 case xm > ym: 1410 return +1 1411 } 1412 } 1413 1414 return 0 1415 } 1416 1417 // Handling of sign bit as defined by IEEE 754-2008, section 6.3: 1418 // 1419 // When neither the inputs nor result are NaN, the sign of a product or 1420 // quotient is the exclusive OR of the operands’ signs; the sign of a sum, 1421 // or of a difference x−y regarded as a sum x+(−y), differs from at most 1422 // one of the addends’ signs; and the sign of the result of conversions, 1423 // the quantize operation, the roundToIntegral operations, and the 1424 // roundToIntegralExact (see 5.3.1) is the sign of the first or only operand. 1425 // These rules shall apply even when operands or results are zero or infinite. 1426 // 1427 // When the sum of two operands with opposite signs (or the difference of 1428 // two operands with like signs) is exactly zero, the sign of that sum (or 1429 // difference) shall be +0 in all rounding-direction attributes except 1430 // roundTowardNegative; under that attribute, the sign of an exact zero 1431 // sum (or difference) shall be −0. However, x+x = x−(−x) retains the same 1432 // sign as x even when x is zero. 1433 // 1434 // See also: https://play.golang.org/p/RtH3UCt5IH 1435 1436 // Add sets z to the rounded sum x+y and returns z. If z's precision is 0, 1437 // it is changed to the larger of x's or y's precision before the operation. 1438 // Rounding is performed according to z's precision and rounding mode; and 1439 // z's accuracy reports the result error relative to the exact (not rounded) 1440 // result. Add panics with ErrNaN if x and y are infinities with opposite 1441 // signs. The value of z is undefined in that case. 1442 func (z *Float) Add(x, y *Float) *Float { 1443 if debugFloat { 1444 x.validate() 1445 y.validate() 1446 } 1447 1448 if z.prec == 0 { 1449 z.prec = umax32(x.prec, y.prec) 1450 } 1451 1452 if x.form == finite && y.form == finite { 1453 // x + y (common case) 1454 1455 // Below we set z.neg = x.neg, and when z aliases y this will 1456 // change the y operand's sign. This is fine, because if an 1457 // operand aliases the receiver it'll be overwritten, but we still 1458 // want the original x.neg and y.neg values when we evaluate 1459 // x.neg != y.neg, so we need to save y.neg before setting z.neg. 1460 yneg := y.neg 1461 1462 z.neg = x.neg 1463 if x.neg == yneg { 1464 // x + y == x + y 1465 // (-x) + (-y) == -(x + y) 1466 z.uadd(x, y) 1467 } else { 1468 // x + (-y) == x - y == -(y - x) 1469 // (-x) + y == y - x == -(x - y) 1470 if x.ucmp(y) > 0 { 1471 z.usub(x, y) 1472 } else { 1473 z.neg = !z.neg 1474 z.usub(y, x) 1475 } 1476 } 1477 if z.form == zero && z.mode == ToNegativeInf && z.acc == Exact { 1478 z.neg = true 1479 } 1480 return z 1481 } 1482 1483 if x.form == inf && y.form == inf && x.neg != y.neg { 1484 // +Inf + -Inf 1485 // -Inf + +Inf 1486 // value of z is undefined but make sure it's valid 1487 z.acc = Exact 1488 z.form = zero 1489 z.neg = false 1490 panic(ErrNaN{"addition of infinities with opposite signs"}) 1491 } 1492 1493 if x.form == zero && y.form == zero { 1494 // ±0 + ±0 1495 z.acc = Exact 1496 z.form = zero 1497 z.neg = x.neg && y.neg // -0 + -0 == -0 1498 return z 1499 } 1500 1501 if x.form == inf || y.form == zero { 1502 // ±Inf + y 1503 // x + ±0 1504 return z.Set(x) 1505 } 1506 1507 // ±0 + y 1508 // x + ±Inf 1509 return z.Set(y) 1510 } 1511 1512 // Sub sets z to the rounded difference x-y and returns z. 1513 // Precision, rounding, and accuracy reporting are as for Add. 1514 // Sub panics with ErrNaN if x and y are infinities with equal 1515 // signs. The value of z is undefined in that case. 1516 func (z *Float) Sub(x, y *Float) *Float { 1517 if debugFloat { 1518 x.validate() 1519 y.validate() 1520 } 1521 1522 if z.prec == 0 { 1523 z.prec = umax32(x.prec, y.prec) 1524 } 1525 1526 if x.form == finite && y.form == finite { 1527 // x - y (common case) 1528 yneg := y.neg 1529 z.neg = x.neg 1530 if x.neg != yneg { 1531 // x - (-y) == x + y 1532 // (-x) - y == -(x + y) 1533 z.uadd(x, y) 1534 } else { 1535 // x - y == x - y == -(y - x) 1536 // (-x) - (-y) == y - x == -(x - y) 1537 if x.ucmp(y) > 0 { 1538 z.usub(x, y) 1539 } else { 1540 z.neg = !z.neg 1541 z.usub(y, x) 1542 } 1543 } 1544 if z.form == zero && z.mode == ToNegativeInf && z.acc == Exact { 1545 z.neg = true 1546 } 1547 return z 1548 } 1549 1550 if x.form == inf && y.form == inf && x.neg == y.neg { 1551 // +Inf - +Inf 1552 // -Inf - -Inf 1553 // value of z is undefined but make sure it's valid 1554 z.acc = Exact 1555 z.form = zero 1556 z.neg = false 1557 panic(ErrNaN{"subtraction of infinities with equal signs"}) 1558 } 1559 1560 if x.form == zero && y.form == zero { 1561 // ±0 - ±0 1562 z.acc = Exact 1563 z.form = zero 1564 z.neg = x.neg && !y.neg // -0 - +0 == -0 1565 return z 1566 } 1567 1568 if x.form == inf || y.form == zero { 1569 // ±Inf - y 1570 // x - ±0 1571 return z.Set(x) 1572 } 1573 1574 // ±0 - y 1575 // x - ±Inf 1576 return z.Neg(y) 1577 } 1578 1579 // Mul sets z to the rounded product x*y and returns z. 1580 // Precision, rounding, and accuracy reporting are as for Add. 1581 // Mul panics with ErrNaN if one operand is zero and the other 1582 // operand an infinity. The value of z is undefined in that case. 1583 func (z *Float) Mul(x, y *Float) *Float { 1584 if debugFloat { 1585 x.validate() 1586 y.validate() 1587 } 1588 1589 if z.prec == 0 { 1590 z.prec = umax32(x.prec, y.prec) 1591 } 1592 1593 z.neg = x.neg != y.neg 1594 1595 if x.form == finite && y.form == finite { 1596 // x * y (common case) 1597 z.umul(x, y) 1598 return z 1599 } 1600 1601 z.acc = Exact 1602 if x.form == zero && y.form == inf || x.form == inf && y.form == zero { 1603 // ±0 * ±Inf 1604 // ±Inf * ±0 1605 // value of z is undefined but make sure it's valid 1606 z.form = zero 1607 z.neg = false 1608 panic(ErrNaN{"multiplication of zero with infinity"}) 1609 } 1610 1611 if x.form == inf || y.form == inf { 1612 // ±Inf * y 1613 // x * ±Inf 1614 z.form = inf 1615 return z 1616 } 1617 1618 // ±0 * y 1619 // x * ±0 1620 z.form = zero 1621 return z 1622 } 1623 1624 // Quo sets z to the rounded quotient x/y and returns z. 1625 // Precision, rounding, and accuracy reporting are as for Add. 1626 // Quo panics with ErrNaN if both operands are zero or infinities. 1627 // The value of z is undefined in that case. 1628 func (z *Float) Quo(x, y *Float) *Float { 1629 if debugFloat { 1630 x.validate() 1631 y.validate() 1632 } 1633 1634 if z.prec == 0 { 1635 z.prec = umax32(x.prec, y.prec) 1636 } 1637 1638 z.neg = x.neg != y.neg 1639 1640 if x.form == finite && y.form == finite { 1641 // x / y (common case) 1642 z.uquo(x, y) 1643 return z 1644 } 1645 1646 z.acc = Exact 1647 if x.form == zero && y.form == zero || x.form == inf && y.form == inf { 1648 // ±0 / ±0 1649 // ±Inf / ±Inf 1650 // value of z is undefined but make sure it's valid 1651 z.form = zero 1652 z.neg = false 1653 panic(ErrNaN{"division of zero by zero or infinity by infinity"}) 1654 } 1655 1656 if x.form == zero || y.form == inf { 1657 // ±0 / y 1658 // x / ±Inf 1659 z.form = zero 1660 return z 1661 } 1662 1663 // x / ±0 1664 // ±Inf / y 1665 z.form = inf 1666 return z 1667 } 1668 1669 // Cmp compares x and y and returns: 1670 // 1671 // -1 if x < y 1672 // 0 if x == y (incl. -0 == 0, -Inf == -Inf, and +Inf == +Inf) 1673 // +1 if x > y 1674 func (x *Float) Cmp(y *Float) int { 1675 if debugFloat { 1676 x.validate() 1677 y.validate() 1678 } 1679 1680 mx := x.ord() 1681 my := y.ord() 1682 switch { 1683 case mx < my: 1684 return -1 1685 case mx > my: 1686 return +1 1687 } 1688 // mx == my 1689 1690 // only if |mx| == 1 we have to compare the mantissae 1691 switch mx { 1692 case -1: 1693 return y.ucmp(x) 1694 case +1: 1695 return x.ucmp(y) 1696 } 1697 1698 return 0 1699 } 1700 1701 // ord classifies x and returns: 1702 // 1703 // -2 if -Inf == x 1704 // -1 if -Inf < x < 0 1705 // 0 if x == 0 (signed or unsigned) 1706 // +1 if 0 < x < +Inf 1707 // +2 if x == +Inf 1708 func (x *Float) ord() int { 1709 var m int 1710 switch x.form { 1711 case finite: 1712 m = 1 1713 case zero: 1714 return 0 1715 case inf: 1716 m = 2 1717 } 1718 if x.neg { 1719 m = -m 1720 } 1721 return m 1722 } 1723 1724 func umax32(x, y uint32) uint32 { 1725 if x > y { 1726 return x 1727 } 1728 return y 1729 } 1730