我的程序是关于一个给出浮点数的方法，在这个方法中我想要乘以或添加那些浮点数。但不要像a * b那样相乘，我想将那些浮点数分解为它们的结构，如符号的位，指数的8位和其余的位作为尾数。
我想实现/模拟软件浮点加法和乘法（以了解更多关于FP硬件必须做什么）。
在计划的负责人有故障：
#define SIGN(x) (x>>31);
#define MANT(x) (x&0x7FFFFF);
#define EXPO(x) ((x>>23)&0xFF);
#define SPLIT(x, s, m, e) do { \
s = SIGN(x); \
m = MANT(x); \
e = EXPO(x); \
if ( e != 0x00 && e != 0xFF ) { \
m |= 0x800000; \
} \
} while ( 0 )
#define BUILD(x, s, m, e) do { \
x = (s << 31) | (e<<23) | (m&0x7FFFFF); \
} while ( 0 )
主要看起来如下：
float f = 2.3;
float g = 1.8;
float h = foo(&f, &g);
计算方法如下：
float foo(float *a, float *b) {
uint32_t ia = *(unsigned int *)a;
uint32_t ib = *(unsigned int *)b;
uint32_t result = 0;
uint32_t signa, signb, signr;
uint32_t manta, mantb, mantr;
uint32_t expoa, expob, expor;
SPLIT(ia, signa, manta, expoa);
SPLIT(ib, signb, mantb, expob);
我已经通过添加指数并乘以它们的尾数来尝试乘法如下：
expor = (expoa -127) + (expob -127) + 127;
mantr = (manta) * (mantb);
signr = signa ^ signb;
新浮点数的返回和重建：
BUILD(result, signr, mantr, expor);
return *(float *)&result;
问题是现在，结果是错误的。 mantr甚至需要一个非常低的负数（如果foo得到1.5和2.4 mantr需要-838860800，结果是2.0000000）。
你不能只截取尾数乘的结果，你需要取前24位（在使用低半部分进行舍入后）并重新规范化（调整指数）。
浮点运算保留顶部有效位。整数乘积中最重要的部分是高位;低位是小数点后的其他位置。 （术语：它是“二进制点”，而不是“小数点”，因为二进制浮点数使用基数2（二进制），而不是10（十进制）。）
对于归一化输入，输入有效数中的隐式前导1意味着用于实现24 x 24 => 48位尾数乘法的32x32 => 64位uint64_t乘积将在2个可能位置之一中具有其高位，因此你不需要进行位扫描就可以找到它。比较或单比特测试会做。
对于次正常输入，这不能保证，因此您需要检查MSB的位置，例如，与GNU C __builtin_clzll。 （有很多特殊情况需要处理一个或两个输入是否正常，和/或输出是否正常。）
有关IEEE-754 binary32格式的更多信息，请参阅https://en.wikipedia.org/wiki/Single-precision_floating-point_format，包括有效数字的隐含前导1。
请参阅@ njuffa的答案，了解实际经过测试的+工作实现，由于某种原因将64位操作作为两个32位的一半，而不是让C有效地执行此操作。
此外，return *(float *)&result;违反严格的别名。它在MSVC上是唯一安全的。在C99 / C11中使用union或memcpy进行类型惩罚。
模拟两个IEEE-754（2008）binary32操作数的乘法比问题暗示的要复杂一些。通常，我们必须区分以下操作数类：零，次正规（0 <| x | <2-126），法线（2126≤| x | <2128），无穷大，NaN。法线在[1,254]中使用有偏差的指数，而任何特殊操作数类在{0,255}中使用有偏差的指数。以下假设我们希望实现浮点乘法，并屏蔽所有浮点异常，并使用舍入到最接近舍入舍入模式。
首先，我们检查任何参数是否属于特殊操作数类。如果是这样，我们按顺序检查特殊情况。如果其中一个参数是NaN，我们返回将Nan转换为QNaN并返回它。如果其中一个操作数为零，则返回一个适当的有符号零，除非另一个参数是无穷大，在这种情况下，我们返回一个特殊的QNaN INDEFINITE，因为这是一个无效的操作。之后我们检查无限的任何参数，返回一个适当签名的无穷大。这留下了我们正常化的次正规值。如果有两个非正规参数，我们只需要对其中一个进行归一化，因为结果将是非正规或零。
法线的乘法按照提问者在问题中设想的那样进行。结果的符号是参数符号的异或，结果的指数是参数的指数的总和（根据指数偏差调整），结果的有效数是从产品的乘积生成的。重要的论点。我们需要完整的产品进行四舍五入。我们可以使用64位类型，或者用一对32位数字表示它。在下面的代码中，我选择了后一种表示法。舍入到最近或偶数是很简单的：如果我们有一个平局（结果恰好在最接近的两个binary32数之间的中间），我们需要向上舍入，如果尾数的最低有效位是1.否则如果最重要的丢弃位（圆位）为1，我们需要向上舍入。
根据舍入前的结果指数，需要考虑结果的三种情况：指数处于正常范围，结果溢出（幅度太大），或者下溢（幅度太小）。在第一种情况下，如果在舍入期间发生溢出，则结果为正常或无穷大。在第二种情况下，结果是无穷大。在最后的情况下，结果是零（严重的下溢），低于正常，或最小的正常（如果发生轮时）。
下面的代码，通过随机测试案例和几千个有趣模式的轻量级测试的简单框架，显示了在几个小时内编写的示例性ISO-C实现，以获得合理的清晰度和合理的性能。我让测试框架在x64平台上运行了一个小时左右，并且没有报告错误。如果您计划在生产中使用代码，则需要构建更严格的测试框架，并且可能需要进行额外的性能调整。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#define FLOAT_MANT_BITS (23)
#define FLOAT_EXPO_BITS (8)
#define FLOAT_EXPO_BIAS (127)
#define FLOAT_MANT_MASK (~((~0) << (FLOAT_MANT_BITS+1))) /* incl. integer bit */
#define EXPO_ADJUST (1) /* adjustment for performance reasons */
#define MIN_NORM_EXPO (1) /* minimum biased exponent of normals */
#define MAX_NORM_EXPO (254) /* maximum biased exponent of normals */
#define INF_EXPO (255) /* biased exponent of infinities */
#define EXPO_MASK (~((~0) << FLOAT_EXPO_BITS))
#define FLOAT_SIGN_MASK (0x80000000u)
#define FLOAT_IMPLICIT_BIT (1 << FLOAT_MANT_BITS)
#define RND_BIT_SHIFT (31)
#define RND_BIT_MASK (1 << RND_BIT_SHIFT)
#define FLOAT_INFINITY (0x7f800000)
#define FLOAT_INDEFINITE (0xffc00000u)
#define MANT_LSB (0x00000001)
#define FLOAT_QNAN_BIT (0x00400000)
#define MAX_SHIFT (FLOAT_MANT_BITS + 2)
uint32_t fp32_mul_core (uint32_t a, uint32_t b)
{
uint64_t prod;
uint32_t expoa, expob, manta, mantb, shift;
uint32_t r, signr, expor, mantr_hi, mantr_lo;
/* split arguments into sign, exponent, significand */
expoa = ((a >> FLOAT_MANT_BITS) & EXPO_MASK) - EXPO_ADJUST;
expob = ((b >> FLOAT_MANT_BITS) & EXPO_MASK) - EXPO_ADJUST;
manta = (a | FLOAT_IMPLICIT_BIT) & FLOAT_MANT_MASK;
mantb = (b | FLOAT_IMPLICIT_BIT) & FLOAT_MANT_MASK;
/* result sign bit: XOR sign argument signs */
signr = (a ^ b) & FLOAT_SIGN_MASK;
if ((expoa >= (MAX_NORM_EXPO - EXPO_ADJUST)) || /* at least one argument is special */
(expob >= (MAX_NORM_EXPO - EXPO_ADJUST))) {
if ((a & ~FLOAT_SIGN_MASK) > FLOAT_INFINITY) { /* a is NaN */
/* return quietened NaN */
return a | FLOAT_QNAN_BIT;
}
if ((b & ~FLOAT_SIGN_MASK) > FLOAT_INFINITY) { /* b is NaN */
/* return quietened NaN */
return b | FLOAT_QNAN_BIT;
}
if ((a & ~FLOAT_SIGN_MASK) == 0) { /* a is zero */
/* return NaN if b is infinity, else zero */
return (expob != (INF_EXPO - EXPO_ADJUST)) ? signr : FLOAT_INDEFINITE;
}
if ((b & ~FLOAT_SIGN_MASK) == 0) { /* b is zero */
/* return NaN if a is infinity, else zero */
return (expoa != (INF_EXPO - EXPO_ADJUST)) ? signr : FLOAT_INDEFINITE;
}
if (((a & ~FLOAT_SIGN_MASK) == FLOAT_INFINITY) || /* a or b infinity */
((b & ~FLOAT_SIGN_MASK) == FLOAT_INFINITY)) {
return signr | FLOAT_INFINITY;
}
if ((int32_t)expoa < (MIN_NORM_EXPO - EXPO_ADJUST)) { /* a is subnormal */
/* normalize significand of a */
manta = a & FLOAT_MANT_MASK;
expoa++;
do {
manta = 2 * manta;
expoa--;
} while (manta < FLOAT_IMPLICIT_BIT);
} else if ((int32_t)expob < (MIN_NORM_EXPO - EXPO_ADJUST)) { /* b is subnormal */
/* normalize significand of b */
mantb = b & FLOAT_MANT_MASK;
expob++;
do {
mantb = 2 * mantb;
expob--;
} while (mantb < FLOAT_IMPLICIT_BIT);
}
}
/* result exponent: add argument exponents and adjust for biasing */
expor = expoa + expob - FLOAT_EXPO_BIAS + 2 * EXPO_ADJUST;
mantb = mantb << FLOAT_EXPO_BITS; /* preshift to align result signficand */
/* result significand: multiply argument signficands */
prod = (uint64_t)manta * mantb;
mantr_hi = (uint32_t)(prod >> 32);
mantr_lo = (uint32_t)(prod >> 0);
/* normalize significand */
if (mantr_hi < FLOAT_IMPLICIT_BIT) {
mantr_hi = (mantr_hi << 1) | (mantr_lo >> (32 - 1));
mantr_lo = (mantr_lo << 1);
expor--;
}
if (expor <= (MAX_NORM_EXPO - EXPO_ADJUST)) { /* normal, may overflow to infinity during rounding */
/* combine biased exponent, sign and signficand */
r = (expor << FLOAT_MANT_BITS) + signr + mantr_hi;
/* round result to nearest or even; overflow to infinity possible */
r = r + ((mantr_lo == RND_BIT_MASK) ? (mantr_hi & MANT_LSB) : (mantr_lo >> RND_BIT_SHIFT));
} else if ((int32_t)expor > (MAX_NORM_EXPO - EXPO_ADJUST)) { /* overflow */
/* return infinity */
r = signr | FLOAT_INFINITY;
} else { /* underflow */
/* return zero, normal, or smallest subnormal */
shift = 0 - expor;
if (shift > MAX_SHIFT) shift = MAX_SHIFT;
/* denormalize significand */
mantr_lo = mantr_hi << (32 - shift) | (mantr_lo ? 1 : 0);
mantr_hi = mantr_hi >> shift;
/* combine sign and signficand; biased exponent known to be zero */
r = mantr_hi + signr;
/* round result to nearest or even */
r = r + ((mantr_lo == RND_BIT_MASK) ? (mantr_hi & MANT_LSB) : (mantr_lo >> RND_BIT_SHIFT));
}
return r;
}
uint32_t float_as_uint (float a)
{
uint32_t r;
memcpy (&r, &a, sizeof r);
return r;
}
float uint_as_float (uint32_t a)
{
float r;
memcpy (&r, &a, sizeof r);
return r;
}
float fp32_mul (float a, float b)
{
return uint_as_float (fp32_mul_core (float_as_uint (a), float_as_uint (b)));
}
/* Fixes via: Greg Rose, KISS: A Bit Too Simple. http://eprint.iacr.org/2011/007 */
static unsigned int z=362436069,w=521288629,jsr=362436069,jcong=123456789;
#define znew (z=36969*(z&0xffff)+(z>>16))
#define wnew (w=18000*(w&0xffff)+(w>>16))
#define MWC ((znew<<16)+wnew)
#define SHR3 (jsr^=(jsr<<13),jsr^=(jsr>>17),jsr^=(jsr<<5)) /* 2^32-1 */
#define CONG (jcong=69069*jcong+13579) /* 2^32 */
#define KISS ((MWC^CONG)+SHR3)
#define ISNAN(x) ((float_as_uint (x) << 1) > 0xff000000)
#define QNAN(x) (x | FLOAT_QNAN_BIT)
#define PURELY_RANDOM (0)
#define PATTERN_BASED (1)
#define TEST_MODE (PURELY_RANDOM)
uint32_t v[8192];
int main (void)
{
uint64_t count = 0;
float a, b, res, ref;
uint32_t i, j, patterns, idx = 0, nbrBits = sizeof (uint32_t) * CHAR_BIT;
/* pattern class 1: 2**i */
for (i = 0; i < nbrBits; i++) {
v [idx] = ((uint32_t)1 << i);
idx++;
}
/* pattern class 2: 2**i-1 */
for (i = 0; i < nbrBits; i++) {
v [idx] = (((uint32_t)1 << i) - 1);
idx++;
}
/* pattern class 3: 2**i+1 */
for (i = 0; i < nbrBits; i++) {
v [idx] = (((uint32_t)1 << i) + 1);
idx++;
}
/* pattern class 4: 2**i + 2**j */
for (i = 0; i < nbrBits; i++) {
for (j = 0; j < nbrBits; j++) {
v [idx] = (((uint32_t)1 << i) + ((uint32_t)1 << j));
idx++;
}
}
/* pattern class 5: 2**i - 2**j */
for (i = 0; i < nbrBits; i++) {
for (j = 0; j < nbrBits; j++) {
v [idx] = (((uint32_t)1 << i) - ((uint32_t)1 << j));
idx++;
}
}
/* pattern class 6: MAX_UINT/(2**i+1) rep. blocks of i zeros an i ones */
for (i = 0; i < nbrBits; i++) {
v [idx] = ((~(uint32_t)0) / (((uint32_t)1 << i) + 1));
idx++;
}
patterns = idx;
/* pattern class 6: one's complement of pattern classes 1 through 5 */
for (i = 0; i < patterns; i++) {
v [idx] = ~v [i];
idx++;
}
/* pattern class 7: two's complement of pattern classes 1 through 5 */
for (i = 0; i < patterns; i++) {
v [idx] = ~v [i] + 1;
idx++;
}
patterns = idx;
#if TEST_MODE == PURELY_RANDOM
printf ("using purely random test vectors\n");
#elif TEST_MODE == PATTERN_BASED
printf ("using pattern-based test vectors\n");
printf ("#patterns = %u\n", patterns);
#endif // TEST_MODE
do {
#if TEST_MODE == PURELY_RANDOM
a = uint_as_float (KISS);
b = uint_as_float (KISS);
#elif TEST_MODE == PATTERN_BASED
a = uint_as_float ((v[KISS%patterns] & 0x7fffff) | (KISS & ~0x7fffff));
b = uint_as_float ((v[KISS%patterns] & 0x7fffff) | (KISS & ~0x7fffff));
#endif // TEST_MODE
res = fp32_mul (a, b);
ref = a * b;
/* check for bit pattern mismatch between result and reference */
if (float_as_uint (res) != float_as_uint (ref)) {
/* if both a and b are NaNs, either could be returned quietened */
if (! (ISNAN (a) && ISNAN (b) &&
((QNAN (float_as_uint (a)) == float_as_uint (res)) ||
(QNAN (float_as_uint (b)) == float_as_uint (res))))) {
printf ("err: a=% 15.8e (%08x) b=% 15.8e (%08x) res=% 15.8e (%08x) ref=%15.8e (%08x)\n",
a, float_as_uint(a), b, float_as_uint (b), res, float_as_uint (res), ref, float_as_uint (ref));
return EXIT_FAILURE;
}
}
count++;
if (!(count & 0xffffff)) printf ("\r%llu", count);
} while (1);
return EXIT_SUCCESS;
}
它要复杂得多。看一下softmath库的源代码（例如https://github.com/riscv/riscv-pk/blob/master/softfloat/f64_mul.c）。克隆并分析。