libgo/go/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 arm port
   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
  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
  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 int, 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
 418 // 64x64 -> 128 multiply.
 419 // adapted from hacker's delight.
 420 func mullu(u, v uint64) (lo, hi uint64) {
 421         const (
 422                 s    = 32
 423                 mask = 1<<s - 1
 424         )
 425         u0 := u & mask
 426         u1 := u >> s
 427         v0 := v & mask
 428         v1 := v >> s
 429         w0 := u0 * v0
 430         t := u1*v0 + w0>>s
 431         w1 := t & mask
 432         w2 := t >> s
 433         w1 += u0 * v1
 434         return u * v, u1*v1 + w2 + w1>>s
 435 }
 436
 437 // 128/64 -> 64 quotient, 64 remainder.
 438 // adapted from hacker's delight
 439 func divlu(u1, u0, v uint64) (q, r uint64) {
 440         const b = 1 << 32
 441
 442         if u1 >= v {
 443                 return 1<<64 - 1, 1<<64 - 1
 444         }
 445
 446         // s = nlz(v); v <<= s
 447         s := uint(0)
 448         for v&(1<<63) == 0 {
 449                 s++
 450                 v <<= 1
 451         }
 452
 453         vn1 := v >> 32
 454         vn0 := v & (1<<32 - 1)
 455         un32 := u1<<s | u0>>(64-s)
 456         un10 := u0 << s
 457         un1 := un10 >> 32
 458         un0 := un10 & (1<<32 - 1)
 459         q1 := un32 / vn1
 460         rhat := un32 - q1*vn1
 461
 462 again1:
 463         if q1 >= b || q1*vn0 > b*rhat+un1 {
 464                 q1--
 465                 rhat += vn1
 466                 if rhat < b {
 467                         goto again1
 468                 }
 469         }
 470
 471         un21 := un32*b + un1 - q1*v
 472         q0 := un21 / vn1
 473         rhat = un21 - q0*vn1
 474
 475 again2:
 476         if q0 >= b || q0*vn0 > b*rhat+un0 {
 477                 q0--
 478                 rhat += vn1
 479                 if rhat < b {
 480                         goto again2
 481                 }
 482         }
 483
 484         return q1*b + q0, (un21*b + un0 - q0*v) >> s
 485 }
 486
 487 // callable from C
 488
 489 func fadd64c(f, g uint64, ret *uint64)            { *ret = fadd64(f, g) }
 490 func fsub64c(f, g uint64, ret *uint64)            { *ret = fsub64(f, g) }
 491 func fmul64c(f, g uint64, ret *uint64)            { *ret = fmul64(f, g) }
 492 func fdiv64c(f, g uint64, ret *uint64)            { *ret = fdiv64(f, g) }
 493 func fneg64c(f uint64, ret *uint64)               { *ret = fneg64(f) }
 494 func f32to64c(f uint32, ret *uint64)              { *ret = f32to64(f) }
 495 func f64to32c(f uint64, ret *uint32)              { *ret = f64to32(f) }
 496 func fcmp64c(f, g uint64, ret *int, retnan *bool) { *ret, *retnan = fcmp64(f, g) }
 497 func fintto64c(val int64, ret *uint64)            { *ret = fintto64(val) }
 498 func f64tointc(f uint64, ret *int64, retok *bool) { *ret, *retok = f64toint(f) }