hon*_*ann 16 algorithm assembly x86-64 mathematical-optimization divide
我在X86-64汇编语言编写的代码库提供所有传统的按位,移位,逻辑,比较,运算和数学函数s0128,s0256,s0512,s1024签署整数类型和f0128,f0256,f0512,f1024浮点类型.到目前为止,我正在处理有符号整数类型,因为浮点函数可能会调用为整数类型编写的一些内部例程.
到目前为止,我已经编写并测试了函数来执行各种一元运算符,比较运算符以及加,减和乘法运算符.
现在我正在尝试决定如何为除法运算符实现函数.
我的第一个想法是,"Newton-Raphson必须是最好的方法".为什么？因为它给出了一个很好的种子(开始猜测)它很快收敛,我想我应该能够弄清楚如何在操作数上执行本机64位除法指令以获得一个很好的种子值.实际上,如果种子值精确到64位,要获得正确的答案应该只采取:
`s0128` : 1~2 iterations : (or 1 iteration plus 1~2 "test subtracts")
`s0256` : 2~3 iterations : (or 2 iterations plus 1~2 "test subtracts")
`s0512` : 3~4 iterations : (or 3 iterations plus 1~2 "test subtracts")
`s1024` : 4~5 iterations : (or 4 iterations plus 1~2 "test subtracts")
然而,更多地思考这个问题让我很奇怪.例如,我记得我编写的核心例程,它对所有大整数类型执行乘法运算:
s0128 : 4 iterations == 4 (128-bit = 64-bit * 64-bit) multiplies + 12 adds
s0256 : 16 iterations == 16 (128-bit = 64-bit * 64-bit) multiplies + 48 adds
s0512 : 64 iterations == 64 (128-bit = 64-bit * 64-bit) multiplies + 192 adds
s1024 : 256 iterations == 256 (128-bit = 64-bit * 64-bit) multiplies + 768 adds
尽管循环相当短且有效(包括缓存方式),但更广泛的数据类型的操作增长是相当可观的.此循环仅将结果的每个64位部分写入一次,并且永远不会读回结果的任何部分以进行进一步处理.
这让我想到更传统的移位和减法类型划分算法是否会更快,特别是对于较大类型.
我的第一个想法是:
result = dividend / divisor // if I remember my terminology
remainder = dividend - (result * divisor) // or something along these lines
#1:为了计算每个比特,通常如果除数小于或等于被除数,则从被除数中减去除数.好吧,通常我们只能通过检查其最重要的64位部分来确定除数肯定小于或绝对大于红利.只有当那些ms64位部分相等时,例程才必须检查下一个较低的64位部分,并且只有当它们相等时,我们才能检查更低的部分,依此类推.因此,在几乎每次迭代(计算结果的每一位)时,我们都可以大大减少为计算此测试而执行的指令.
#2:然而......平均而言,我们会发现大约50%的时间我们需要从被除数中减去除数,所以无论如何我们都需要减去它们的整个宽度.在这种情况下,我们实际执行了比传统方法更多的指令(我们首先减去它们,然后测试标志以确定除数<=被除数).因此,我们实现节约的一半时间,以及我们实现损失的一半时间.对于大型类型s1024(包含-16-64位组件),节省了大量资金并且损失很小,因此这种方法应该实现大量的净节省.在像s0128(包含-2- 64位组件)这样的微小类型上,节省很小,损失很大但不是很大.
所以,我的问题是,"什么是最有效的划分算法",给出:
#1: modern x86-64 CPUs like FX-8350
#2: executing in 64-bit mode only (no 32-bit)
#3: implementation entirely in assembly-language
#4: 128-bit to 1024-bit integer operands (nominally signed, but...)
注意:我的猜测是,实际的实现只能在无符号整数上运行.在乘法的情况下,结果是更容易和更有效(可能)将负操作数转换为正,然后执行无符号乘法,然后如果只有一个原始操作数为负则否定结果.但是,如果它是有效的,我会考虑一个有符号整数算法.
注:如果最佳答案是我的浮点类型(不同f0128,f0256,f0512,f1024),请解释原因.
注意:我的内部核心无符号乘法例程对所有这些整数数据类型执行乘法运算,产生双倍宽度结果.换一种说法:
u0256 = u0128 * u0128 // cannot overflow
u0512 = u0256 * u0256 // cannot overflow
u1024 = u0512 * u0512 // cannot overflow
u2048 = u1024 * u1024 // cannot overflow
我的代码库为每个有符号整数数据类型提供了两个版本的multiply:
s0128 = s0128 * s0128 // can overflow (result not fit in s0128)
s0256 = s0256 * s0256 // can overflow (result not fit in s0256)
s0512 = s0512 * s0512 // can overflow (result not fit in s0512)
s1024 = s1024 * s1024 // can overflow (result not fit in s1024)
s0256 = s0128 * s0128 // cannot overflow
s0512 = s0256 * s0256 // cannot overflow
s1024 = s0512 * s0512 // cannot overflow
s2048 = s1024 * s1024 // cannot overflow
这与我的代码库中的"永不丢失精度"和"永不溢出"的策略一致(当由于精度损失或由于溢出/下溢而导致答案无效时,将返回错误).但是,当调用双宽度返回值函数时,不会发生此类错误.
您肯定知道现有的任意精度包（例如http://gmplib.org/）以及它们是如何运行的吗？它们通常被设计为“尽可能快地”以任意精度运行。
如果您将它们专门用于固定大小（例如，应用[手动]部分评估技术来折叠常量和展开循环），我希望您获得非常好的例程，用于您想要的那种特定的固定大小精度。
此外，如果您还没有看过，请查看 D. Knuth 的Seminumerical Algorithms和oldie but really goodie，它提供了多精度算术的关键算法。（大多数软件包都基于这些想法，但 Knuth 有很好的解释并且非常正确）。
关键思想是将多精度数字视为非常大的基数（例如，基数 2^64），并将标准的三年级算术应用于“数字”（例如 64 位字）。除法包括“估计商（大基数）数字，将估计值乘以除数，从被除数中减去，左移一位，重复”，直到你得到足够多的数字来满足你。对于除法，是的，它都是未签名的（在包装器中进行符号处理）。基本技巧是很好地估计商数（使用处理器为您提供的单精度指令），并快速进行多精度乘以个位数。有关详细信息，请参阅 Knuth。请参阅有关多精度算术的技术研究论文（您可以进行一些研究）以获取奇异的（“最快的”）改进。