最近的一个问题是,编译器是否允许用浮点乘法替换浮点除法,这激发了我提出这个问题.
在严格的要求下,代码转换后的结果应与实际的除法运算在位上相同,可以看出,对于二进制IEEE-754算术,这对于2的幂的除数是可能的.只要除数的倒数可以表示,乘以除数的倒数就可以得到与除数相同的结果.例如,乘法0.5可以代替除法2.0.
然后人们想知道其他除数这样的替换是什么工作,假设我们允许任何短指令序列取代除法但运行速度明显更快,同时提供相同的结果.特别是除了普通乘法之外,还允许融合乘法 - 加法运算.在评论中,我指出了以下相关文件:
Nicolas Brisebarre,Jean-Michel Muller和Saurabh Kumar Raina.当事先知道除数时,加速正确舍入的浮点除法.IEEE Transactions on Computers,Vol.53,第8期,2004年8月,第1069-1072页.
本文作者提出的技术预先计算了除数y的倒数作为归一化的头尾对z h:z l如下:z h = 1/y,z l = fma(-y,z h,1 )/ y.之后,将除数q = x/y计算为q = fma(z h,x,z l*x).本文推导出除数y必须满足的各种条件才能使该算法起作用.正如人们容易观察到的那样,当头尾迹象不同时,该算法存在无穷大和零的问题.更重要的是,它将无法为数量非常小的股息x提供正确的结果,因为商尾的计算z l*x受到下溢的影响.
本文还提到了另一种基于FMA的划分算法,该算法由Peter Markstein在IBM工作时开创.相关参考是:
PW Markstein.在IBM RISC System/6000处理器上计算基本功能.IBM Journal of Research&Development,Vol.1990年1月34日第1号,第111-119页
在Markstein算法中,首先计算倒数rc,从中形成初始商q = x*rc.然后,用FMA精确地计算除法的余数r = fma(-y,q,x),并且最终计算出改进的,更准确的商,其为q …
rootn (float_t x, int_t n)是一个计算第n个根x 1/n的函数,并且受一些编程语言(如 OpenCL)的支持.当使用IEEE-754浮点数时n,假设只x需要处理归一化的操作数,则可以基于对基础位模式的简单操作生成任何有效的低精度起始近似值.
二进制指数root (x, n)将是二进制指数的1/n x.IEEE-754浮点数的指数字段是有偏差的.我们可以简单地将偏置指数除以n,然后应用偏移量来补偿先前忽略的偏差,而不是对指数进行不偏置,将其除去,并重新偏置结果.此外,我们可以简单地划分整个操作数x,重新解释为整数,而不是提取,然后除去指数字段.找到所需的偏移是微不足道的,因为1的参数将返回1的结果n.
如果我们有两个辅助函数可供我们使用,__int_as_float()它将IEEE-754重新解释binary32为int32,并且__float_as_int()重新解释int32操作数为binary32,则我们以rootn (x, n)直接的方式得出以下低精度近似:
rootn (x, n) ~= __int_as_float((int)(__float_as_int(1.0f)*(1.0-1.0/n)) + __float_as_int(x)/n)
通过恒定除数的整数除法的公知优化,__float_as_int (x) / n可以将整数除法简化为移位或乘法加移位.一些有用的例子是:
rootn (x, 2) ~= __int_as_float (0x1fc00000 + __float_as_int (x) / 2) // sqrt (x)
rootn (x, 3) ~= __int_as_float (0x2a555556 + …有一个现有的问题“3 个长整数的平均值”，它特别关注三个有符号整数的平均值的有效计算。
然而，无符号整数的使用允许额外的优化不适用于上一个问题中涵盖的场景。这个问题是关于三个无符号整数的平均值的有效计算，其中平均值向零舍入，即在数学术语中我想计算？(a + b + c) / 3 ?。
计算此平均值的一种直接方法是
avg = a / 3 + b / 3 + c / 3 + (a % 3 + b % 3 + c % 3) / 3;
首先，现代优化编译器会将除法转换为具有倒数加移位的乘法，并将模运算转换为反向乘法和减法，其中反向乘法可能使用许多体系结构上可用的scale_add习语，例如leax86_64的，add用lsl #n在ARM，iscadd在NVIDIA GPU。
在尝试以适用于许多常见平台的通用方式优化上述内容时，我观察到整数运算的成本通常在逻辑关系中？（添加|子）？转移？规模_添加？MUL。这里的成本是指所有延迟、吞吐量限制和功耗。当处理的整数类型比本地寄存器宽度更宽时，例如在uint64_t32 位处理器上处理数据时，任何此类差异都会变得更加明显。
因此，我的优化策略是尽量减少指令数量，并在可能的情况下用“廉价”操作替换“昂贵”操作，同时不增加寄存器压力并为宽无序处理器保留可利用的并行性。
第一个观察结果是，我们可以通过首先应用产生和值和进位值的 CSA（进位保存加法器）将三个操作数的和减少为两个操作数的和，其中进位值的权重是和的两倍价值。在大多数处理器上，基于软件的 CSA 的成本是五个逻辑s。一些处理器，如 …
c algorithm bit-manipulation micro-optimization extended-precision
为了以合理的精度简单有效地实现快速数学函数,多项式极小极大近似通常是选择的方法.Minimax近似通常使用Remez算法的变体生成.各种广泛使用的工具,如Maple和Mathematica,都具有内置功能.通常使用高精度算法计算生成的系数.众所周知,简单地将这些系数舍入到机器精度导致所得实施中的次优精度.
相反,人们搜索密切相关的系数集合,这些系数可以精确地表示为机器编号,以生成机器优化的近似值.两篇相关论文是:
Nicolas Brisebarre,Jean-Michel Muller和Arnaud Tisserand,"计算机高效多项式近似",ACM数学软件交易,Vol.32,第2期,2006年6月,第236-256页.
Nicolas Brisebarre和Sylvain Chevillard,"高效多项式L∞近似",第18届IEEE计算机算术研讨会(ARITH-18),蒙彼利埃(法国),2007年6月,第169-176页.
来自后一篇论文的LLL算法的实现可作为Sollya工具的fpminimax()命令获得.我的理解是,所有用于生成机器优化近似的算法都是基于启发式算法,因此通常不知道通过最佳近似可以实现什么样的精度.我不清楚用于评估近似的FMA(融合乘法 - 加法)的可用性是否对该问题的答案有影响.在我看来它应该是天真的.
我目前正在研究[-1,1]上的反正切的简单多项式近似,使用Horner方案和FMA在IEEE-754单精度算法中进行评估.请参阅atan_poly()下面的C99代码中的功能.由于目前无法访问Linux机器,我没有使用Sollya来生成这些系数,而是使用了我自己的启发式算法,可以将其简单地描述为最陡的体面和模拟退火的混合(以避免陷入局部最小值) .机器优化多项式的最大误差非常接近1 ulp,但理想情况下我希望最大ulp误差低于1 ulp.
我知道我可以改变我的计算以提高精度,例如使用表示超过单精度精度的前导系数,但我想保持代码完全一样(即尽可能简单)仅调整系数以提供最准确的结果.
"经过验证的"最佳系数集将是理想的,欢迎指向相关文献的指针.我进行了一次文献检索,但找不到任何有助于超越Sollya的艺术发展水平的论文,也没有一篇能够fpminimax()在这个问题上研究FMA(如果有的话)的作用.
// max ulp err = 1.03143
float atan_poly (float a)
{
float r, s;
s = a * a;
r = 0x1.7ed1ccp-9f;
r = fmaf (r, s, -0x1.0c2c08p-6f);
r = fmaf (r, s, 0x1.61fdd0p-5f);
r = fmaf (r, s, -0x1.3556b2p-4f);
r = fmaf (r, s, 0x1.b4e128p-4f);
r = fmaf (r, s, -0x1.230ad2p-3f);
r = …在各种情况下,例如对于数学函数的参数减少,需要计算(a - K) / (a + K),其中a是正变量参数并且K是常数.在许多情况下,K是2的幂,这是与我的工作相关的用例.我正在寻找比直接划分更准确地计算这个商的有效方法.可以假设对融合乘法 - 加法(FMA)的硬件支持,因为此操作由所有主要的CPU和GPU架构提供,并且可以通过函数fma()和C/C++获得fmaf().
为了便于探索,我正在尝试float算术.由于我计划将该方法移植到double算术,因此不能使用高于参数和结果的本机精度的操作.到目前为止我的最佳解决方案是
/* Compute q = (a - K) / (a + K) with improved accuracy. Variant 1 */
m = a - K;
p = a + K;
r = 1.0f / p;
q = m * r;
t = fmaf (q, -2.0f*K, m);
e = fmaf (q, -m, t);
q = fmaf (r, e, …我试图找到最有效的方法来计算 32 位无符号整数的模 255。我的主要重点是找到一种可以在 x86 和 ARM 平台上运行良好的算法，并着眼于除此之外的适用性。首先，我试图避免内存操作（这可能很昂贵），所以我正在寻找有点复杂的方法，同时避免使用表格。我还试图避免可能昂贵的操作，例如分支和乘法，并尽量减少使用的操作和寄存器的数量。
下面的 ISO-C99 代码捕获了我迄今为止尝试过的八个变体。它包括一个用于详尽测试的框架。我对这个粗略的执行时间测量进行了猛烈抨击，这似乎工作得很好，可以获得第一次性能印象。在一些平台上我试过（全部具有快速整数倍）的变种WARREN_MUL_SHR_2，WARREN_MUL_SHR_1和DIGIT_SUM_CARRY_OUT_1似乎是最高效的。我的实验表明，我在Compiler Explorer 中尝试的 x86、ARM、PowerPC 和 MIPS 编译器都很好地利用了特定于平台的功能，例如三输入LEA、字节扩展指令、乘法累加和指令预测。
该变体NAIVE_USING_DIV使用整数除法，与除数反乘，然后减法。这是基线情况。现代编译器知道如何有效地实现 255 的无符号整数除法（通过乘法），并将在适当的情况下使用离散替换反乘。要计算模数，base-1可以对base数字求和，然后折叠结果。例如3334 mod 9: sum 3+3+3+4 = 13, fold 1+3 = 4. 如果折叠后的结果是base-1，我们需要生成0来代替。DIGIT_SUM_THEN_FOLD使用这种方法。
A. Cockburn，“使用 8/16 位算法有效实现 OSI 传输协议校验和算法”，ACM SIGCOMM 计算机通信评论，卷。17, No. 3, 七月/八月 1987 年，第 13-20 页
展示了base-1在校验和计算模 255 的上下文中有效地添加数字模数的不同方法。计算数字的逐字节总和，并且在每次添加之后，也添加来自加法的任何进位。所以这将是一个ADD a, b …
我目前正在研究如何使用各种现代处理器的快速单精度浮点倒数功能来计算基于定点Newton-Raphson迭代的64位无符号整数除法的起始近似.它需要尽可能精确地计算2 64 /除数,其中初始近似必须小于或等于数学结果,基于以下定点迭代的要求.这意味着这种计算需要低估.我目前有以下代码,基于广泛的测试,效果很好:
#include <stdint.h> // import uint64_t
#include <math.h> // import nextafterf()
uint64_t divisor, recip;
float r, s, t;
t = uint64_to_float_ru (divisor); // ensure t >= divisor
r = 1.0f / t;
s = 0x1.0p64f * nextafterf (r, 0.0f);
recip = (uint64_t)s; // underestimate of 2**64 / divisor
虽然这段代码很实用,但在大多数平台上并不是很快.一个显而易见的改进,需要一些特定于机器的代码,是r = 1.0f / t用代码来替换除法,该代码利用硬件提供的快速浮点倒数.这可以通过迭代来增强,以产生在数学结果的1 ulp内的结果,因此在现有代码的上下文中产生低估.x86_64的示例实现将是:
#include <xmmintrin.h>
/* Compute 1.0f/a almost correctly rounded. Halley iteration with cubic convergence */
inline float fast_recip_f32 (float a)
{
__m128 t; …我需要一个高效的算法在两个浮点数之间做数学::幂函数,你知道怎么做,(我需要算法不使用函数本身)
我正在研究基于双输入最小/最大操作的九个元素的排序和中值选择网络.Knuth,TAOCP Vol.3,第2版.状态(第226页)九元素排序网络需要至少25次比较,这转换为相同数量的SWAP()基元或50分钟/最大值操作.显然,通过消除冗余操作,可以将分拣网络转换为中值选择网络.传统观点似乎是,这不会导致最佳的中值选择网络.虽然这似乎在经验上是正确的,但我在文献中找不到证据证明这一定是必然的.
LukáŝSekanina,"中位电路的进化设计空间探索".在:EvoWorkshops,2004年3月,第240-249页,给出了最佳九输入中值选择网络所需的最小/最大操作次数为30(表1).我通过John L. Smith给出的众所周知的中值选择网络"在XC4000E FPGA中实现中值滤波器"来验证这一点.XCELL杂志,Vol.23,1996,p.来自Chaitali Chakrabarti和Li-Yu Wang早期工作的"9"中间网络,"用于排序过滤器的基于网络的新型排序过滤器".IEEE超大规模集成系统交易,Vol.2,No.4(1994),pp.502-507,其中后者通过简单地消除冗余分量转换成前者.请参阅以下代码中的变体4和5.
通过消除冗余操作,检查公布的最佳九元排序网络是否适合转换为有效的中间选择网络,我设法找到的最佳版本来自John M. Gamble的在线生成器,它需要32分钟/最大操作,所以两个害羞的最佳操作计数.这在下面的代码中显示为变体1.其他最佳分拣网络分别减少到36分钟/最大操作(变体2)和38分钟/最大操作(变体3).
是否有任何已知的九元素分拣网络(即50个双输入最小/最大操作)通过单独消除冗余操作,减少到最佳九输入中值选择网络(具有30个双输入最小/最大操作) ？
下面的代码使用float数据作为测试用例,因为许多处理器为浮点数据提供最小/最大操作,但不提供整数数据,GPU是一个例外.由于特殊浮点操作数的问题(在我的实际用例中没有出现),最佳代码序列通常需要使用编译器提供的"快速数学"模式,例如在Godbolt测试平台中.
#include <cstdlib>
#include <cstdio>
#include <algorithm>
#define VARIANT 1
#define FULL_SORT 0
typedef float T;
#define MIN(a,b) std::min(a,b)
#define MAX(a,b) std::max(a,b)
#define SWAP(i,j) do { T s = MIN(a##i,a##j); T t = MAX(a##i,a##j); a##i = s; a##j = t; } while (0)
#define MIN3(x,y,z) MIN(a##x,MIN(a##y,a##z))
#define MAX3(x,y,z) MAX(a##x,MAX(a##y,a##z))
#define MED3(x,y,z) …在JavaScript中,Math.cbrt(1728)评估确切的结果12.
然而,看似等效的表达式Math.pow(1728, 1/3)评估为11.999999999999998.
为什么这些结果的精度会有所不同？