仅使用浮点乘法和除法对double进行按位转换

用汤姆·莱勒（Tom Lehrer）的话说，“这完全没有意义，但也许有一天可能会在某些奇怪的情况下对某些人有用。”

问题如下：假设您正在一个仅提供IEEE-754双精度浮点（“ double”）类型的编程环境中工作，并且没有可以访问该类型表示形式的操作（例如C ++按位转换操作） ，或Javascript的DataView对象）。您有一个double，并且想要将其转换为两个无符号的32位整数（存储为double）的位表示形式，反之亦然。这个问题有时会出现，但是我很好奇另一个问题：您的编程环境有多严格？您能只用浮点乘法和加法来做到吗？

使用浮点运算的按位转换在有限的解释语言或C ++ constexpr上下文等情况下很有用。通常，可以使用二进制搜索来完成从double到int的转换，将其与2的幂进行比较以找出指数的位。可以通过更多二进制搜索或使用指数知识将分数位缩放到整数范围，从中提取分数位。

但是，可以在没有按位运算，分支，求幂，除法或浮点比较的情况下完成此操作吗？

乍一看似乎不太可能，但是我发现答案是肯定的，尽管有一些值得注意的警告，但乘法和加法大多数就足够了。即使没有这些限制，在大多数环境中也无法区分或生成不同的NaN值（无逐位转换），但是仅使用乘法和加法就不可能将NaN，Infinity或-Infinity转换为无符号的32位值。另一个有问题的值是“负零”，使用加法和乘法无法将其与“正零”区分开。我的所有代码都使用减法，尽管可以通过用+（b * -1）替换a-b来删除。最后，这依赖于IEEE-754操作（在通常的舍入模式下，“四舍五入，并保持偶数”），因此它在使用不安全的数学优化（着色器编译器中的默认设置或启用）的环境中不起作用在许多其他编译器中通过/ fp：fast之类的标志来表示）。

因此，如果您只需要一个解决方案，就在这里，但否则请继续寻求解释：

function double_as_uint32s（double）{//不处理NaN，Infinity或-Infinity。将-0视为0。var a = double，b，c，d，e，f，g，h，i，j，k，l，m，n，低，高； f = 2.2250738585072014e-308 + a; j = 5e-324; b = j + f; b- = f; m = -5e-324； d = m + b; b = 4.4989137945431964e + 161； d = b * d; d = b * d; g = d * d; d = 1.0； g = d-g; h = m + f; f = h-f; f = j + f; f = b * f; f = b * f; f * = f; f = d-f; f * = g; g = -2.2250738585072014e-308 + a; h = j + g; h- = g; h = m + h; h = b * h; h = b * h; h * = h; h = d-h; c = m + g； c- = g; c = j + c; c = b * c; c = b * c; c * = c; c = d-c; c * = h; k = c * f; c = 5.562684646268003e-309 * a; g = j + c; g- = c; g = m + g； g = b * g; g = b * g; g * = g; g = d-g; h = m + c; h- = c; h = j + h; h = b * h; h = b * h; h * = h; h = d-h; g = h * g; h = a * g; g = d-g; c = g * c; g = 1024.0 * g； f = 2.0 + g； c + = h; h = 7.458340731200207e-155 * c； l = 1.0000000000000002; g ＝ l * h； g = m + g； e = j + g; e- = g; e = b * e; e = b * e; c = e * c; e = d-e; g = e * h; c = g + c; e = 512.0 * e; g = 8.636168555094445e-78 * c； e + = f; f ＝ l * g； f = m + f; h = j + f; f = h-f; f = b * f; f = b * f; c = f * c; f = d-f; g = f * g; f = 256.0 * f； c = g + c; e = f + e; f = 2.938735877055719e-39 * c; g = l * f; g = m + g； h = j + g; g = h-g; g = b * g; g = b * g; c = g * c; g = d-g; f = g * f; c = f + c; f = 128.0 * g； g = 5.421010862427522e-20 * c； e = f + e; f ＝ l * g； f = m + f; h = j + f; f = h-f; f = b * f; f = b * f; c = f * c; f = d-f; g = f * g; f = 64.0 * f； c = g + c; e = f + e; i = 2.3283064365386963e-10; f = i * c; g = l * f; g = m + g； h = j + g; g = h-g; g = b * g; g = b * g; c = g * c; g = d-g; f = g * f; c = f + c; f = 32.0 * g； g = 1.52587890625e-05 * c; e = f + e; f ＝ l * g； f = m + f; h = j + f; f = h-f; f = b * f; f = b * f; c = f * c; f = d-f; g = f * g; f = 16.0 * f； c = g + c; e = f + e; f = 0.00390625 * c； g = l * f; g = m + g； h = j + g; g = h-g; g = b * g; g = b * g; c = g * c; g = d-g; f = g * f; c = f + c; f = 8.0 * g； g = 0.0625 * c； e = f + e; f ＝ l * g； f = m + f; h = j + f; f = h-f; f = b * f; f = b * f; c = f * c; f = d-f; g = f * g; f = 4.0 * f; c = g + c; e = f + e; f ＝ 0.25 * c； g = l * f; g = m + g； h = j + g; g = h-g; g = b * g; g = b * g; c = g * c; g = d-g; f = g * f; c = f + c; f = g + g； e = f + e; n = 0.5； f = n * c; g = l * f; g = m + g； h = j + g; g = h-g; g = b * g; g = b * g; c = g * c; g = d-g; f = g * f; c = f + c; e = g + e; f = d-k; g = j + a; g- = a; g = m + g； g = b * g; g = b * g; g * = g; g = d-g; h = m + a; a = h-a; a = j + a; a = b * a; a = b * a; a * = a; a = d-a; a * = g; g = f * a; a = d-a; a = e * a; a + = g; e = l * c; e = m + e; g = j + e; e = g-e; e = b * e; e = b * e; g = n * c; c = e * c; e = d-e; e * = g; c = e + c; e = 4.450147717014403e-308 + c; g = j + e; g- = e; g = m + g； g = b * g; g = b * g; g * = g; g = d-g; h = m + e; e = h-e; e = j + e; e = b * e; e = b * e; e * = e; e = d-e; e * = g; g = e + e; d- = g; c = d * c; c = b * c; b * = c; c = -4503599627370496.0 * f; c + = b; b = i * c; b = -0.4999999998835847 + b; b = 4503599627370497.0 + b; d = -4503599627370497.0 + b; b = 2147483648.0 * e; a = 1048576.0 * a; a = b + a; b = d + a; a = -4294967296.0 * d; a + = c;低= a;高= b; return [low，high];}函数uint32s_as_double（low，high）{var a =低，b =高，c，d，e，f，g，h，i，j，k，l，m; b = 9.5367431640625e-07 * b; f = -0.4999999998835847; c = f + b; g = 4503599627370497.0； c = g + c; e = -4503599627370497.0; c = e + c; d = b-c； c = 0.00048828125 * c; b = f + c; b = g + b; k = e + b; l = c-k； j = 2.2250738585072014e-308; c = j + 1; c- = l; i = 4.49423283715579e + 307; b = i * c; c = 1.0； b = c-b; a = 2.220446049250313e-16 * a; h = -0.00048828125 + 1； a = d + a; d = b * h; d + = d; h = f + d; h = g + h； h = e + h; d- = h; b + = a; b = j * b; m = 1.3407807929942597e + 154; h = m * h; h = c + h; b = h * b; b * = h; d + = d; h = f + d; h = g + h； h = e + h; d- = h; h = m * h; h = c + h; b = h * b; d + = d; h = f + d; h = g + h； h = e + h; d- = h; h = 1.157920892373162e + 77 * h; h = c + h; b = h * b; d + = d; h = f + d; h = g + h； h = e + h; d- = h; h = 3.402823669209385e + 38 * h; h = c + h; b = h * b; d + = d; h = f + d; h = g + h； h = e + h; d- = h; h = 1.8446744073709552e + 19 * h; h = c + h; b = h * b; d + = d; h = f + d; h = g + h； h = e + h; d- = h; h = 4294967295.0 * h; h = c + h; b = h * b; d + = d; h = f + d; h = g + h； h = e + h; d- = h; h = 65535.0 * h; h = c + h; b = h * b; d + = d; h = f + d; h = g + h； h = e + h; d- = h; h = 255.0 * h; h = c + h; b = h * b; d + = d; h = f + d; h = g + h； h = e + h; d- = h; h = 15.0 * h; h = c + h; b = h * b; d + = d; f + = d; f = g + f; e + = f; d- = e; e = 3.0 * e; e = c + e; b = e * b; d + = d; d = c + d; b = d * b; d = -0.99951171875 + 1； e = j + d; d = e-d; d = i * d; e = j + a; a = e-a; a = i * a; a = c-a; a = d * a; a = m * a; a = m * a; a- = a; a = b + a; b = k + k; b = c-b; a * = b;返回一个；}

（我主要是在开玩笑，但是总有一天会发现代码被埋在库中，这很有趣，而且移植到几乎任何语言都应该很容易。）

这旨在对您需要了解的有关double的所有内容进行简要说明，以理解本文的其余部分。如果您已经知道，请跳过。

double是64位值。从最高有效位到最低有效位，它由一个1位符号，一个11位指数和一个52位分数组成。这些位将被解释为特殊值或数值，如以下伪代码中所述。伪代码中的运算和值具有无限的精度，并且**是幂运算。

if（sign == 1）s = -1; else s = 1; if（exponent == 0x7FFF）{//最大指数表示一个特殊值，如果（fraction == 0）返回NaN; //不是数字，否则（sign == 1）返回-Infinity;否则返回Infinity;}否则，如果（exponent == 0）{//零指数表示次标准值。 return s *（0.0 +分数*（2 ** -52））*（2 ** -1022）;} else {//其他都是正常值。 return s *（1.0 +分数*（2 ** -52））*（2 **（指数1023））;}

普通值具有隐式前导1，可以将其视为“ 1.fraction”。次正态不是，因此可以认为是“ 0.fraction”，但在其他方面它们与指数== 1相同。

隐式前导确保每个值都有唯一的表示形式（0 / -0和NaN除外）。

次规范确保可表示数字之间的距离只会随着您接近零而不断减小，因此两个连续值（也称为“最后一个单位”或ULP）之间的差异始终可以精确表示。

对于正数，浮点值以其64位整数表示形式增加，因此可以将它们作为整数进行比较，或者可以通过在int64表示形式中加1来找到下一个可表示的值。

双精度数的加法和乘法被定义为精确，无限精确的数学加法和乘法。如果确切的结果可以用一个double表示，则该double是结果，否则将进行舍入。 IEEE-754指定了几种可以使用的舍入模式，但我将重点介绍使用最广泛的一种“舍入到最近，并保持偶数”。这意味着将使用最接近的可表示值，或者如果精确结果位于两个可表示值之间的一半，则将使用其最低有效分数位为零的值。如果无限精确的结果太大或太小，则将其舍入为Infinity或-Infinity（有关正式定义，请参阅IEEE-754-2008第4.3.1节）。

最后，我们应该考虑特殊值。如果NaN是加法或乘法的输入，则结果将始终为NaN。与Infinity或-Infinity的乘法和加法将导致其他Infinity或-Infinity值，但将Infinity乘以零或从Infinity中减去Infinity除外，这两者都会导致NaN。

从这一点开始，这是对诸如识字编程之类的尝试，基本上按照我创建它的顺序进行介绍，从乘，加，减开始，然后逐步构建更强大的功能。该代码是用C ++编写的，并且经过重构以简化说明。我确实使用了循环和函数，但仅在编译器可以完全展开或内联的情况下使用。

我已省略了double p2（int e）函数，该函数提供2的幂–在使用该函数的任何地方，它都会内联为一个常量，但是最简单的方法是使用带有2098个值的查找表。

宏CONSTEXPR的定义如下，主要是为了允许调整内联或轻松地从所有内容中删除constexpr关键字：

在整篇文章中，我都使用指数来表示双精度的已编码指数位，而不是无偏/已解码指数（指数-1023）。希望那不会太令人困惑。

我首先研究仅加法和乘法可以做什么。假设“ true”为1.0，“ false”为0.0，我实现了一些逻辑运算：

CONSTEXPR double double_and（double a，double b）{返回a * b;} CONSTEXPR double double_not（double a）{return 1-a;} CONSTEXPR double double_or（double a，double b）{返回a + b-a * b ;} CONSTEXPR double select（double condition，double if_true，double if_false）{返回条件* if_true + double_not（condition）* if_false;}

这些内容大多在不进行进一步评论的情况下进行介绍，因为可以对其进行详尽的测试。但是select是一件棘手的事情。因为Infinity * 0 = NaN且NaN +任何事物= NaN，所以我们永远不能忽略Infinity值，并且对于从不执行可能创建它们的操作必须谨慎。

鉴于我想将任意浮点数转换为其按位表示，因此我必须首先弄清楚我可以在任何浮点数上执行哪些操作而不会冒创建无穷大的风险。

这里的一种选择是将数值乘以1.0到-1.0之间（含1.0和-1.0），因为结果永远不会增加。这适用于任何舍入模式。

我们还可以在p2（969）和-p2（969）互斥之间添加任何常量值，因为当添加到幅度最大的正值或负值时，该值不会四舍五入。但是，这仅在取整（round-to-nearest）或取整（round-toward-zero）模式下有效，因为当加上最小的非零值时，取整（round-toward-positive）和取整（round-toward-negative）可能取整为Infinity。

我认为我需要构造（x == y）和（x y时，结果将为正；当x = 1;测试/ = 2）{双重测试= tmp * p2（-test）; double too_small = is_exp_0_or_1（trial）; tmp =选择（too_small，tmp，试用）; e + =选择（too_small，0，测试）; } return select（is_exp_0_or_1（v），double_not（is_exp_0（v）），e + 2）;}

这将检查编码后的指数是否小于2 + 1024，如果不小于2 + 1024，则会从编码后的指数中减去1024（乘以p2（-1024）），然后将1024.0添加到我们的指数值中。以较小的2的幂进行重复，直到我们知道剩余的编码指数为0、1或2，并且e变量将包含相减的量。最后，它使用is_exp_0_or_1和is_exp_0函数显式处理零和一的情况。

这是迈向按位转换的重要一步，但get_encoded_exponent中的tmp很有趣。在函数结束时，我们保留了其符号位和分数位，但其指数仅转换为0、1或2。这使得测试x = 1;测试/ = 2）{双重测试= res * p2（-test）; res = select（is_exp_0_or_1（trial），res，trial）; } return select（is_exp_0_or_1（res），res，res * p2（-1））;}

现在，我们可以添加一个常数，以将所有非负值移出零或一的指数范围，这样只有小于零的值才能通过is_exp_0_or_1测试。

拼凑出按位投射之前，我们需要的最终工具是地板。为此，我们将仅考虑零到p2（52）之间的数字，并使用我过去见过的技巧（例如，在musl libc的floor.c中）。诀窍是加减p2（52）。在p2（52）到p2（53）的范围内，ULP恰好为1，因此x + p2（52）-p2（52）执行舍入到最接近整数运算。从这里，我们可以简单地检查它是否四舍五入，如果是，则减去1：

CONSTEXPR double small_positive_floor（double v）{//警告：对于负数和p2（52）上的某些值不正确//（尽管名称有效，但对零有效）double r = v + p2（52）-p2（52 ）;返回select（is_less_than（v，r），r-1，r）;}

这使我们可以从浮点整数中提取特定位。具体来说，我使用以下惯用法从整数x拆分n个低位：high_part = floor（x * p2（-n））; low_part = x-high_part * p2（n）;

那么，我们将double转换为bit的距离有多近？ get_encoded_exponent给我们提供了指数位。 is_less_than_zero给我们符号位。

对于小数，make_exp_0_or_1已为我们提供了所有小数位，但保留了符号，如果数字不是非正规数，则保留隐式前导1。

如果值为负，我们可以乘以-1来清除符号位。如果值不是非正规的，我们可以减去隐含的前导1，只剩下小数位，然后将其放大p2（1047），使小数1为1.0：

CONSTEXPR double get_fraction（double v）{double result = make_exp_0_or_1（v）* select（is_less_than_zero（v），-1，1）;结果-=选择（is_exp_0（v），0，p2（-1022））;结果=结果* p2（1074/2）* p2（1074/2）;返回结果；}

这为我们提供了一个1位符号值，一个11位指数值和一个52位分数值（所有值均以整数形式存储在double中），因此我们只需要将其拆分为两个32位值即可。

这些传统的按位运算是使用乘以2的幂乘以恒定移位（用底数来截断结果），再加上设置的位（而不是按位“或”

......