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 9package runtime 10 11const ( 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 29func 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 61func 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 93func 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 144func 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 195func 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 250func fsub64(f, g uint64) uint64 { 251 return fadd64(f, fneg64(g)) 252} 253 254func fneg64(f uint64) uint64 { 255 return f ^ (1 << (mantbits64 + expbits64)) 256} 257 258func 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 288func 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 317func 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 330func 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 343func 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 374func 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 409func 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. 420func 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 439func 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 462again1: 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 475again2: 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 487func fadd32(x, y uint32) uint32 { 488 return f64to32(fadd64(f32to64(x), f32to64(y))) 489} 490 491func fmul32(x, y uint32) uint32 { 492 return f64to32(fmul64(f32to64(x), f32to64(y))) 493} 494 495func fdiv32(x, y uint32) uint32 { 496 return f64to32(fdiv64(f32to64(x), f32to64(y))) 497} 498 499func feq32(x, y uint32) bool { 500 cmp, nan := fcmp64(f32to64(x), f32to64(y)) 501 return cmp == 0 && !nan 502} 503 504func fgt32(x, y uint32) bool { 505 cmp, nan := fcmp64(f32to64(x), f32to64(y)) 506 return cmp >= 1 && !nan 507} 508 509func fge32(x, y uint32) bool { 510 cmp, nan := fcmp64(f32to64(x), f32to64(y)) 511 return cmp >= 0 && !nan 512} 513 514func feq64(x, y uint64) bool { 515 cmp, nan := fcmp64(x, y) 516 return cmp == 0 && !nan 517} 518 519func fgt64(x, y uint64) bool { 520 cmp, nan := fcmp64(x, y) 521 return cmp >= 1 && !nan 522} 523 524func fge64(x, y uint64) bool { 525 cmp, nan := fcmp64(x, y) 526 return cmp >= 0 && !nan 527} 528 529func fint32to32(x int32) uint32 { 530 return f64to32(fintto64(int64(x))) 531} 532 533func fint32to64(x int32) uint64 { 534 return fintto64(int64(x)) 535} 536 537func fint64to32(x int64) uint32 { 538 return f64to32(fintto64(x)) 539} 540 541func fint64to64(x int64) uint64 { 542 return fintto64(x) 543} 544 545func f32toint32(x uint32) int32 { 546 val, _ := f64toint(f32to64(x)) 547 return int32(val) 548} 549 550func f32toint64(x uint32) int64 { 551 val, _ := f64toint(f32to64(x)) 552 return val 553} 554 555func f64toint32(x uint64) int32 { 556 val, _ := f64toint(x) 557 return int32(val) 558} 559 560func f64toint64(x uint64) int64 { 561 val, _ := f64toint(x) 562 return val 563} 564 565func f64touint64(x float64) uint64 { 566 if x < float64(1<<63) { 567 return uint64(int64(x)) 568 } 569 y := x - float64(1<<63) 570 z := uint64(int64(y)) 571 return z | (1 << 63) 572} 573 574func f32touint64(x float32) uint64 { 575 if x < float32(1<<63) { 576 return uint64(int64(x)) 577 } 578 y := x - float32(1<<63) 579 z := uint64(int64(y)) 580 return z | (1 << 63) 581} 582 583func fuint64to64(x uint64) float64 { 584 if int64(x) >= 0 { 585 return float64(int64(x)) 586 } 587 // See ../cmd/compile/internal/gc/ssa.go:uint64Tofloat 588 y := x & 1 589 z := x >> 1 590 z = z | y 591 r := float64(int64(z)) 592 return r + r 593} 594 595func fuint64to32(x uint64) float32 { 596 return float32(fuint64to64(x)) 597} 598