我的程序是关于一个给定浮点的方法，在这种方法中，我想乘以或添加这些浮点数。但不是像 * 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 二进制 32 格式的更多信息，包括有效数的隐含前导 1，请参阅 https://en.wikipedia.org/wiki/Single-precision_floating-point_format。
并查看@njuffa对实际测试+工作实现的答案，由于某种原因，该实现将64位操作作为两个32位的一半，而不是让C有效地做到这一点。
此外，return *(float *)&result;违反了严格的别名。 它只在 MSVC 上是安全的。 在 C99/C11 中使用联合或 memcpy 进行类型双关。
模拟两个 IEEE-754 (2008)binary32操作数的乘法比问题所暗示的要复杂一些。一般来说，我们必须区分以下操作数类：零、次正态 (0 <|x| <2-126)、法线 (2 126 ≤ |x| <2128)、无穷大、NaN。法态在 [1， 254] 中使用有偏差的指数，而任何特殊操作数类在 {0， 255} 中使用有偏差的指数。下面假设我们要实现浮点乘法，屏蔽所有浮点异常，并使用舍入到最接近到偶数舍入模式。
首先，我们检查是否有任何参数属于特殊操作数类。如果是这样，我们将按顺序检查特殊情况。如果其中一个参数是 NaN，我们将该 NaN 转换为 QNaN 并返回它。如果其中一个操作数为零，我们返回一个适当签名的零，除非另一个参数是无穷大，在这种情况下，我们返回一个特殊的 QNaNINDEFINITE因为这是一个无效的操作。之后，我们检查无穷大的任何参数，返回一个适当签名的无穷大。这留下了次正态，我们将其归一化。如果有两个次正规参数，我们只需要规范化其中一个，因为结果将下溢到零。
法线的乘法按照提问者在问题中的设想进行。结果的符号是参数符号的异或，结果的指数是参数指数的总和(根据指数偏差进行调整)，结果的显着性由参数的显著性乘积生成。我们需要完整的产品进行四舍五入。我们可以为此使用 64 位类型，也可以用一对 32 位数字表示它。在下面的代码中，我选择了后一种表示形式。四舍五入到最接近或偶数很简单：如果我们有一个平局(结果正好在最接近的两个binary32数之间的中间)，我们需要四舍五入，如果尾数的最低有效位是 1。否则，如果最高有效丢弃位(舍入位)为 1，我们需要四舍五入。
根据舍入前的结果指数，结果需要考虑三种情况：指数在正常范围内、结果溢出(幅度太大)或下溢(幅度太小)。在第一种情况下，如果在舍入期间发生溢出，则结果为正常或无穷大。在第二种情况下，结果是无穷大。在最后一种情况下，结果为零(严重下溢)、次正态或最小正态(如果发生舍入)。
下面的代码，带有一个简单的框架，通过随机测试用例和数千个有趣的模式进行轻度测试，展示了一个在几个小时内编写的示例性 ISO-C 实现，以实现合理的清晰度和合理的性能。我让测试框架在 x64 平台上运行了一个小时左右，没有报告任何错误。如果计划在生产环境中使用代码，则需要构造更严格的测试框架，并且可能需要额外的性能调优。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <limits.h>
#define FLOAT_MANT_BITS (23)
#define FLOAT_EXPO_BITS (8)
#define FLOAT_EXPO_BIAS (127)
#define FLOAT_MANT_MASK (~((~0u) << (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 (~((~0u) << 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 (1u << 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)
{
unsigned long long 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 vectorsn");
#elif TEST_MODE == PATTERN_BASED
printf ("using pattern-based test vectorsn");
printf ("#patterns = %un", 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
i = KISS % patterns;
j = KISS % patterns;
a = uint_as_float ((v[i] & 0x7fffff) | (KISS & ~0x7fffff));
b = uint_as_float ((v[j] & 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)。克隆它并进行分析。