关键字 |
| 快速傅里叶变换;大整数乘法;GPGPU;CUDA编程;cuFFT |
我的介绍。 |
| 大整数相乘在计算科学中有很多应用。许多密码算法都需要对非常大的整数子集进行操作。一个常见的应用是公钥密码学，其算法通常使用具有数百位数字[1]的整数的算术。任意精确算术也被用于计算基本的数学常数，如π到数百万或更多的数字，并分析数字串的性质，或更一般地研究函数的精确行为，如黎曼ζ函数，其中某些问题难以通过分析方法探索。另一个例子是用极高的放大倍率渲染分形图像，例如在Mandelbrot集合[2]中发现的那些图像。 |
| 三种最流行的大整数乘法算法是Karatsuba-Ofman [3]， Toom-Cook[4]和FFT乘法[5]算法。经典的乘法运算有O(n)2)复杂度，其中n是位数。通过使用FFT的多项式乘法，它的时间复杂度为O(nlogn)，我们可以显著减少乘法所需的时间。对于这种基于大数据量的应用，基于图形处理单元(GPU)的算法是一种经济有效的解决方案。GPU可以在单指令多数据(SIMD)模式下并行处理大量数据。2006年11月，NVIDIA[6]发布了专门用于计算密集型高度并行计算的计算统一设备架构(CUDA)。NVIDIA CUDA快速傅里叶变换库(cuFFT)为计算fft提供了一个简单的接口，最快可达10倍。通过在NVIDIA GPU内部使用数百个处理器核心，cuFFT提供了GPU的floatingâ ' Â点性能，而无需开发自己的自定义GPU FFT实现。cuFFT使用基于著名的Cooley-Tukey和Bluestein算法的算法，因此您可以自信地获得比以往更快的准确结果[8]。 |
2相关工作 |
| 首先让我们讨论一些库和框架，它们可以执行任意精确算术，特别是大整数乘法。[9]是微软在。net 4.0中引入的，是用于大整数的。net类型。BigInteger类型是一种不可变类型，它表示任意大的整数，其值理论上没有上界或下界。BigInteger类型的成员与其他整型类型的成员非常相似。因为BigInteger类型是不可变的(参见可变性和BigInteger结构)，并且因为它没有上界或下界，所以任何导致BigInteger值增长过大的操作都可以抛出OutOfMemoryException。其中一个是IntX[10]。IntX是一个任意精度整数库，使用快速哈特利变换可以快速乘法大整数。GNU MP库[11]-世界上最快的著名的Bignum库之一。GNU MP(或GMP)是用普通C语言和汇编语言编写的，所以它可以编译成优化的本机代码，它还使用了快速计算算法——这就是为什么它非常快。GMP是用于任意精确算术的免费库，可以对有符号整数、有理数和浮点数进行操作。 There is no practical limit to the precision except the ones implied by the available memory in the machine GMP runs on. Among other multiplication algorithms, GMP also uses FFT multiplication (depending on operand size). In [12] implemented FFT multiplication algorithm and done experiments by comparing FFT multiplications with normal multiplications at various bases. Also there is an existing large number multiplication implementation on CUDA by K. Zhao [13], where implemented multipleprecision Karatsuba and Montgomery multiplication algorithms. |
3离散傅里叶变换 |
| n点离散时间信号f(x)的一维离散傅里叶变换(DFT)由以下公式给出: |
 (1） |
| 对于u = 0,1,2，…， n−1。 |
| 类似地，对于给定的F(u)，我们可以通过DFT逆得到原始离散函数F(x): |
 （2） |
| 对于x = 0,1,2，…， n−1。 |
| 离散傅里叶变换经常对每个数据样本进行评估，可以视为从信号中提取特定的频率成分。 |
四，快速傅里叶变换 |
| 用c的元素计算F的简单方法n需要O (n2)操作。然而，有一些算法可以在O(nlogn)步中计算F。我们计算z的傅里叶变换= (z1,……, z2 n)．让F2 n表示相对于C的傅里叶变换2 n令Fn表示相对于C的傅里叶变换n．定义 |
 （3） |
| 让[x]k表示kth向量x的坐标 |
 （4） |
| 我们定义 |
 （5） |
| 这里ZE和Zo分别是由z的偶分量和奇分量构成的向量。用这种表示法，式4可以更简洁地写成如下: |
 (6） |
| 方程6适用于所有k = 0，…，2n - 1，但我们只需要计算Ek和Okk = 0,…，n－1because |
 (7） |
| 假设需要T(n)次运算才能计算C中向量的傅里叶变换n．上面的技巧表明我们可以计算C中向量的傅里叶变换2 n使用2 t (n) + 8 n。下面是计算的细目。 |
| •我们可以计算{Ok}和{Ek}对于k=0，…，n－1in 2T(n) steps. |
| •我们可以计算ω2 nkk = 0,…，2n－1in 2n steps. |
| •计算公式6的每个实例需要3个以上步骤。所以总共是6n个额外的步骤。 |
| 显然T(2°)≤8。一个简单的归纳论证表明 |
 （8） |
| 这表明，当n =2时k向量C的傅里叶变换n可以在O(n log (n))步内计算。值得一提的是，这里的“step”指的是比简单浮点运算更复杂的操作。例如，一个典型的步骤是将大小为log (n)的复数与单位的n次方根相乘。计算傅里叶变换所需的浮点运算数量的实际分析将取决于这些单独步骤的效率。 |
V. gpu架构和cuda FFT (cufft)库 |
| gpu是大型多线程多核芯片。NVIDIA Tesla产品拥有多达128个标量处理器，超过12,000个并发线程，超过470个GFLOPS的持续性能。NVidia显卡架构由多个所谓的流多处理器(SM)组成。每一个都包括8个着色处理器(SP)核心，一个由所有SP共享的本地内存，16384个寄存器，以及用于超越函数硬件加速的快速ALU单元。全局内存由所有SMs共享，最大容量为4gb，内存带宽为144gb /s。FERMI架构引入了配备32个sp和32768寄存器的新型SMs，用于快速双精度浮点性能的改进ALU单元和L1缓存[14]。 |
| CUDA是一种可扩展的并行编程模型和并行计算的软件环境。它非常便于程序员使用，为C语言引入少量扩展，以提供并行执行。另一个重要的特性是数据结构的灵活性，对GPU不同物理内存级别的显式访问，以及一个适合程序员的良好框架，包括编译器、CUDA软件开发工具包(CUDA SDK)、调试器、分析器以及cuFFT和cuBLAS科学库。GPU以SIMT(单指令多线程)方式执行指令。一个典型CUDA程序的执行如图1所示。 |
| cuFFT的关键功能是[8] |
| •复杂和真实数据类型的1D, 2D, 3D转换 |
| •一维变换大小可达1.28亿个元素 |
| •灵活的数据布局，允许单个元素和数组尺寸之间的任意步幅 |
| •基于Cooley-Tukey和Bluestein的FFT算法 |
| 熟悉类似FFTW高级接口的API |
| •流异步执行 |
| •单和双精度转换 |
| •批处理执行多个转换 |
| •就地变换和就地变换 |
| •灵活的输入和输出数据布局，类似tp FFTW“高级界面” |
| •线程安全&可从多个主机线程调用 |
六、cuda FFT大整数乘法算法 |
| 我们现在给出了用cuFFT乘大整数的算法。基本思想是使用快速多项式乘法来执行快速整数乘法。图2所示。cuFFT乘法算法流程图。让我们一步一步地讨论这张图的每一点。 |
| (颜色编码:白色主机代码，黄色内存复制操作，绿色设备代码) |
| 1.输入值(a和b)可以以不同的方式给予程序输入。我们可以从文件或数据库中读取它们，或者从程序接口中设置它们。在我们的例子中，用户可以从接口插入输入值，但出于测试和实验的目的，我们可以选择生成随机的大整数。对于数字生成，我们使用由用户设置的位数。 |
| 2.对于FFT乘法，首先我们必须用多项式表示一个整数。多项式的默认符号是它的系数形式。以系数形式表示的多项式p可用其系数向量A = {A来描述0,一个1……,n - 1}如下: |
 （9) |
| 我们称x为多项式的基，p(x)是多项式的求值，由它的系数向量a定义，以x为基。两个多项式相乘会得到第三个多项式，这个过程称为向量卷积。与整数相乘一样，向量卷积需要O(n2)时间。但是卷积在离散傅里叶变换下变成了乘法(时域卷积可以通过频域乘法实现): |
| F(c) = F(a) F(b) |
| c = F－1F (F (a) (b)) |
| 当我们用多项式表示一个整数时，我们可以选择使用什么底数。任何正整数都可以用作底数，但为了简单起见，我们将底数限制在10。考虑整数12345，它的多项式形式使用base = 10是a ={5,4,3,2,1}。由于FFT处理复数，我们必须得到复数的向量。为此，通过将虚数部分设为零，我们得到一个={[5,0]，[4,0]，[3,0]，[2,0]，[1,0]}向量。 |
| 3.CUDA内核blockSize和gridSize的确定。我们通过实验找到了CUDA内核的最佳块大小，并根据块大小计算网格大小。 |
| 4.这些传输是GPU计算中数据移动最慢的部分。实际的传输速度(带宽)取决于您所使用的硬件类型，但不管这一点，它仍然是最慢的。设备内存和设备处理器之间的峰值理论带宽明显高于主机内存和设备内存之间的峰值理论带宽。因此，为了在我们的应用程序中获得最大的收益，我们确实需要最小化这些主机设备数据传输。如果在整个应用程序中发生了多次传输，请尝试减少这个数量并观察结果。在我们的例子中，我们有两个“从主机复制数据到设备”操作: |
| cudaMemcpy(dev_a, a, numBytes, cudaMemcpyHostToDevice); |
| (dev_b, b, numBytes, cudaMemcpyHostToDevice); |
| 5.这里我们对数组a和数组b执行两个正向fft。结果arrayA和arrayB的乘法多项式，arrayC的阶比最高阶的arrayA和arrayB的阶大两倍。因此，度数组c将被计算为: |
| signalSize = 2(a;长度>长度?a.长度: |
| 在进行前向傅里叶变换之前，我们将arrayA和arrayB的系数填充为零，直到signalSize。我们使用signalSize和cufftype.c2c执行cuFFT计划。 |
| 6.此时，我们在设备上有两个传输的数组，我们必须执行成对乘法。因为我们要处理复杂的数据，我们有: |
| 复杂的c; |
| C.x = a.x * b.x - a.y * b.y; |
| C.y = a.x * b.x + a.y * b.y; |
| 7.这里我们对数组c执行逆FFT |
| 8.cuFFT使用的变换公式是非正交的根号因子(signalSize)。通过signalSize和1/signalSize归一化是使用FFTs计算傅立叶级数系数时所需要的，因此我们将点7的结果除到设备上的signalSize。 |
| 9.这一点与第4点类似，只是这里我们将内存从设备复制回主机 |
| 10.携带传播是得到结果多项式的最后一步。我们将给出这个操作的伪代码: |
| 伪代码1。大整数乘法的进位传播 |
| 输入:整数数组sourceArray |
| 输出:整数数组resultArray |
| 1.CarryPropagation () |
| 2.进位:= 0 |
| 3.for (i from 0 to sourceArray.length) |
| 4.sum:= carry + sourcarray [i] |
| 5.Mod:= sum % 10 |
| 6.resultArray.Add (mod) |
| 7.进位:= sum / 10 |
| 8.结束了 |
| 9.While(进位> 0) |
| 10.如果携带。长度> 1) |
| 11.指数:= carry。长度- 1 |
| 12.其他的 |
| 13.指数:= 0 |
| 14.结束 |
| 15.resultArray.Add(携带(指数)) |
| 16.Carry:= Carry / 10 |
| 17.结束时 |
| 11.在最后一步中，我们将多项式表示为整数，并将结果发送到输出。 |
7实验与结果 |
| 在采用Intel(R) Core(TM) i3-2100 3.10GHz处理器和4gb主存的台式电脑上，我们实现了CUDA 5.5的cuFFT乘法算法，并使用GeForce GT 630显卡进行了一系列基准测试。该显卡具有2.1的计算能力，由96个CUDA组成，用于整数和单精度浮点运算。我们首先改变块大小，并根据块大小计算网格大小，以发现它们对性能的影响。测试中随机生成大整数。我们生成具有不同位数的整数。表1展示了块大小对性能的影响。我们以毫秒为单位计算一次乘法运算所需的时间。在进行100次测试乘法运算后，我们对每个数据大小取平均值。计算时不考虑大整数生成的时间。 |
| 图3所示。实验结果示意图见表1。正如我们所看到的，当内核的块大小为256时，我们有最好的乘法性能。我们用下面的公式计算网格大小 |
| gridSize = vectorSize/ blockSize + (vectorSize % blockSize != 0 ?1: 0) |
| 其中vectorSize是乘法多项式的长度，“/”和“%”是div和mod操作。图3横轴表示整数相乘的位数，纵轴为乘法时间，单位为毫秒。块大小决定了在单个块中可以并行执行的线程数。块大小的最大值基于GPU。对于具有计算能力1的gpu，每个块最多可以有512个线程。对于具有计算能力的gpu，每块1024个线程。X和3.x。 |
| Ando Emerencia在他的论文中介绍了FFT乘法算法及其与正常乘法的比较，在不同的基数[12]。图4给出了以10为基数的[12]的结果。我们将讨论以10为基底，因为我们的实验也是以10为基底进行的。所采取的个别测量值也列在下表中。我们在这里看到，对于超过200的输入大小，FFT方法变得越来越快于普通的乘法。对于输入大小为10的5， FFT乘法算法需要258ms，而普通乘法需要28秒以上，因此所需时间相差超过100倍。 |
| 我们在这里看到，对于超过200的输入大小，FFT方法变得越来越快于普通的乘法。对于输入大小为10的5， FFT乘法算法需要258ms，而普通乘法需要28秒以上，因此所需时间相差超过100倍。如果我们比较图4和图3，我们可以看到cuFFT乘法快了很多倍，输入大小为10时大约快了5倍5． |
| 在大多数计算机软件中，任意精度算术是通过调用一个外部库来实现的，该外部库提供数据类型和子例程，以存储具有所要求精度的数字并执行计算。不同的库有不同的表示任意精度数的方法。有许多任意精度的算术库可以执行大整数乘法。我们选择其中三个来进行实验: |
| •微软。net system . numbers . biginteger [9] |
| •IntX [10] |
| GNU MP(或GMP) [11] |
| 在表2中给出了cuFFT乘法与。net BigInteger、IntX和GMP库乘法的比较。对不同大小的输入整数(位数)进行实验。所提供的实验结果是一次乘法运算所需的时间，单位为毫秒。我们为每个数据大小和库做了100次乘法，并取平均值。图5为表2的图表表示。 |
| 从图中我们可以看到，对于5000位的整数，cuFFT比。net BigInteger和IntX快，但GMP是最快的。cuFFT从50000位整数开始比GMP快。数据大小为105．cuFFT比GMP快2倍左右。 |
8结论 |
| 正如我们所看到的，我们使用基于并行GPU的乘法算法获得性能，并使用CUDA的cuFFT库。实验表明，当我们处理大约2时，cuFFT乘法变得比所有其他测试方法都快15位数的整数。 |
表格一览 |
 |
|
| 表1 |
表2 |
|
|
数字一览 |
|
|
|
参考文献 |
- R. L. Rivest, A. Shamir和L. Adleman，“一种获取数字签名和公钥密码系统的方法”，ACM通讯，Vol. 21, pp. 120-126, 1978。
- R. Brooks和P. Matelski，“PSL(2,C)的2-产生子群的动力学”，黎曼曲面及其相关主题，第97卷。pp.65 - 71, 1997。
- A. Karatsuba和Yu。Ofman，“自动计算机多数字乘法”，DokladyAkad。Nauk SSSR卷145，页。293 - 294年,1962年。
- “图姆-库克乘法”，http://en.wikipedia.org/wiki/Toom-Cook_multiplication
- A. Schönhage和V. Strassen，“schnellemultiplikationgo ß erzahlen”，计算，Vol. 7, pp. 281-292, 1971。
- NVIDIA CUDA C编程指南，英伟达公司，[在线]。可用:www.nvidia.com, 2012。
- J. Sanders和E. Kandrot，CUDA实例， Addison-Wesley, 2010。
- “cuFFT”,https://developer.nvidia.com/cufft
- ”。Net BigInteger”，http://msdn.microsoft.com/en-us/library/system.numerics.biginteger
- “IntX”,https://intx.codeplex.com/
- “GMP”,https://gmplib.org/
- A. Emerencia，“使用傅里叶变换乘法巨大整数”，2007，[在线]。可用:http://www.cs.rug.nl/ ~安藤/ pdf / Ando_Emerencia_multiplying_huge_integers_using_fourier_transforms_paper.pdf
- 赵凯，“GPU上的多精度模块化乘法的实现”，技术研究，2009。
- D. Luebke，“并行计算的民主化”，NVIDIA公司，2007年。
|
菜单
(1）
（2）
（3）
（4）
（5）
(6）
(7）
（8）
（9)
