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