Source file src/pkg/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) }