将 IEEE 32 位浮点 [1:2) 的范围映射到任意 [a:b)
map range of IEEE 32bit float [1:2) to some arbitrary [a:b)
背景故事：具有任意端点的统一 PRNG
我有一个快速统一的伪随机数生成器，可以在 [1:2) 范围内创建统一的 float32 数字，即 u : 1 <= u <= 2-eps。不幸的是，将端点 [1:2) 映射到任意范围 [a:b) 的端点在浮点数学中并不简单。我想通过简单的仿射计算来精确匹配端点。
正式说明
我想为 1<=x<2 和 精确 映射的任意 a,b 创建一个 IEEE-754 32 位浮点仿射函数 f(x,a,b)
1 -> a 和 nextlower(2) -> nextlower(b)
其中 nextlower(q) 是下一个较低的 FP 可表示数（例如在 C++ 中 std::nextafter(float(q),float(q-1))）
我试过的
简单映射f(x,a,b) = (x-1)*(b-a) + a总是满足f(1)条件但有时由于浮点舍入而无法满足f(2)条件。
我尝试用自由设计参数替换 1，以本着 Kahan summation 的精神取消 FP 错误。
即与
f(x,c0,c1,c2) = (x-c0)*c1 + c2
一个数学解是c0=1,c1=(b-a),c2=a（上面的简单映射），
但是额外的参数让我可以使用常量 c0,c1,c2 来匹配端点。我不确定我是否充分理解 Kahan 求和背后的原理以应用它们来确定参数或什至确信存在解决方案。感觉就像我在黑暗中颠簸，其他人可能已经找到了光。
旁白：我可以假设以下内容
- a < b
- a 和 b 都远离零，即可以忽略次正规
- a 和 b 足够远（以可表示的 FP 值测量）以减轻非均匀量化并避免退化情况
更新
我正在使用 Chux 答案的修改形式来避免分裂。
虽然我不能 100% 确定我的重构保留了所有的魔力，但它在我所有的测试用例中仍然有效。
float lerp12(float x,float a,float b)
{
const float scale = 1.0000001f;
// scale = 1/(nextlower(2) - 1);
const float ascale = a*scale;
const float bscale = nextlower(b)*scale;
return (nextlower(2) - x)*ascale + (x - 1.0f)*bscale;
}
请注意，只有最后一行 (5 FLOPS) 取决于 x，因此如果 (a,b) 保持不变，则可以重复使用其他行。
基于融合乘加的简单 lerping 可以可靠地命中插值因子 0 和 1 的端点。对于 [1, 2) 中的 x，插值因子 x - 1 未达到统一，可以通过将 x-1 乘以 (2.0f / nextlower(2.0f)) 进行轻微拉伸来修复。显然端点也需要调整为端点nextlower(b)。对于下面的 C 代码，我使用了问题中提供的 nextlower() 的定义，这可能不是提问者想要的，因为对于浮点数 q 足够大，q == (q - 1)。
Asker 在评论中指出，据了解，这种映射不会导致伪随机数在区间 [a, b) 中完全均匀分布，只是近似如此，而且病态当 a 和 b 非常接近时，可能会发生映射。我还没有从数学上证明下面 map() 的实现可以保证期望的行为，但它似乎对大量随机测试用例是这样的。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
float nextlowerf (float q)
{
return nextafterf (q, q - 1);
}
float map (float a, float b, float x)
{
float t = (x - 1.0f) * (2.0f / nextlowerf (2.0f));
return fmaf (t, nextlowerf (b), fmaf (-t, a, a));
}
float uint32_as_float (uint32_t a)
{
float r;
memcpy (&r, &a, sizeof(r));
return r;
}
// George Marsaglia's KISS PRNG, period 2**123. Newsgroup sci.math, 21 Jan 1999
// Bug fix: Greg Rose, "KISS: A Bit Too Simple" http://eprint.iacr.org/2011/007
static uint32_t kiss_z=362436069, kiss_w=521288629;
static uint32_t kiss_jsr=123456789, kiss_jcong=380116160;
#define znew (kiss_z=36969*(kiss_z&65535)+(kiss_z>>16))
#define wnew (kiss_w=18000*(kiss_w&65535)+(kiss_w>>16))
#define MWC ((znew<<16)+wnew )
#define SHR3 (kiss_jsr^=(kiss_jsr<<13),kiss_jsr^=(kiss_jsr>>17), \
kiss_jsr^=(kiss_jsr<<5))
#define CONG (kiss_jcong=69069*kiss_jcong+1234567)
#define KISS ((MWC^CONG)+SHR3)
int main (void)
{
float a, b, x, r;
float FP32_MIN_NORM = 0x1.000000p-126f;
float FP32_MAX_NORM = 0x1.fffffep+127f;
do {
do {
a = uint32_as_float (KISS);
} while ((fabsf (a) < FP32_MIN_NORM) || (fabsf (a) > FP32_MAX_NORM) || isnan (a));
do {
b = uint32_as_float (KISS);
} while ((fabsf (b) < FP32_MIN_NORM) || (fabsf (b) > FP32_MAX_NORM) || isnan (b) || (b < a));
x = 1.0f;
r = map (a, b, x);
if (r != a) {
printf ("lower bound failed: a=%12.6a b=%12.6a map=%12.6a\n", a, b, r);
return EXIT_FAILURE;
}
x = nextlowerf (2.0f);
r = map (a, b, x);
if (r != nextlowerf (b)) {
printf ("upper bound failed: a=%12.6a b=%12.6a map=%12.6a\n", a, b, r);
return EXIT_FAILURE;
}
} while (1);
return EXIT_SUCCESS;
}
OP的目标
I want to make an IEEE-754 32 bit floating point affine function f(x,a,b) for 1<=x<2 and arbitrary a,b that exactly maps 1 -> a and nextlower(2) -> nextlower(b)
这与“将 IEEE 32 位浮点数 [1:2) 映射到某个任意 [a:b)”略有不同。
一般情况
将 x0 映射到 y0，将 x1 映射到 y1 以及介于两者之间的各种 x 到 y :
m = (y1 - y0)/(x1 - x0);
y = m*(x - x0) + y0;
OP案例
// x0 = 1.0f;
// x1 = nextafterf(2.0f, 1.0f);
// y0 = a;
// y1 = nextafterf(b, a);
#include <math.h> // for nextafterf()
float x = random_number_1_to_almost_2();
float m = (nextafterf(b, a) - a)/(nextafterf(2.0f, 1.0f) - 1.0f);
float y = m*(x - 1.0f) + a;
nextafterf(2.0f, 1.0f) - 1.0f、x - 1.0f和nextafterf(b, a)是准确的，没有计算错误。
nextafterf(2.0f, 1.0f) - 1.0f 是一个略小于 1.0f 的值。
推荐
在端点具有更好的对称性和数值稳定性的其他重组是可能的。
float x = random_number_1_to_almost_2();
float afactor = nextafterf(2.0f, 1.0f) - x; // exact
float bfactor = x - 1.0f; // exact
float xwidth = nextafterf(2.0f, 1.0f) - 1.0f; // exact
// Do not re-order next line of code, perform 2 divisions
float y = (afactor/xwidth)*a + (bfactor/xwidth)*nextafterf(b, a);
注意 afactor/xwidth 和 bfactor/xwidth 在端点处都恰好为 0.0 或 1.0，因此满足“映射 1 -> a 和 nextlower(2) -> nextlower(b)”。不需要扩展精度。
OP 的 (x-c0)*c1 + c2 在将 (x-c0)*c1 除以 (2.0 - 1.0) 或 1.0（隐含）时出现问题，而它应该除以 nextafterf(2.0f, 1.0f) - 1.0f。
背景故事：具有任意端点的统一 PRNG
我有一个快速统一的伪随机数生成器，可以在 [1:2) 范围内创建统一的 float32 数字，即 u : 1 <= u <= 2-eps。不幸的是，将端点 [1:2) 映射到任意范围 [a:b) 的端点在浮点数学中并不简单。我想通过简单的仿射计算来精确匹配端点。
正式说明
我想为 1<=x<2 和 精确 映射的任意 a,b 创建一个 IEEE-754 32 位浮点仿射函数 f(x,a,b)
1 -> a 和 nextlower(2) -> nextlower(b)
其中 nextlower(q) 是下一个较低的 FP 可表示数（例如在 C++ 中 std::nextafter(float(q),float(q-1))）
我试过的
简单映射f(x,a,b) = (x-1)*(b-a) + a总是满足f(1)条件但有时由于浮点舍入而无法满足f(2)条件。
我尝试用自由设计参数替换 1，以本着 Kahan summation 的精神取消 FP 错误。
即与
f(x,c0,c1,c2) = (x-c0)*c1 + c2
一个数学解是c0=1,c1=(b-a),c2=a（上面的简单映射），
但是额外的参数让我可以使用常量 c0,c1,c2 来匹配端点。我不确定我是否充分理解 Kahan 求和背后的原理以应用它们来确定参数或什至确信存在解决方案。感觉就像我在黑暗中颠簸，其他人可能已经找到了光。
旁白：我可以假设以下内容
- a < b
- a 和 b 都远离零，即可以忽略次正规
- a 和 b 足够远（以可表示的 FP 值测量）以减轻非均匀量化并避免退化情况
更新
我正在使用 Chux 答案的修改形式来避免分裂。 虽然我不能 100% 确定我的重构保留了所有的魔力，但它在我所有的测试用例中仍然有效。
float lerp12(float x,float a,float b)
{
const float scale = 1.0000001f;
// scale = 1/(nextlower(2) - 1);
const float ascale = a*scale;
const float bscale = nextlower(b)*scale;
return (nextlower(2) - x)*ascale + (x - 1.0f)*bscale;
}
请注意，只有最后一行 (5 FLOPS) 取决于 x，因此如果 (a,b) 保持不变，则可以重复使用其他行。
基于融合乘加的简单 lerping 可以可靠地命中插值因子 0 和 1 的端点。对于 [1, 2) 中的 x，插值因子 x - 1 未达到统一，可以通过将 x-1 乘以 (2.0f / nextlower(2.0f)) 进行轻微拉伸来修复。显然端点也需要调整为端点nextlower(b)。对于下面的 C 代码，我使用了问题中提供的 nextlower() 的定义，这可能不是提问者想要的，因为对于浮点数 q 足够大，q == (q - 1)。
Asker 在评论中指出，据了解，这种映射不会导致伪随机数在区间 [a, b) 中完全均匀分布，只是近似如此，而且病态当 a 和 b 非常接近时，可能会发生映射。我还没有从数学上证明下面 map() 的实现可以保证期望的行为，但它似乎对大量随机测试用例是这样的。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
float nextlowerf (float q)
{
return nextafterf (q, q - 1);
}
float map (float a, float b, float x)
{
float t = (x - 1.0f) * (2.0f / nextlowerf (2.0f));
return fmaf (t, nextlowerf (b), fmaf (-t, a, a));
}
float uint32_as_float (uint32_t a)
{
float r;
memcpy (&r, &a, sizeof(r));
return r;
}
// George Marsaglia's KISS PRNG, period 2**123. Newsgroup sci.math, 21 Jan 1999
// Bug fix: Greg Rose, "KISS: A Bit Too Simple" http://eprint.iacr.org/2011/007
static uint32_t kiss_z=362436069, kiss_w=521288629;
static uint32_t kiss_jsr=123456789, kiss_jcong=380116160;
#define znew (kiss_z=36969*(kiss_z&65535)+(kiss_z>>16))
#define wnew (kiss_w=18000*(kiss_w&65535)+(kiss_w>>16))
#define MWC ((znew<<16)+wnew )
#define SHR3 (kiss_jsr^=(kiss_jsr<<13),kiss_jsr^=(kiss_jsr>>17), \
kiss_jsr^=(kiss_jsr<<5))
#define CONG (kiss_jcong=69069*kiss_jcong+1234567)
#define KISS ((MWC^CONG)+SHR3)
int main (void)
{
float a, b, x, r;
float FP32_MIN_NORM = 0x1.000000p-126f;
float FP32_MAX_NORM = 0x1.fffffep+127f;
do {
do {
a = uint32_as_float (KISS);
} while ((fabsf (a) < FP32_MIN_NORM) || (fabsf (a) > FP32_MAX_NORM) || isnan (a));
do {
b = uint32_as_float (KISS);
} while ((fabsf (b) < FP32_MIN_NORM) || (fabsf (b) > FP32_MAX_NORM) || isnan (b) || (b < a));
x = 1.0f;
r = map (a, b, x);
if (r != a) {
printf ("lower bound failed: a=%12.6a b=%12.6a map=%12.6a\n", a, b, r);
return EXIT_FAILURE;
}
x = nextlowerf (2.0f);
r = map (a, b, x);
if (r != nextlowerf (b)) {
printf ("upper bound failed: a=%12.6a b=%12.6a map=%12.6a\n", a, b, r);
return EXIT_FAILURE;
}
} while (1);
return EXIT_SUCCESS;
}
OP的目标
I want to make an IEEE-754 32 bit floating point affine function f(x,a,b) for 1<=x<2 and arbitrary a,b that exactly maps 1 -> a and nextlower(2) -> nextlower(b)
这与“将 IEEE 32 位浮点数 [1:2) 映射到某个任意 [a:b)”略有不同。
一般情况
将 x0 映射到 y0，将 x1 映射到 y1 以及介于两者之间的各种 x 到 y :
m = (y1 - y0)/(x1 - x0);
y = m*(x - x0) + y0;
OP案例
// x0 = 1.0f;
// x1 = nextafterf(2.0f, 1.0f);
// y0 = a;
// y1 = nextafterf(b, a);
#include <math.h> // for nextafterf()
float x = random_number_1_to_almost_2();
float m = (nextafterf(b, a) - a)/(nextafterf(2.0f, 1.0f) - 1.0f);
float y = m*(x - 1.0f) + a;
nextafterf(2.0f, 1.0f) - 1.0f、x - 1.0f和nextafterf(b, a)是准确的，没有计算错误。
nextafterf(2.0f, 1.0f) - 1.0f 是一个略小于 1.0f 的值。
推荐
在端点具有更好的对称性和数值稳定性的其他重组是可能的。
float x = random_number_1_to_almost_2();
float afactor = nextafterf(2.0f, 1.0f) - x; // exact
float bfactor = x - 1.0f; // exact
float xwidth = nextafterf(2.0f, 1.0f) - 1.0f; // exact
// Do not re-order next line of code, perform 2 divisions
float y = (afactor/xwidth)*a + (bfactor/xwidth)*nextafterf(b, a);
注意 afactor/xwidth 和 bfactor/xwidth 在端点处都恰好为 0.0 或 1.0，因此满足“映射 1 -> a 和 nextlower(2) -> nextlower(b)”。不需要扩展精度。
OP 的 (x-c0)*c1 + c2 在将 (x-c0)*c1 除以 (2.0 - 1.0) 或 1.0（隐含）时出现问题，而它应该除以 nextafterf(2.0f, 1.0f) - 1.0f。