浮点范围缩减
Floating point range reduction
我正在使用 Mono 在 C# 中实现一些 32 位浮点三角函数,希望利用 Mono.Simd。我目前只缺少固体范围减少。
我现在很困惑,因为显然 Mono's SIMD extensions 不包括浮点数和整数之间的转换,这意味着我无法访问通常的方法 rounding/truncation 。但是我可以在整数和浮点数之间按位转换。
这样的事情能做吗?如果需要,我可以上下缩放域,但理想情况下,范围缩小应该导致域为 [0, 2 pi] 或 [-pi, pi]。我有一种预感,如果定义域是 2 的幂,就可以用指数做一些 IEEE 魔术,但我真的不确定该怎么做。
编辑:
好吧,我试过弄乱这个 C 代码,感觉就像我在某事的边缘(它不起作用但小数部分总是正确的,至少在十进制/ base10 中......)。核心原则似乎是获取你的域和输入指数之间的指数差异,并用移位的尾数和调整后的指数组成一个新的浮点数。但它不适用于负数,我不知道如何处理2 的非幂(或任何小数 - 事实上,除 2 以外的任何东西都不起作用!)。
// here's another more correct attempt:
float fmodulus(float val, int domain)
{
const int mantissaMask = 0x7FFFFF;
const int exponentMask = 0x7F800000;
int ival = *(int*)&val;
int mantissa = ival & mantissaMask;
int rawExponent = ival & exponentMask;
int exponent = (rawExponent >> 23) - (129 - domain);
// powers over one:
int p = exponent;
mantissa <<= p;
rawExponent = exponent >> p;
rawExponent += 127;
rawExponent <<= 23;
int newVal = rawExponent & exponentMask;
newVal |= mantissa & mantissaMask;
float ret = *(float*)&newVal;
return ret;
}
float range_reduce(float value, int range )
{
const int mantissaMask = 0x7FFFFF;
const int exponentMask = 0x7F800000;
int ival = *(int*)&value;
// grab exponent:
unsigned exponent = (ival & exponentMask) >> 23;
// grab mantissa:
unsigned mantissa = ival & mantissaMask;
// remove bias, and see how much the exponent is over range/domain
unsigned char erange = (unsigned char)(exponent - (125 + range));
// check if sign bit is set - that is, the exponent is under our range
if (erange & 0x80)
{
// don't do anything then.
erange = 0;
}
// shift mantissa (and chop off bits) by the reduced amount
int inewVal = (mantissa << (erange)) & mantissaMask;
// add exponent, and subtract the amount we reduced the argument with
inewVal |= ((exponent - erange) << 23) & exponentMask;
// reinterpret
float newValue = *(float*)&inewVal;
return newValue;
//return newValue - ((erange) & 0x1 ? 1.0f : 0.0f);
}
int main()
{
float val = 2.687f;
int ival = *(int*)&val;
float correct = fmod(val, 2);
float own = range_reduce(val, 2);
getc(stdin);
}
编辑 2:
好的,我真的想从 ieee 二进制系统的角度来理解这一点。如果我们这样写模运算:
output = input % 2
[exponent] + [mantissa_bit_n_times_exponent]
3.5 = [2] + [1 + 0.5] ->[1] + [0.5] = 1.5
4.5 = [4] + [0 + 0 + 0.5] ->[0.5] + [0] = 0.5
5.5 = [4] + [0 + 1 + 0.5] ->[1] + [0.5] = 1.5
2.5 = [2] + [0 + 0.5] ->[0.5] + [0] = 0.5
2.25 = [2] + [0 + 0 + 0.25] ->[0.25] = 0.25
2.375 = [2] + [0 + 0 + 0.25 + 0.125] ->[0.25] + [0.125] = 0.375
13.5 = [8] + [4 + 0 + 1 + 0.5] ->[1] + [0.5] = 1.5
56.5 = [32] + [16 + 8 + 0 + 0 + 0 + 0.5] ->[0.5] = 0.5
我们可以看到所有情况下的输出都是一个新数,没有原始指数,尾数移动了一个量(这是基于指数和尾数的第一个非零位在尾数的第一个指数位被忽略后 ) 进入指数。但我不确定这是否是正确的方法,它在纸面上效果很好。
编辑3:
我卡在单声道版本 2.0.50727.1433
检查您的单声道版本,因为 ConvertToInt
和 ConvertToIntTruncated
have been added 4 years ago 应该在 2.10 版本后就存在。
您可以将问题简化为取一个浮点数 mod 1. 为简化这一点,您可以使用位运算计算浮点数的底数,然后使用浮点数减法。以下是这些操作的(不安全)C# 代码:
// domain is assumed to be positive
// returns value in [0,domain)
public float fmodulus(float val, float domain)
{
if (val < 0)
{
float negative = fmodulus(-val, domain);
if (domain - negative == domain)
return 0;
else
return domain-negative;
}
if (val < domain)
return val; // this avoids losing accuracy
return fmodOne(val / domain) * domain;
}
// assumes val >= 1, so val is positive and the exponent is at least 0
unsafe public float fmodOne(float val)
{
int iVal = *(int*)&val;
int uncenteredExponent = iVal >> 23;
int exponent = uncenteredExponent - 127; // 127 corresponds to 2^0 times the mantissa
if (exponent >= 23)
return 0; // not enough precision to distinguish val from an integer
int unneededBits = 23 - exponent; // between 0 and 23
int iFloorVal = (iVal >> unneededBits) << unneededBits; // equivalent to using a mask to zero the bottom bits of the mantissa
float floorVal = *(float*)&iFloorVal; // convert the bit pattern back to a float
return val-floorVal;
}
例如,fmodulus(100.1f, 1) 为 0.09999847。 100.1f 的位模式是
0 10000101 10010000011001100110011
floorVal (100f) 的位模式是
0 10000101 10010000000000000000000
浮点数减法给出接近 0.1f 的值:
0 01111011 10011001100110000000000
实际上,我很惊讶最后8位被清零了。我认为只有 0.1f 的最后 6 位应该被替换为 0。也许可以比依赖浮点减法做得更好。
我正在使用 Mono 在 C# 中实现一些 32 位浮点三角函数,希望利用 Mono.Simd。我目前只缺少固体范围减少。 我现在很困惑,因为显然 Mono's SIMD extensions 不包括浮点数和整数之间的转换,这意味着我无法访问通常的方法 rounding/truncation 。但是我可以在整数和浮点数之间按位转换。
这样的事情能做吗?如果需要,我可以上下缩放域,但理想情况下,范围缩小应该导致域为 [0, 2 pi] 或 [-pi, pi]。我有一种预感,如果定义域是 2 的幂,就可以用指数做一些 IEEE 魔术,但我真的不确定该怎么做。
编辑: 好吧,我试过弄乱这个 C 代码,感觉就像我在某事的边缘(它不起作用但小数部分总是正确的,至少在十进制/ base10 中......)。核心原则似乎是获取你的域和输入指数之间的指数差异,并用移位的尾数和调整后的指数组成一个新的浮点数。但它不适用于负数,我不知道如何处理2 的非幂(或任何小数 - 事实上,除 2 以外的任何东西都不起作用!)。
// here's another more correct attempt:
float fmodulus(float val, int domain)
{
const int mantissaMask = 0x7FFFFF;
const int exponentMask = 0x7F800000;
int ival = *(int*)&val;
int mantissa = ival & mantissaMask;
int rawExponent = ival & exponentMask;
int exponent = (rawExponent >> 23) - (129 - domain);
// powers over one:
int p = exponent;
mantissa <<= p;
rawExponent = exponent >> p;
rawExponent += 127;
rawExponent <<= 23;
int newVal = rawExponent & exponentMask;
newVal |= mantissa & mantissaMask;
float ret = *(float*)&newVal;
return ret;
}
float range_reduce(float value, int range )
{
const int mantissaMask = 0x7FFFFF;
const int exponentMask = 0x7F800000;
int ival = *(int*)&value;
// grab exponent:
unsigned exponent = (ival & exponentMask) >> 23;
// grab mantissa:
unsigned mantissa = ival & mantissaMask;
// remove bias, and see how much the exponent is over range/domain
unsigned char erange = (unsigned char)(exponent - (125 + range));
// check if sign bit is set - that is, the exponent is under our range
if (erange & 0x80)
{
// don't do anything then.
erange = 0;
}
// shift mantissa (and chop off bits) by the reduced amount
int inewVal = (mantissa << (erange)) & mantissaMask;
// add exponent, and subtract the amount we reduced the argument with
inewVal |= ((exponent - erange) << 23) & exponentMask;
// reinterpret
float newValue = *(float*)&inewVal;
return newValue;
//return newValue - ((erange) & 0x1 ? 1.0f : 0.0f);
}
int main()
{
float val = 2.687f;
int ival = *(int*)&val;
float correct = fmod(val, 2);
float own = range_reduce(val, 2);
getc(stdin);
}
编辑 2:
好的,我真的想从 ieee 二进制系统的角度来理解这一点。如果我们这样写模运算:
output = input % 2
[exponent] + [mantissa_bit_n_times_exponent]
3.5 = [2] + [1 + 0.5] ->[1] + [0.5] = 1.5
4.5 = [4] + [0 + 0 + 0.5] ->[0.5] + [0] = 0.5
5.5 = [4] + [0 + 1 + 0.5] ->[1] + [0.5] = 1.5
2.5 = [2] + [0 + 0.5] ->[0.5] + [0] = 0.5
2.25 = [2] + [0 + 0 + 0.25] ->[0.25] = 0.25
2.375 = [2] + [0 + 0 + 0.25 + 0.125] ->[0.25] + [0.125] = 0.375
13.5 = [8] + [4 + 0 + 1 + 0.5] ->[1] + [0.5] = 1.5
56.5 = [32] + [16 + 8 + 0 + 0 + 0 + 0.5] ->[0.5] = 0.5
我们可以看到所有情况下的输出都是一个新数,没有原始指数,尾数移动了一个量(这是基于指数和尾数的第一个非零位在尾数的第一个指数位被忽略后 ) 进入指数。但我不确定这是否是正确的方法,它在纸面上效果很好。
编辑3: 我卡在单声道版本 2.0.50727.1433
检查您的单声道版本,因为 ConvertToInt
和 ConvertToIntTruncated
have been added 4 years ago 应该在 2.10 版本后就存在。
您可以将问题简化为取一个浮点数 mod 1. 为简化这一点,您可以使用位运算计算浮点数的底数,然后使用浮点数减法。以下是这些操作的(不安全)C# 代码:
// domain is assumed to be positive
// returns value in [0,domain)
public float fmodulus(float val, float domain)
{
if (val < 0)
{
float negative = fmodulus(-val, domain);
if (domain - negative == domain)
return 0;
else
return domain-negative;
}
if (val < domain)
return val; // this avoids losing accuracy
return fmodOne(val / domain) * domain;
}
// assumes val >= 1, so val is positive and the exponent is at least 0
unsafe public float fmodOne(float val)
{
int iVal = *(int*)&val;
int uncenteredExponent = iVal >> 23;
int exponent = uncenteredExponent - 127; // 127 corresponds to 2^0 times the mantissa
if (exponent >= 23)
return 0; // not enough precision to distinguish val from an integer
int unneededBits = 23 - exponent; // between 0 and 23
int iFloorVal = (iVal >> unneededBits) << unneededBits; // equivalent to using a mask to zero the bottom bits of the mantissa
float floorVal = *(float*)&iFloorVal; // convert the bit pattern back to a float
return val-floorVal;
}
例如,fmodulus(100.1f, 1) 为 0.09999847。 100.1f 的位模式是
0 10000101 10010000011001100110011
floorVal (100f) 的位模式是
0 10000101 10010000000000000000000
浮点数减法给出接近 0.1f 的值:
0 01111011 10011001100110000000000
实际上,我很惊讶最后8位被清零了。我认为只有 0.1f 的最后 6 位应该被替换为 0。也许可以比依赖浮点减法做得更好。