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