如何正确实现浮点数乘法（软件 FP）
How to correctly implement multiply for floating point numbers (software FP)
我的程序是关于一个给定浮点数的方法，在这个方法中我想乘或加这些浮点数。但不像 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;
新浮动的return和重建：
BUILD(result, signr, mantr, expor);
return *(float *)&result;
现在的问题是，结果是错误的。 mantr 甚至取一个非常低的负数（如果 foo 得到 1.5 和 2.4 mantr 取 -838860800 并且结果为 2.0000000）。
要复杂得多。查看 softmath 库的源代码（例如 https://github.com/riscv/riscv-pk/blob/master/softfloat/f64_mul.c）。克隆它并分析。
尾数相乘的结果不能直接truncate，需要把top24位（用低半部分四舍五入后）重新归一化（调整指数）。
浮点运算保留最高有效位。整数乘积的最高有效部分是高位；低位在小数点后更远的位置。 （术语：它是 "binary point"，而不是 "decimal point"，因为二进制浮点数使用基数 2（二进制），而不是 10（十进制）。）
对于归一化输入，输入尾数中的隐式前导 1 表示用于实现 24 x 24 => 48 位尾数乘法的 32x32 => 64 位 uint64_t 乘积将在 2 个可能的位置之一有它的高位，所以你不需要位扫描来找到它。比较或单比特测试就可以了。
对于次正规输入，不能保证这一点，因此您需要检查 MSB 的位置，例如使用 GNU C __builtin_clzll。 （有许多特殊情况需要处理，一个或两个输入是次正规的，and/or输出是次正规的。）
有关 IEEE-754 binary32 格式的更多信息，包括尾数的隐含前导 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 操作数的乘法比问题所暗示的要复杂一些。一般来说，我们要区分以下操作数classes：zeros, subnormals (0 < |x| < 2-126), normals (2 126 ≤ |x| < 2128)，无穷大，NaN。法线在 [1, 254] 中使用有偏指数，而任何特殊操作数 classes 在 {0, 255} 中使用有偏指数。下面假设我们要实现浮点乘法，屏蔽所有浮点异常，并使用舍入到最近到偶数舍入模式。
首先，我们检查是否有任何参数属于特殊操作数class。如果是这样，我们按顺序检查特殊情况。如果参数之一是 NaN，我们将该 NaN 转换为 QNaN 并 return 它。如果其中一个操作数是零，我们 return 一个适当签名的零， 除非 另一个参数是无穷大，在这种情况下我们 return 一个特殊的 QNaN INDEFINITE 因为这是一个无效的操作。之后我们检查无穷大的任何参数，returning 一个适当签名的无穷大。这留下了次正规，我们将其归一化。如果有两个次正规参数，我们只需要对其中一个进行归一化，因为结果将下溢为零。
法线的乘法按照提问者在问题中的设想进行。结果的符号是参数符号的异或，结果的指数是参数的指数之和（针对指数偏差进行了调整），结果的有效数是从乘积生成的论据的重要性。我们需要四舍五入的完整产品。我们可以为此使用 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 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
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;
}
我的程序是关于一个给定浮点数的方法，在这个方法中我想乘或加这些浮点数。但不像 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;
新浮动的return和重建：
BUILD(result, signr, mantr, expor);
return *(float *)&result;
现在的问题是，结果是错误的。 mantr 甚至取一个非常低的负数（如果 foo 得到 1.5 和 2.4 mantr 取 -838860800 并且结果为 2.0000000）。
要复杂得多。查看 softmath 库的源代码（例如 https://github.com/riscv/riscv-pk/blob/master/softfloat/f64_mul.c）。克隆它并分析。
尾数相乘的结果不能直接truncate，需要把top24位（用低半部分四舍五入后）重新归一化（调整指数）。
浮点运算保留最高有效位。整数乘积的最高有效部分是高位；低位在小数点后更远的位置。 （术语：它是 "binary point"，而不是 "decimal point"，因为二进制浮点数使用基数 2（二进制），而不是 10（十进制）。）
对于归一化输入，输入尾数中的隐式前导 1 表示用于实现 24 x 24 => 48 位尾数乘法的 32x32 => 64 位 uint64_t 乘积将在 2 个可能的位置之一有它的高位，所以你不需要位扫描来找到它。比较或单比特测试就可以了。
对于次正规输入，不能保证这一点，因此您需要检查 MSB 的位置，例如使用 GNU C __builtin_clzll。 （有许多特殊情况需要处理，一个或两个输入是次正规的，and/or输出是次正规的。）
有关 IEEE-754 binary32 格式的更多信息，包括尾数的隐含前导 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 操作数的乘法比问题所暗示的要复杂一些。一般来说，我们要区分以下操作数classes：zeros, subnormals (0 < |x| < 2-126), normals (2 126 ≤ |x| < 2128)，无穷大，NaN。法线在 [1, 254] 中使用有偏指数，而任何特殊操作数 classes 在 {0, 255} 中使用有偏指数。下面假设我们要实现浮点乘法，屏蔽所有浮点异常，并使用舍入到最近到偶数舍入模式。
首先，我们检查是否有任何参数属于特殊操作数class。如果是这样，我们按顺序检查特殊情况。如果参数之一是 NaN，我们将该 NaN 转换为 QNaN 并 return 它。如果其中一个操作数是零，我们 return 一个适当签名的零， 除非 另一个参数是无穷大，在这种情况下我们 return 一个特殊的 QNaN INDEFINITE 因为这是一个无效的操作。之后我们检查无穷大的任何参数，returning 一个适当签名的无穷大。这留下了次正规，我们将其归一化。如果有两个次正规参数，我们只需要对其中一个进行归一化，因为结果将下溢为零。
法线的乘法按照提问者在问题中的设想进行。结果的符号是参数符号的异或，结果的指数是参数的指数之和（针对指数偏差进行了调整），结果的有效数是从乘积生成的论据的重要性。我们需要四舍五入的完整产品。我们可以为此使用 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 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
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;
}