Home | History | Annotate | Download | only in runtime
      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 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 
    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