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  }
   628
View as plain text