src/runtime/softfloat64.go

   1 // Copyright 2010 The Go Authors. All rights reserved.
   2 // Use of this source code is governed by a BSD-style
   3 // license that can be found in the LICENSE file.
   4
   5 // Software IEEE754 64-bit floating point.
   6 // Only referred to (and thus linked in) by softfloat targets
   7 // and by tests in this directory.
   8
   9 package runtime
  10
  11 const (
  12         mantbits64 uint = 52
  13         expbits64  uint = 11
  14         bias64          = -1<<(expbits64-1) + 1
  15
  16         nan64 uint64 = (1<<expbits64-1)<<mantbits64 + 1<<(mantbits64-1) // quiet NaN, 0 payload
  17         inf64 uint64 = (1<<expbits64 - 1) << mantbits64
  18         neg64 uint64 = 1 << (expbits64 + mantbits64)
  19
  20         mantbits32 uint = 23
  21         expbits32  uint = 8
  22         bias32          = -1<<(expbits32-1) + 1
  23
  24         nan32 uint32 = (1<<expbits32-1)<<mantbits32 + 1<<(mantbits32-1) // quiet NaN, 0 payload
  25         inf32 uint32 = (1<<expbits32 - 1) << mantbits32
  26         neg32 uint32 = 1 << (expbits32 + mantbits32)
  27 )
  28
  29 func funpack64(f uint64) (sign, mant uint64, exp int, inf, nan bool) {
  30         sign = f & (1 << (mantbits64 + expbits64))
  31         mant = f & (1<<mantbits64 - 1)
  32         exp = int(f>>mantbits64) & (1<<expbits64 - 1)
  33
  34         switch exp {
  35         case 1<<expbits64 - 1:
  36                 if mant != 0 {
  37                         nan = true
  38                         return
  39                 }
  40                 inf = true
  41                 return
  42
  43         case 0:
  44                 // denormalized
  45                 if mant != 0 {
  46                         exp += bias64 + 1
  47                         for mant < 1<<mantbits64 {
  48                                 mant <<= 1
  49                                 exp--
  50                         }
  51                 }
  52
  53         default:
  54                 // add implicit top bit
  55                 mant |= 1 << mantbits64
  56                 exp += bias64
  57         }
  58         return
  59 }
  60
  61 func funpack32(f uint32) (sign, mant uint32, exp int, inf, nan bool) {
  62         sign = f & (1 << (mantbits32 + expbits32))
  63         mant = f & (1<<mantbits32 - 1)
  64         exp = int(f>>mantbits32) & (1<<expbits32 - 1)
  65
  66         switch exp {
  67         case 1<<expbits32 - 1:
  68                 if mant != 0 {
  69                         nan = true
  70                         return
  71                 }
  72                 inf = true
  73                 return
  74
  75         case 0:
  76                 // denormalized
  77                 if mant != 0 {
  78                         exp += bias32 + 1
  79                         for mant < 1<<mantbits32 {
  80                                 mant <<= 1
  81                                 exp--
  82                         }
  83                 }
  84
  85         default:
  86                 // add implicit top bit
  87                 mant |= 1 << mantbits32
  88                 exp += bias32
  89         }
  90         return
  91 }
  92
  93 func fpack64(sign, mant uint64, exp int, trunc uint64) uint64 {
  94         mant0, exp0, trunc0 := mant, exp, trunc
  95         if mant == 0 {
  96                 return sign
  97         }
  98         for mant < 1<<mantbits64 {
  99                 mant <<= 1
 100                 exp--
 101         }
 102         for mant >= 4<<mantbits64 {
 103                 trunc |= mant & 1
 104                 mant >>= 1
 105                 exp++
 106         }
 107         if mant >= 2<<mantbits64 {
 108                 if mant&1 != 0 && (trunc != 0 || mant&2 != 0) {
 109                         mant++
 110                         if mant >= 4<<mantbits64 {
 111                                 mant >>= 1
 112                                 exp++
 113                         }
 114                 }
 115                 mant >>= 1
 116                 exp++
 117         }
 118         if exp >= 1<<expbits64-1+bias64 {
 119                 return sign ^ inf64
 120         }
 121         if exp < bias64+1 {
 122                 if exp < bias64-int(mantbits64) {
 123                         return sign | 0
 124                 }
 125                 // repeat expecting denormal
 126                 mant, exp, trunc = mant0, exp0, trunc0
 127                 for exp < bias64 {
 128                         trunc |= mant & 1
 129                         mant >>= 1
 130                         exp++
 131                 }
 132                 if mant&1 != 0 && (trunc != 0 || mant&2 != 0) {
 133                         mant++
 134                 }
 135                 mant >>= 1
 136                 exp++
 137                 if mant < 1<<mantbits64 {
 138                         return sign | mant
 139                 }
 140         }
 141         return sign | uint64(exp-bias64)<<mantbits64 | mant&(1<<mantbits64-1)
 142 }
 143
 144 func fpack32(sign, mant uint32, exp int, trunc uint32) uint32 {
 145         mant0, exp0, trunc0 := mant, exp, trunc
 146         if mant == 0 {
 147                 return sign
 148         }
 149         for mant < 1<<mantbits32 {
 150                 mant <<= 1
 151                 exp--
 152         }
 153         for mant >= 4<<mantbits32 {
 154                 trunc |= mant & 1
 155                 mant >>= 1
 156                 exp++
 157         }
 158         if mant >= 2<<mantbits32 {
 159                 if mant&1 != 0 && (trunc != 0 || mant&2 != 0) {
 160                         mant++
 161                         if mant >= 4<<mantbits32 {
 162                                 mant >>= 1
 163                                 exp++
 164                         }
 165                 }
 166                 mant >>= 1
 167                 exp++
 168         }
 169         if exp >= 1<<expbits32-1+bias32 {
 170                 return sign ^ inf32
 171         }
 172         if exp < bias32+1 {
 173                 if exp < bias32-int(mantbits32) {
 174                         return sign | 0
 175                 }
 176                 // repeat expecting denormal
 177                 mant, exp, trunc = mant0, exp0, trunc0
 178                 for exp < bias32 {
 179                         trunc |= mant & 1
 180                         mant >>= 1
 181                         exp++
 182                 }
 183                 if mant&1 != 0 && (trunc != 0 || mant&2 != 0) {
 184                         mant++
 185                 }
 186                 mant >>= 1
 187                 exp++
 188                 if mant < 1<<mantbits32 {
 189                         return sign | mant
 190                 }
 191         }
 192         return sign | uint32(exp-bias32)<<mantbits32 | mant&(1<<mantbits32-1)
 193 }
 194
 195 func fadd64(f, g uint64) uint64 {
 196         fs, fm, fe, fi, fn := funpack64(f)
 197         gs, gm, ge, gi, gn := funpack64(g)
 198
 199         // Special cases.
 200         switch {
 201         case fn || gn: // NaN + x or x + NaN = NaN
 202                 return nan64
 203
 204         case fi && gi && fs != gs: // +Inf + -Inf or -Inf + +Inf = NaN
 205                 return nan64
 206
 207         case fi: // ±Inf + g = ±Inf
 208                 return f
 209
 210         case gi: // f + ±Inf = ±Inf
 211                 return g
 212
 213         case fm == 0 && gm == 0 && fs != 0 && gs != 0: // -0 + -0 = -0
 214                 return f
 215
 216         case fm == 0: // 0 + g = g but 0 + -0 = +0
 217                 if gm == 0 {
 218                         g ^= gs
 219                 }
 220                 return g
 221
 222         case gm == 0: // f + 0 = f
 223                 return f
 224
 225         }
 226
 227         if fe < ge || fe == ge && fm < gm {
 228                 f, g, fs, fm, fe, gs, gm, ge = g, f, gs, gm, ge, fs, fm, fe
 229         }
 230
 231         shift := uint(fe - ge)
 232         fm <<= 2
 233         gm <<= 2
 234         trunc := gm & (1<<shift - 1)
 235         gm >>= shift
 236         if fs == gs {
 237                 fm += gm
 238         } else {
 239                 fm -= gm
 240                 if trunc != 0 {
 241                         fm--
 242                 }
 243         }
 244         if fm == 0 {
 245                 fs = 0
 246         }
 247         return fpack64(fs, fm, fe-2, trunc)
 248 }
 249
 250 func fsub64(f, g uint64) uint64 {
 251         return fadd64(f, fneg64(g))
 252 }
 253
 254 func fneg64(f uint64) uint64 {
 255         return f ^ (1 << (mantbits64 + expbits64))
 256 }
 257
 258 func fmul64(f, g uint64) uint64 {
 259         fs, fm, fe, fi, fn := funpack64(f)
 260         gs, gm, ge, gi, gn := funpack64(g)
 261
 262         // Special cases.
 263         switch {
 264         case fn || gn: // NaN * g or f * NaN = NaN
 265                 return nan64
 266
 267         case fi && gi: // Inf * Inf = Inf (with sign adjusted)
 268                 return f ^ gs
 269
 270         case fi && gm == 0, fm == 0 && gi: // 0 * Inf = Inf * 0 = NaN
 271                 return nan64
 272
 273         case fm == 0: // 0 * x = 0 (with sign adjusted)
 274                 return f ^ gs
 275
 276         case gm == 0: // x * 0 = 0 (with sign adjusted)
 277                 return g ^ fs
 278         }
 279
 280         // 53-bit * 53-bit = 107- or 108-bit
 281         lo, hi := mullu(fm, gm)
 282         shift := mantbits64 - 1
 283         trunc := lo & (1<<shift - 1)
 284         mant := hi<<(64-shift) | lo>>shift
 285         return fpack64(fs^gs, mant, fe+ge-1, trunc)
 286 }
 287
 288 func fdiv64(f, g uint64) uint64 {
 289         fs, fm, fe, fi, fn := funpack64(f)
 290         gs, gm, ge, gi, gn := funpack64(g)
 291
 292         // Special cases.
 293         switch {
 294         case fn || gn: // NaN / g = f / NaN = NaN
 295                 return nan64
 296
 297         case fi && gi: // ±Inf / ±Inf = NaN
 298                 return nan64
 299
 300         case !fi && !gi && fm == 0 && gm == 0: // 0 / 0 = NaN
 301                 return nan64
 302
 303         case fi, !gi && gm == 0: // Inf / g = f / 0 = Inf
 304                 return fs ^ gs ^ inf64
 305
 306         case gi, fm == 0: // f / Inf = 0 / g = Inf
 307                 return fs ^ gs ^ 0
 308         }
 309         _, _, _, _ = fi, fn, gi, gn
 310
 311         // 53-bit<<54 / 53-bit = 53- or 54-bit.
 312         shift := mantbits64 + 2
 313         q, r := divlu(fm>>(64-shift), fm<<shift, gm)
 314         return fpack64(fs^gs, q, fe-ge-2, r)
 315 }
 316
 317 func f64to32(f uint64) uint32 {
 318         fs, fm, fe, fi, fn := funpack64(f)
 319         if fn {
 320                 return nan32
 321         }
 322         fs32 := uint32(fs >> 32)
 323         if fi {
 324                 return fs32 ^ inf32
 325         }
 326         const d = mantbits64 - mantbits32 - 1
 327         return fpack32(fs32, uint32(fm>>d), fe-1, uint32(fm&(1<<d-1)))
 328 }
 329
 330 func f32to64(f uint32) uint64 {
 331         const d = mantbits64 - mantbits32
 332         fs, fm, fe, fi, fn := funpack32(f)
 333         if fn {
 334                 return nan64
 335         }
 336         fs64 := uint64(fs) << 32
 337         if fi {
 338                 return fs64 ^ inf64
 339         }
 340         return fpack64(fs64, uint64(fm)<<d, fe, 0)
 341 }
 342
 343 func fcmp64(f, g uint64) (cmp int32, isnan bool) {
 344         fs, fm, _, fi, fn := funpack64(f)
 345         gs, gm, _, gi, gn := funpack64(g)
 346
 347         switch {
 348         case fn, gn: // flag NaN
 349                 return 0, true
 350
 351         case !fi && !gi && fm == 0 && gm == 0: // ±0 == ±0
 352                 return 0, false
 353
 354         case fs > gs: // f < 0, g > 0
 355                 return -1, false
 356
 357         case fs < gs: // f > 0, g < 0
 358                 return +1, false
 359
 360         // Same sign, not NaN.
 361         // Can compare encodings directly now.
 362         // Reverse for sign.
 363         case fs == 0 && f < g, fs != 0 && f > g:
 364                 return -1, false
 365
 366         case fs == 0 && f > g, fs != 0 && f < g:
 367                 return +1, false
 368         }
 369
 370         // f == g
 371         return 0, false
 372 }
 373
 374 func f64toint(f uint64) (val int64, ok bool) {
 375         fs, fm, fe, fi, fn := funpack64(f)
 376
 377         switch {
 378         case fi, fn: // NaN
 379                 return 0, false
 380
 381         case fe < -1: // f < 0.5
 382                 return 0, false
 383
 384         case fe > 63: // f >= 2^63
 385                 if fs != 0 && fm == 0 { // f == -2^63
 386                         return -1 << 63, true
 387                 }
 388                 if fs != 0 {
 389                         return 0, false
 390                 }
 391                 return 0, false
 392         }
 393
 394         for fe > int(mantbits64) {
 395                 fe--
 396                 fm <<= 1
 397         }
 398         for fe < int(mantbits64) {
 399                 fe++
 400                 fm >>= 1
 401         }
 402         val = int64(fm)
 403         if fs != 0 {
 404                 val = -val
 405         }
 406         return val, true
 407 }
 408
 409 func fintto64(val int64) (f uint64) {
 410         fs := uint64(val) & (1 << 63)
 411         mant := uint64(val)
 412         if fs != 0 {
 413                 mant = -mant
 414         }
 415         return fpack64(fs, mant, int(mantbits64), 0)
 416 }
 417 func fintto32(val int64) (f uint32) {
 418         fs := uint64(val) & (1 << 63)
 419         mant := uint64(val)
 420         if fs != 0 {
 421                 mant = -mant
 422         }
 423         // Reduce mantissa size until it fits into a uint32.
 424         // Keep track of the bits we throw away, and if any are
 425         // nonzero or them into the lowest bit.
 426         exp := int(mantbits32)
 427         var trunc uint32
 428         for mant >= 1<<32 {
 429                 trunc |= uint32(mant) & 1
 430                 mant >>= 1
 431                 exp++
 432         }
 433
 434         return fpack32(uint32(fs>>32), uint32(mant), exp, trunc)
 435 }
 436
 437 // 64x64 -> 128 multiply.
 438 // adapted from hacker's delight.
 439 func mullu(u, v uint64) (lo, hi uint64) {
 440         const (
 441                 s    = 32
 442                 mask = 1<<s - 1
 443         )
 444         u0 := u & mask
 445         u1 := u >> s
 446         v0 := v & mask
 447         v1 := v >> s
 448         w0 := u0 * v0
 449         t := u1*v0 + w0>>s
 450         w1 := t & mask
 451         w2 := t >> s
 452         w1 += u0 * v1
 453         return u * v, u1*v1 + w2 + w1>>s
 454 }
 455
 456 // 128/64 -> 64 quotient, 64 remainder.
 457 // adapted from hacker's delight
 458 func divlu(u1, u0, v uint64) (q, r uint64) {
 459         const b = 1 << 32
 460
 461         if u1 >= v {
 462                 return 1<<64 - 1, 1<<64 - 1
 463         }
 464
 465         // s = nlz(v); v <<= s
 466         s := uint(0)
 467         for v&(1<<63) == 0 {
 468                 s++
 469                 v <<= 1
 470         }
 471
 472         vn1 := v >> 32
 473         vn0 := v & (1<<32 - 1)
 474         un32 := u1<<s | u0>>(64-s)
 475         un10 := u0 << s
 476         un1 := un10 >> 32
 477         un0 := un10 & (1<<32 - 1)
 478         q1 := un32 / vn1
 479         rhat := un32 - q1*vn1
 480
 481 again1:
 482         if q1 >= b || q1*vn0 > b*rhat+un1 {
 483                 q1--
 484                 rhat += vn1
 485                 if rhat < b {
 486                         goto again1
 487                 }
 488         }
 489
 490         un21 := un32*b + un1 - q1*v
 491         q0 := un21 / vn1
 492         rhat = un21 - q0*vn1
 493
 494 again2:
 495         if q0 >= b || q0*vn0 > b*rhat+un0 {
 496                 q0--
 497                 rhat += vn1
 498                 if rhat < b {
 499                         goto again2
 500                 }
 501         }
 502
 503         return q1*b + q0, (un21*b + un0 - q0*v) >> s
 504 }
 505
 506 func fadd32(x, y uint32) uint32 {
 507         return f64to32(fadd64(f32to64(x), f32to64(y)))
 508 }
 509
 510 func fmul32(x, y uint32) uint32 {
 511         return f64to32(fmul64(f32to64(x), f32to64(y)))
 512 }
 513
 514 func fdiv32(x, y uint32) uint32 {
 515         // TODO: are there double-rounding problems here? See issue 48807.
 516         return f64to32(fdiv64(f32to64(x), f32to64(y)))
 517 }
 518
 519 func feq32(x, y uint32) bool {
 520         cmp, nan := fcmp64(f32to64(x), f32to64(y))
 521         return cmp == 0 && !nan
 522 }
 523
 524 func fgt32(x, y uint32) bool {
 525         cmp, nan := fcmp64(f32to64(x), f32to64(y))
 526         return cmp >= 1 && !nan
 527 }
 528
 529 func fge32(x, y uint32) bool {
 530         cmp, nan := fcmp64(f32to64(x), f32to64(y))
 531         return cmp >= 0 && !nan
 532 }
 533
 534 func feq64(x, y uint64) bool {
 535         cmp, nan := fcmp64(x, y)
 536         return cmp == 0 && !nan
 537 }
 538
 539 func fgt64(x, y uint64) bool {
 540         cmp, nan := fcmp64(x, y)
 541         return cmp >= 1 && !nan
 542 }
 543
 544 func fge64(x, y uint64) bool {
 545         cmp, nan := fcmp64(x, y)
 546         return cmp >= 0 && !nan
 547 }
 548
 549 func fint32to32(x int32) uint32 {
 550         return fintto32(int64(x))
 551 }
 552
 553 func fint32to64(x int32) uint64 {
 554         return fintto64(int64(x))
 555 }
 556
 557 func fint64to32(x int64) uint32 {
 558         return fintto32(x)
 559 }
 560
 561 func fint64to64(x int64) uint64 {
 562         return fintto64(x)
 563 }
 564
 565 func f32toint32(x uint32) int32 {
 566         val, _ := f64toint(f32to64(x))
 567         return int32(val)
 568 }
 569
 570 func f32toint64(x uint32) int64 {
 571         val, _ := f64toint(f32to64(x))
 572         return val
 573 }
 574
 575 func f64toint32(x uint64) int32 {
 576         val, _ := f64toint(x)
 577         return int32(val)
 578 }
 579
 580 func f64toint64(x uint64) int64 {
 581         val, _ := f64toint(x)
 582         return val
 583 }
 584
 585 func f64touint64(x uint64) uint64 {
 586         var m uint64 = 0x43e0000000000000 // float64 1<<63
 587         if fgt64(m, x) {
 588                 return uint64(f64toint64(x))
 589         }
 590         y := fadd64(x, -m)
 591         z := uint64(f64toint64(y))
 592         return z | (1 << 63)
 593 }
 594
 595 func f32touint64(x uint32) uint64 {
 596         var m uint32 = 0x5f000000 // float32 1<<63
 597         if fgt32(m, x) {
 598                 return uint64(f32toint64(x))
 599         }
 600         y := fadd32(x, -m)
 601         z := uint64(f32toint64(y))
 602         return z | (1 << 63)
 603 }
 604
 605 func fuint64to64(x uint64) uint64 {
 606         if int64(x) >= 0 {
 607                 return fint64to64(int64(x))
 608         }
 609         // See ../cmd/compile/internal/ssagen/ssa.go:uint64Tofloat
 610         y := x & 1
 611         z := x >> 1
 612         z = z | y
 613         r := fint64to64(int64(z))
 614         return fadd64(r, r)
 615 }
 616
 617 func fuint64to32(x uint64) uint32 {
 618         if int64(x) >= 0 {
 619                 return fint64to32(int64(x))
 620         }
 621         // See ../cmd/compile/internal/ssagen/ssa.go:uint64Tofloat
 622         y := x & 1
 623         z := x >> 1
 624         z = z | y
 625         r := fint64to32(int64(z))
 626         return fadd32(r, r)
 627 }