专利名称：使用单指令多数据(simd)运算来计算超越函数的制作方法
背景技术：
本发明涉及超越函数的计算。在许多领域中非常需要诸如指数、对数和三角函数及它们的反函数之类的超越函数的快速而准确的求值。为了更快的求值，软件实现在计算中通常使用查找表来逼近一个或多个中间值。
例如，实现浮点数学函数的标准方法是使用预先计算的值表并用基于表条目和较小的“被归约的”自变量的简单重构公式在它们之间进行内插。例如，浮点数x的正弦(sin)(x)可以使用下列重构公式用预先计算的各个“断点”的sin(A)和余弦(cos)(A)值表来计算sin(x)＝sin(A)+sin(A)[cos(r)-1]+cos(A)sin(r) [1]其中r＝x-A。通常，断点均匀间隔一定距离d(例如，对于sin是π/32)，因此对于n∈□，A＝nd。在断点间隔距离d的情况下，直接的余项运算能找到满足|r|≤d/2的被归约的自变量。如果此边界相当小，例如在2-5的数量级，则可通过多项式来逼近sin(r)和cos(r)-1，从而收敛迅速且不需要多项式有许多项，并且与总的结果的大小相比，该多项式的大小较小。
后一特性意味着与总的结果相比，多项式中的舍入误差相对较小，总的结果是由单个表条目(在上述例子中为sin(A))主导的。因此，计算可以被组织成表条目和相对较小项的最终相加，这使得总误差接近0.5个理想最小单位(ulp)。
在浮点超越函数的许多应用中，通常同时需要sin(x)和cos(x)。虽然提供能在单独的计算中有效率地计算这两者的组合sincos例程是合乎需要的，但上述的表驱动的技术引起严重的问题。因为当A较小时(例如当断点为最小的非零值±d且r≈±d/2时)，表条目主导的属性趋向于崩溃，所以将执行使用最初几个表条目的较小的输入的单独路径指令。此单独路径通常是纯多项式，且常常相当长，因为求值是对远远大于d/2的x来求值的。
在两条路径之间有分支要选择是相当不利的，因为难以通过重叠多个调用来实现软件流水线处理，并且还可引起严重的误预测惩罚。更严重的是，对于sin和cos的组合的单指令多数据(SIMD)实现困难将加剧，因为在两种情形中对于不同种类的值运用特殊分支。对于sin，它发生在输入接近π/2的偶数倍时，而对于cos，它发生在输入接近π/2的奇数倍时。因此，特别是在SIMD实现中，需要计算超越函数的无分支方式。
附图简要说明
图1为根据本发明的一个实施例的方法的流程图。
图2为根据本发明的一个实施例的确定sin(x)和cos(x)的方法的流程图。
图3为可配合本发明的实施例使用的计算机系统的框图。
详细说明可能需要大约同时为相同的x计算诸如sin(x)和cos(x)之类的浮点超越函数。在各种实施例中，能以与单个正弦或余弦的计算几乎相同的效率来一起计算正弦和余弦。
在某些实现中，可以使用SIMD浮点运算。在某些此类实现中，可以使用包括对压缩数据格式的运算并提供提高的SIMD计算性能的SIMD流扩展2(SSE2)指令。这些指令可以是IntelPENTIUM 4(英特尔奔腾4)处理器指令集或其它此类处理器指令集的一部分。
以这种方式，可以使用同一指令流分别在并行操作的一半中计算sin和cos。为了维持此并行性，根据本发明的一个实施例的算法可以使用“无分支”技术来避免要为小自变量提供专用代码，不然它会在sin和cos指令流之间产生不对称。结果，可以减少分支误预测。
在本发明的各种实施例中，可以用三个基本步骤来计算超越函数归约、逼近和重构。归约可用于根据预定等式来变换输入自变量x以将其限制于预定范围。接着，逼近是通过计算该归约的被归约的自变量的逼近多项式来执行。最后，重构使用该逼近多项式的结果和多项式余项来得到原始函数的最终结果。
现参见图1，所示的是根据本发明的一个实施例的方法的流程图。如图1所示，方法10始于归约给定函数的输入自变量x(框20)。在一个实施例中，归约可以取r＝x-A的形式。接着，可以用具有主导项f(A)+σr的多项式来逼近已被归约的自变量(框30)。在各种实施例中，不论输入自变量的大小如何，这两项总是主导最终结果。最终，可以通过对逼近结果和多项式余项求和来执行重构以得到最终结果(框40)。
本发明的实施例可适用于在x＝0附近斜率大小接近2的幂的数学函数f(x)。此类函数包括例如在x＝0处均具有接近1的斜率的sin(x)和正切(tan)(x)，并通过使用cos(x)＝sin(x+π/2)而包括cos(x)。
在这些实施例中，可以执行归约来得到用于计算逼近的范围被归约的自变量。在一个实施例中，逼近可以表示成 其中，对于某个α，|o|＝±2α。尽管α可以变化，但是在某些实施例中它可以在大约-3和1之间，并且在特定的实施例中可以在大约1/8和1之间。在上述公式2中，f(A)和f′(A)可以从查找表中合适的断点得到。在某些实施例中，α可以在x的范围上变化，且可以制成类似于f(A)的查找表的形式的表格。
作为一个例子，对于正弦函数，核心逼近可以采用下列形式sin(x)＝(sin(A)+σr)+(cos(A)-σ)□r+sin(A)[cos(r)-1]+cos(A)[sin(r)-r][3]其中，σ为舍入到1位精度的cos(A)。sin(A)和cos(A)都可以通过找到存储在查找表中的合适的断点来获得。其中A非常小，σ＝±1。在其它实施例中，σ可以等于最接近的2的幂。
此逼近的重构具有以下特性即使是对于很小的x，最前面的两项f(A)+σr(在上述例子中，sin(A)+σr)总是构成最终答案的主导部分。在多项式的低端|(f′(A)-σ)·r|远远小于|σr|，而在高端，f(A)大到足以主导该重构。
因为乘以2的幂是准确的，所以总是可通过简单的浮点乘法来准确地计算±σr。f(A)+σr的和则可以通过准确求和的技术分两部分来计算。因为通常或者f(A)＝0，或者|σr|≤|f(A)|，所以可以通过进行下列三个连续的加/减运算来获得准确的和Hi＝f(A)+σr [4]med＝Hi-f(A) [5]Lo＝σr-Med[6]这些运算准确地产生Hi+Lo＝f(A)+σr，且Hi用作总的结果的高部分，而Lo可以被加入多项式和其它部分中。虽然上述求和需要几次浮点运算，但是其等待时间通常大大低于完全多项式的等待时间，因此，对总的等待时间具有最小的影响。
在一个特定实施例中，上述一般方法可以理想地适用于sin和cos的组合实现。在这一实施例中，除了异常小或异常大的输入的非常罕见的情形以外，除单个常数以外算法的两“侧”可以完全相同。现在参见图2，图2示出根据本发明的一个实施例的确定sin(x)和cos(x)的方法的流程图。如图2所示，方法100始于接收对sin(x)和cos(x)的请求(框110)。例如，在某些实施例中，未编译的程序可包括执行sin(x)和cos(x)的计算的函数调用。在编译期间，编译器可使函数调用被对这里讨论的组合sincos运算的函数调用所取代，因为该程序很可能会在对sin(x)的函数调用附近的代码中包括对cos(x)的函数调用。
仍参见图2，接着可以执行x的归约，例如，r＝x-A(框120)。然后，可以根据多项式逼近来并行地逼近sin(A)和sin(A+π/2)使得f(A)+σr为该逼近的两个主导项(框130)。最后，可以通过用逼近结果和多项式余项的求和来并行地重构sin(x)和cos(x)。以这种方式，可以在与获取sin(x)或cos(x)所需的时间量基本相同的时间量中获得sin(x)和cos(x)(框140)。另外，这些结果可以通过利用使用SIMD指令的指令级并行性以无分支的方式取得。
因此，根据方法100的流程图，从x至r的初始范围归约可以如下执行x≈Nπ32+r---(7)]]>因此，|r|≤π64+TM,]]>其中TM为机器的单位舍入，例如，对于单精度是2-24或对于双精度是2-53。在此特定实施例中，输入可限于在|N|≤932560的情况下的输入，因为在此以外，范围归约可能不够精确。因此，如果输入超过该值，可以使用具有更精确的范围归约的替换算法。然而，应理解在通常的应用中预期这些值不常出现。
另外，在此特定实施例中，在所产生的近似为x4/7！的最小中间结果在双精度下可能下溢的情况下的输入也可能由此对|x|≤2-252引起走向专用代码的分支。可以通过查看输入的指数和最高几位有效位来测试很小及很大自变量的不测事件。因此，对于2-252≤|x|≤90112可以取主路径，它基本上可涵盖所有这些输入。
然而，对于异常输入，放弃和使用替换算法是唯一需要的分支。根据此特定实施例的下列算法是无分支的，并且可以按需要计算正弦和余弦。虽然这里讨论的算法是就正弦而给出的，但是也可以通过将N加上16(即，x加上 )来得到余弦。
为了避免分支，每次可以最高精度地执行范围归约r＝x-N(P1+P2+P3)[8]其中，P1和P2为32位的数(所以乘以N是精确的)而P3为53位的数，每个数都是表示π/32的值的机器数。这些近似的π一起足以应付受限制的范围内的所有情形。在此特定实施例的其它实现中，执行下列两个步骤r＝x-N(P1+P2)[9]上式为多项式计算给出足够好的r，并且甚至简单的x-NP1做最高项也已足够。因此，可以隐藏部分归约的等待时间。
对于根据此特定实施例的算法，主归约序列为·y=32πx]]>·N＝integer(y)·m1＝NP1m2＝NP2
·r1＝x-m1·r＝r1-m2(它可用于大部分计算)·c1＝r1-rm3＝NP3·c2＝c1-m2·c＝c2m3可以用“移位器”法来舍入到整数，即，N＝(y+s)-s，其中，s＝252+251。
接着，使用范围被归约的值，可以根据B＝M{π/32}用查表来逼近sin(B)，其中M＝N mod64(注意，为了将此讨论与上述一般实施例相关，B＝A)。在此特定实施例中，所存储的值为σ，它是最接近cos(B)的2的幂；Chl，它是cos(B)-σ的53位的值；以及Shi和Slo，它们分别是sin(B)的(53和24)位的值。
所存储的这些值可以被组织成4*64双精度的数。即，可以在64个断点处计算每个值(例如，Nπ/64，其中N＝1到64)。然而，Slo和σ均可表示成单精度数，所以在某些实施例中，这些值可被存储为3*64个双精度的数。
核心逼近的多项式可以如下组织sin(B+r+c)＝[sin(B)+σr]+r(cos(B)-σ)+sin(B)[cos(r+c)-1]+cos(B)[sin(r+c)-r] [10]该式近似为[Shi+σr]+Chlr+Slo+Shi[(cos(r)-1)-rc]+(Chl+σ)[sin(r)-r+c] [11]实际所计算的可以是此多项式逼近。和可以分成四个部分hi+med+pols+corr，其中，hi＝Shi+σr[12]med＝Chlrpols＝Shi(cos(r)-1)+(Chl+σ)(sin(r)-r)[13]corr＝Slo+c□((Chl+σ)-Shl□r) [14]
应注意，与最终结果相比，pols和corr非常小，而乘以σ是精确的，因为它是2的幂。因此，假设对各分量求和是精确的，只有med中有实质性误差，该误差由Chl中定标的逼近误差和乘法中的舍入误差构成。然而，Chl·r占最终结果的比例不大，因为此项中的误差在最终结果中从未超过约0.02ulp。
然而，在对各分量求和时应避免舍入误差，因为它们可能会对最终误差产生实质性的影响。通常，σr相对于Shi可能非常大；对于B＝{π\32}且r≈-π/64，有σr≈B/2。因此，Shi不是结果的主导部分，并且必须精确地进行Shi+σr求和。
实际上，等待时间临界部分是多项式计算，因此，在其被计算时，可以执行两次连续的补偿求和，即，Shl+σr的第一次相加，以及其高部分和Chlr的下一次相加。在一些实施例中，后者不是必需的，但可能是适合的，因为它显著提高准确度而不明显影响总的等待时间。事实上，在某些实施例中，这种扩展的精度和并行性一起提高了逼近的性能，因为多项式的求值顺序变得不重要。当能以任意顺序来对多项式求值时，就可以充分地利用并行性，从而，甚至长多项式也可以用最小等待时间来求值。
当A变大时，不再需要如此介意f′(A)-σ应很接近2的幂。在这一实施例中，可以使用σ＝0。或者，当A很大并可接受σr中的舍入误差时，可以用标准长度的浮点数替换σ。
在其它实施例中，如果已知r不具有完整个数的有效位，则可以使用更多位(例如2位或3位)而不是1位的σ的逼近而不会在乘积σr中引起舍入误差。如果通过典型的余项运算来计算r，则可能出现这种情形。例如，如果r＝x-Nd′被设置，其中 且d′为设计成允许精确地乘以N的d的短版本，则随着N增大，r中的有效位将减少。因此，在更远离0时，σ中的有效位数可能增加，这极佳地补偿了f′(A)不能再被2的幂很好地逼近的事实。
实施例可以在代码中实现，并可以被存储在其上已存储有指令的存储介质上，这些指令能用于将计算机系统编程以执行这些指令。该存储介质可包括但不限于任何类型的盘片，包括软盘、光盘、光盘只读存储器(CD-ROM)、可重写光盘(CD-RW)和磁光盘；半导体器件，例如只读存储器(ROM)、随机存取存储器(RAM)、可擦除可编程只读存储器(EPROM)、闪存、电可擦除可编程只读存储器(EEPROM)；磁或光卡或任何类型的适合存储电子指令的介质。
示例性实施例可以在用于由用硬件设备的合适组合配置的合适的计算机系统执行的软件中实现。图3为可以配合本发明的实施例使用的计算机系统400的框图。
现参见图3，在一个实施例中，计算机系统400包括处理器410，该处理器可包括通用或专用处理器，例如微处理器、微控制器、可编程门阵列(PGA)等。如这里所使用的，“计算机系统”一词可指任何类型的基于处理器的系统，例如，台式计算机、服务器计算机、膝上型计算机等。
在一个实施例中，处理器410可以通过主机总线415与存储器集线器430耦合，该存储器集线器可以通过存储器总线425与系统存储器420(例如，动态RAM)耦合。存储器集线器430还可以通过高级图形端口(AGP)总线433与视频控制器435耦合，该视频控制器可以与显示器437耦合。AGP总线433可以符合由加利福尼亚的圣克拉拉的英特尔公司于1998年5月4日公布的加速图形端口接口规范修订版2.0。
存储器集线器430还可以(通过集线器链路438)被耦合至与输入/输出(I/O)集线器440，而输入/输出(I/O)集线器440可与输入/输出(I/O)扩展总线442和如由PCI局部总线规范产品版本1995年6月的修订版2.1所定义的外围部件互连(PCI)总线444耦合。I/O扩展总线442可以与控制对一个或多个I/O设备的访问的I/O控制器446耦合。如图3所示，在一个实施例中这些设备可包括诸如软盘驱动器450之类的存储设备和诸如键盘452和鼠标454之类的输入设备。如图3所示，I/O集线器440还可与例如硬盘驱动器456和光盘(CD)驱动器458耦合。应理解系统中还可以包括其它存储介质。
PCI总线444还可以与各种部件(例如与网络端口(未示出)耦合的网络控制器460)耦合。其它设备可以与I/O扩展总线442和PCI总线444耦合，这些设备有例如与并行端口、串行端口耦合的输入/输出控制电路、非易失性存储器等。
虽然参照系统400的具体部件进行说明，但预期所述及所图示的实施例的许多修改和变更是可能的。特别是，虽然图3示出诸如个人计算机之类的系统的框图，但是应理解可以在诸如蜂窝电话、个人数字助理(PDA)等无线设备中实现本发明的实施例。
在某些实施例中，上述用于计算超越函数的无分支软件方法可以用系统400的处理器410的汇编语言编写。这种代码可以是将以特定源代码编写的较高级的程序编译成处理器410的机器代码的编译揣程序的一部分。
该编译器可包括根据常规技术对源代码进行语法分析并检测对超越函数的引用的操作。然后，编译器可以用合适的实现该超越函数的无分支方法的汇编语言指令序列来代替此高级函数调用的所有实例。特别是在某些实施例中，编译器可检测对正弦或余弦运算的调用，并用上述组合的sincos算法取代该调用。在其它实施例中，代码可以是诸如数学函数库等能用合乎需要的编程语言来调用的软件库的一部分。
虽然就有限数目的实施例说明了本发明，但本领域的技术人员将可理解源自本发明的许多修改和变更。旨在使所附权利要求书覆盖落在本发明的精神和范围内的所有这些修改和变更。
权利要求
1.一种方法，包括根据第一归约序列将函数的输入自变量x归约到范围被归约的值r；逼近具有主导部分f(A)+σr的对应的r的函数的多项式；以及使用所述多项式来得到所述函数的第一结果。
2.如权利要求1所述的方法，其特征在于，所述主导部分包括第一项f(A)和第二项σr，其中A等于x减r，而σ的绝对值为2的幂。
3.如权利要求1所述的方法，其特征在于，逼近所述多项式包括执行多个连续的加/减运算。
4.如权利要求1所述的方法，其特征在于，逼近所述多项式包括使用查找表来得到f(A)的断点。
5.如权利要求1所述的方法，其特征在于，还包括将所述输入自变量x限制于预定窗口内的值。
6.如权利要求1所述的方法，其特征在于，还包括将所述输入自变量x限制于2-252和90112之间的值。
7.如权利要求1所述的方法，其特征在于，得到所述函数的第一结果包括得到sin(x)。
8.如权利要求7所述的方法，其特征在于，还包括使用第二输入y来得到所述函数的第二结果，其中y比x大π/2。
9.如权利要求8所述的方法，其特征在于，得到所述函数的第二结果包括得到cos(x)。
10.如权利要求9所述的方法，其特征在于，还包括使用单指令多数据(SIMD)浮点运算来得到sin(x)利cos(x)。
11.如权利要求9所述的方法，其特征在于，还包括并行地得到所述第一结果和所述第二结果。
12.一种包括机器可访问存储介质的产品，所述机器可访问存储介质包含在被执行的情况下使系统能够执行以下方法的指令根据第一归约序列将函数的输入自变量x归约到范围被归约的值r；逼近具有主导部分f(A)+σr的对应的r的函数的多项式；以及使用所述多项式来得到所述函数的第一结果。
13.如权利要求12所述的产品，其特征在于，还包括在被执行的情况下使所述系统能够逼近所述多项式的指令，在所述多项式中，所述主导部分包括第一项f(A)和第二项σr，其中A等于x减r，而σ的绝对值为2的幂。
14.如权利要求12所述的产品，其特征在于，还包括在被执行的情况下使所述系统能通过使用查表来得到f(A)的断点以逼近所述多项式的指令。
15.如权利要求12所述的产品，其特征在于，还包括在被执行的情况下使系统能得到等于cos(x)的所述函数的第二结果的指令，其中所述第一结果等于sin(x)。
16.如权利要求15所述的产品，其特征在于，还包括在被执行的情况下使所述系统能使用单指令多数据(SIMD)浮点运算来得到sin(x)和cos(x)的指令。
17.如权利要求15所述的产品，其特征在于，还包括在被执行的情况下使所述系统能并行地得到所述第一结果和所述第二结果的指令。
18.一种系统，包括处理器；以及与所述处理器耦合的动态随机存取存储器，它包括在被执行的情况下使所述系统能根据第一归约序列将函数的输入自变量x归约到范围被归约的值r，逼近具有主导部分f(A)+σr的对应的r的函数的多项式，并使用所述多项式来得到所述函数的第一结果的指令。
19.如权利要求18所述的系统，其特征在于，所述动态随机存取存储器还包括在被执行的情况下使所述系统能得到等于cos(x)的所述函数的第二结果的指令，其中所述第一结果等于sin(x)。
20.如权利要求19所述的系统，其特征在于，所述动态随机存取存储器还包括在被执行的情况下使所述系统能使用单指令多数据(SIMD)浮点运算来得到sin(x)和cos(x)的指令。
21.如权利要求20所述的系统，其特征在于，所述动态随机存取存储器还包括在被执行的情况下使所述系统在函数调用请求sin(x)或cos(x)中的任意一个时能使用单指令多数据(SIMD)浮点运算来得到sin(x)和cos(x)的指令。
22.如权利要求20所述的系统，其特征在于，所述动态随机存取存储器还包括在被执行的情况下使所述系统能并行地得到所述第一结果和所述第二结果的指令。
全文摘要
在一个实施例中，本发明包括一种根据第一归约序列将函数的输入自变量x归约成范围被归约的值r，逼近具有主导部分f(A)+σr的对应的r的函数的多项式，并使用该多项式来得到该函数的结果的方法。
文档编号G06F7/548GK1918542SQ200580004840
公开日2007年2月21日 申请日期2005年3月4日 优先权日2004年3月11日
发明者J·哈里森, P·P·T·唐 申请人:英特尔公司