当前位置：文档之家 › QuickMul Practical FFT-based Integer Multiplication

QuickMul Practical FFT-based Integer Multiplication

  • 格式：pdf
  • 大小：94.68 KB
  • 文档页数：4

QuickMul:PracticalFFT-basedIntegerMultiplicationCheeYapandChenLiDepartmentofComputerScienceCourantInstitute,NewYorkUniversityemail:{yap,chenli}@cs.nyu.edu

October6,2000

1IntroductionTheuseofarithmeticpackagesforarbitrarilylargeintegersisgrowinginmanyareasofapplication,beyondtheirtraditionalapplicationswhichismainlyincomputeralgebra.Forinstance,theBigIntegerclasslibrarythatisconsideredastandardpartofthepopularJavaprogramminglanguage.Onerecentareaofapplicationisinrobustgeometricalgorithms(e.g.,[4,2]).Thecriticalalgorithminalltheseapplicationsistheintegermultiplicationalgorithm.ThefastestknownalgorithmhereisduetoSch¨onhageandStrassen(1971)[5],achievingthetimeboundT(N)=O(NlogNloglogN)(1)

formultiplyingtwoN-bitintegers.Inthisnote,weareinterestedinexploringintegermultiplicational-gorithmswhich,liketheSch¨onhage-Strassenalgorithm,arebasedontheFastFourierTransform(FFT).Thesealgorithmsmaynotachievetherecordbound(1),buthasorderT(N)=O(NlogO(1)N)withsmallimplicitconstants.ThehopeisthatthesmallimplicitconstantsmaymakesuchalgorithmsmoreefficientforpracticalvaluesofN.TheoriginalSch¨onhage-Strassenalgorithmisrelativelycomplex;asimplifiedalgorithmisgivenin[8,chap.1]withtimeboundO(Nlog1+ǫN)foranyǫ.Suchsimplifiedalgorithmsarequiteeasytoimplement(asthispaperwillshow).Becausesuchsimplificationsarenotwell-known,manyimportantbignumberpackagescontinuetoavoidFFT-multiplicationsalgorithms.TwonotableexceptionsamongthefreelyavailablepackagesareDavidBailey’sMPFUNpackage(writteninFortran)[1]andBrunoHaible’sCLN(writteninC++)[3].Forinstance,onlyrecently(August2000)didthewidelyusedGnugmppackageimplementFFT-multiplication.

2TheOne-PrimeFFTMultiplicationThepresentpaperimplementsanFFT-multiplicationalgorithmdescribedin[7].ThebasicideaistoperformtheFFTintheringZM={0,1,...,M−1}ofnumbersmoduloM.HereMisspeciallychosenprimenumber.Itwassuggestedin[7]thatsuchanalgorithmshouldhave“small”constants,butinputnumbersmustbelessthansomefinitebutverylargelimit.Thisfinitelimitationplussmallconstantsiswhatwemean1by“practical”inthetitle.Typically,smallconstantsimplythatthealgorithmissimple.Inparticular,someoptimizationstepsintheSch¨onhage-Strassenalgorithmwillbediscarded.Wespecifychoosetwo“gamerules”helptofurtherensuresmallconstants:(1)assumea32-bitmachineand(2)insistthatthemultiplicationalgorithmisnon-recursive.Rule(1)meansthatwewantMatmost32-bitssothatarithmeticinZMcanbeperformedinO(1)time.Rule(2)doesnotimplythatwedonotuserecursionatall–indeed,theFFTalgorithmwhichwewilluseisinherentlyrecursive.However,therecursioninFFTalgorithmsisrelativelysimpleanddonotincurlargeoverhead.So,whatdoesrule(2)reallyexclude?AsimplifiedFFT-basedmultiplicationalgorithmhasthefollowingbasicsteps.AssumewearegiventwoN-bitintegersUandV,andwewanttocomputetheirproduceW=UV.

I.PreparationBreakupUandVintotwoK-vectorsVin(ZM)K(forsuitablenaturalnumbersKandM).

II.DFTComputationComputetheirdiscreteFouriertransforms,DFT(V)∈(ZM)K.III.Component-wiseProductComputethecomponent-wiseproductofDFT(V).Thisresultsinavectorin(ZM)KwhichwemaydenoteasDFT(

W)toobtainWistheconvolutionofV.

V.Re-assemblyItisrelativelystraightforwardtore-assembleW=U·Vfrom

U=(0,0,...,0,UK−1,UK−2,...,U0).Similarly,letLetusestimatethelargestvalueofNwhichcanbeachievedunderrules(1)and(2).Supposewecomputethe2K-vectorWistheconvolutionofV,itiseasytoseethateachWiisthesumof≤Knumbers,eachwith≤Lbits.ThusWihas≤2L+lg(K)bits.[Note:lgislogarithmtobase2.]Soweneed

2L+lg(K)≤lgM=30.9....SinceLandlg(K)areintegers,wehavehave2L+lg(K)≤30.ClearlyN=KLismaximizedbymakingKaslargeaspossible.LetNmaxbethismaximumvalue.SettingL=1,wehavelg(K)≤(lgM)−2orK≤M/4=503,316,480.SinceKisapowerof2,weobtainlg(K)=⌊lg(503,316,480)⌋=28orNmax=228.Sothismethodworksforintegersupto256megabitsor32MBlong(counting“1024”ratherthan“1000”asonekilo).Ingeneral,giventwonumberseachwithN≤Nmaxbits,wechoosepositiveK,LsoastominimizeK(andhencemaximizeL)andsubjectto

KL≥N,2L+lg(K)≤lgM.ThefollowingtablegivestheupperboundonNforanyL=1,2,...,8.Llg(K)814716618520422324226128

Forinstance,forN=1,000,000,wechooseL=6andK=218=262,144.Howfastcanwemultiplytwosuchnumbers?Well,FFTon2K-vectorstakes1.5(2K)lg(2K)(properlyimplemented)arithmeticoperations.Thisgives3·218·19orabout15millionfloatingpointoperations(flops).Countingafactoroffperflops,wehaveabout15fsecondsonamegaflopmachine.Thisfactorfcanbeestimatedforvariousplatforms.Forinstance,iff=20,thisis300seconds.

3ImplementationofQuickMulWehaveimplementedouralgorithminC++,inarelativelystraightforwardimplementation.Thiscodeandexperimentalresultscanbeobtainedfromourwebsitehttp://cs.nyu.edu/exact/.ThefollowingisatablecomparingouralgorithmQMUL(for“quickMul”)tothelatestGnu’sgmp(version3.1),BrunoHaible’sCLN(asincorporatedintoLiDIA).Forhistorical2reasons,wealsoincludeacomparisonwiththeoriginalGnu’sbigInteger.Theplatformforthesetestsisa“SunUltraSPARC-IIi”440MHzCPUmachinewith512MBmainmemory.

合集下载

快速傅里叶变换FFT的C语言实现及应用

快速傅里叶变换FFT的C语言实现及应用

快速傅里叶变换FFT的C语言实现及应用快速傅里叶变换（Fast Fourier Transform, FFT）是一种快速计算离散傅里叶变换（Discrete Fourier Transform, DFT）的算法。

它能够在较短时间内计算出巨大数据集的傅里叶变换，广泛应用于信号处理、图像处理、通信等领域。

C语言是一种广泛应用于嵌入式系统和科学计算的编程语言，拥有高效、灵活和可移植等特点。

下面将介绍FFT的C语言实现及应用。

首先，我们需要了解离散傅里叶变换的定义。

离散傅里叶变换将一组离散的时域信号转换成一组对应的频域信号，可以表示为以下公式：X(k) = Σ[ x(n) * W^(kn) ]其中，X(k)是频域信号，x(n)是时域信号，W是单位复数旋转因子，其计算公式为W=e^(-j*2π/N)，其中j是虚数单位，N是信号的长度。

实现FFT算法的关键在于计算旋转因子的值，一种常用的计算方法是采用蝶形算法（butterfly algorithm）。

蝶形算法将DFT分解成多个子问题的求解，通过递归调用可以快速计算出结果。

以下是一种基于蝶形算法的FFT实现的示例代码：```c#include <stdio.h>#include <math.h>typedef structfloat real;float imag;if (N <= 1)return;}for (int i = 0; i < N/2; i++)even[i] = signal[2*i];odd[i] = signal[2*i + 1];}fft(even, N/2);fft(odd, N/2);for (int k = 0; k < N/2; k++)signal[k].real = even[k].real + temp.real;signal[k].imag = even[k].imag + temp.imag;signal[k + N/2].real = even[k].real - temp.real; signal[k + N/2].imag = even[k].imag - temp.imag; }int maiint N = sizeof(signal) / sizeof(signal[0]);fft(signal, N);printf("频域信号：\n");for (int i = 0; i < N; i++)printf("%f + %fi\n", signal[i].real, signal[i].imag);}return 0;```以上代码实现了一个简单的4点FFT算法，输入时域信号为{1,0,1,0}，输出为对应的频域信号。

《快速傅里叶变换》课件

《快速傅里叶变换》课件

FFT算法的出现极大地推动了数字信号 处理技术的发展和应用。
FFT的历史背景
01
1960年代，Cooley和Tukey提 出了基于“分治”思想的FFT 算法，为快速傅里叶变换的实 用化奠定了基础。
02
随后，出现了多种FFT算法的 变种和优化，如Radix-2、 Radix-4等。
03
随着计算机技术的发展，FFT 算法在硬件实现上也得到了广 泛应用，如FPGA、GPU等。
《快速傅里叶变换》ppt课件
contents
目录
• FFT简介 • FFT基本原理 • FFT实现 • FFT的应用 • FFT的优化与改进 • FFT的挑战与未来发展
01 FFT简介
FFT的定义
快速傅里叶变换（FFT）：一种高效计算离散傅里叶变换（DFT）及其逆变换的 算法。它将复杂度为$O(N^2)$的DFT计算降低到$O(Nlog N)$，大大提高了计 算效率。
详细描述
混合基数FFT算法结合了基数-2和基数-4算法的特点，利用两者在计算过程中的 互补性，减少了计算量，提高了计算效率。同时，该算法在处理大规模数据时 ，能够保持较高的精度。
分段FFT算法
总结词
分段FFT算法将输入数据分成若干段，对每一段进行快速傅里叶变换，以降低计算复杂度和提高计算效率。
详细描述
02 FFT基本原理
离散傅里叶变换（DFT）
定义
应用
DFT是时间域信号到频域的变换，通 过计算信号中各个频率成分的幅度和 相位，可以分析信号的频谱特性。
DFT在信号处理、图像处理、频谱分 析等领域有广泛应用。
计算量
DFT的计算量随着信号长度N的增加 而呈平方关系增长，因此对于长信号 ，计算量巨大。
C语言实现FFT

C语言实现FFT

C语言实现FFTFFT(快速傅里叶变换)是一种常用的算法，用于计算离散傅里叶变换(DFT)的快速算法。

它在信号处理、图像处理、通信等方面有广泛的应用。

C语言提供了一组标准库函数来支持FFT算法的实现。

下面以C语言为例，展示如何实现FFT算法。

1.理解DFT首先，我们需要理解离散傅里叶变换(DFT)的概念。

DFT将时域离散信号转换为频域离散信号，它的计算公式如下：其中N是信号的长度，k表示频域的频率，n表示时域的时间。

离散信号经过DFT变换后，可以得到相应频率的幅度和相位信息。

2. Cooley-Tukey算法FFT算法采用了Cooley-Tukey算法的思想，它的主要思路是将DFT问题递归地分解为更小的DFT问题。

这样可以减少计算量，提高计算效率。

Cooley-Tukey算法具体如下：-如果信号的长度N是2的整数次幂，将原始信号分为偶数索引和奇数索引的两部分。

-对分离出的偶数索引部分和奇数索引部分分别进行DFT变换。

-最后将两个结果进行合并。

下面是一个简单的例子，展示了如何使用C语言实现FFT算法：```c#include <stdio.h>#include <math.h>if (N <= 1) return;for (int i = 0; i < N/2; i++) even[i] = x[2*i];odd[i] = x[2*i + 1];}fft(even, N/2);fft(odd, N/2);for (int k = 0; k < N/2; k++) x[k] = even[k] + t;x[k + N/2] = even[k] - t;}int maiint N = 8; // 信号长度//输入信号x[0]=1;x[1]=1;x[2]=1;x[3]=1;x[4]=0;x[5]=0;x[6]=0;x[7]=0;fft(x, N);//打印频域的幅度信息for (int k = 0; k < N; k++)printf("Amplitude at frequency %d: %f\n", k, cabs(x[k]));}return 0;```在上面的示例中，我们首先定义了一个fft函数，用于实现FFT算法。

数字信号处理 第4章 FFT基本思想和2种基本的FFT

数字信号处理 第4章 FFT基本思想和2种基本的FFT


= −W
W的对称性
W的可约性
2 rk WN rk = WN / 2
长序列变成短序列 若N → 2个N / 2
2 则N 2次复述乘法 →（N / 2）= N 2 / 2次复数乘法 2
从信号的特殊性上考虑
– 如奇、偶、虚、实性
W 0 X (0) X (1) W 0 = X (2) W 0 0 X (3) W
对 N = 2M ， 共可分 M 次，即 m = 0,1,L , M − 1，
8点FFT时间抽取算法信号流图
每一级有 N/2 个如下的“蝶形”单元：
xm ( p )
xm +1 ( p )
W
r N
xm (q)
−1
xm +1 (q )
算法讨论（ “级”的概念、碟形单元、 “组” 的概念、旋转因子的分布、码位倒置）
r =2l ,r =2l +1
A(k ), B(k )
C(k) = D(k) =
N / 4−1 l =0
∑x(4l)W
l =0
lk N/4
, k = 0,1,..., N / 4 −1
N / 4−1
lk x(4l + 2)WN / 4 , k = 0,1,..., N / 4 −1 ∑
k A(k) = C(k) +WN / 2 D(k), k = 0,1,..., N / 4 −1 k A(k + N / 4) = C(k) −WN / 2 D(k), k = 0,1,..., N / 4 −1
x(6)
n N
N n = 0,1,L , 2
由此得到基本 运算单元
g (0) g (1) g (2) g (3)