大數乘法,在演算法上,主要有幾種思路?
看了一些同取對數的和按高地位快速分治或模擬手算的累加型的方法,想知道在實際應用領域,如信息安全、大質數乘法、大數階乘等方面的主流實現演算法及相關數學計算原理為何?除了信息安全方面教材上出現比較多的RSA分治演算法,有別的么?
不需要代碼or偽代碼,給些提示link足以。
1、分治乘法(最簡單的是Karatsuba乘法,一般化以後有Toom-Cook乘法);
2、快速傅里葉變換FFT(為了避免精度問題,可以改用快速數論變換FNTT),時間複雜度O(NlgNlglgN)。具體可參照Sch?nhage;
3、中國剩餘定理(把每個數分解到一些互素的模上,然後每個同餘方程對應乘起來就行);
4、剛看到一個在漸進意義上FNTT還快的演算法Furer"s algorithm,不過好像不太實用。下面的reference[3]給出了模運算的版本。
//按照時域和頻域分別採樣的FFT,無需bit-reverse-copy
#include &
#include &
#include &
typedef long long sll;
#define MAXK 17 //至多容納2^MAXK個百進位數作NTT變換
#define MAXN 100010
#define BASE 100
using namespace std;
sll alpha[2][MAXK+1] = {0}; //alpha[0..1][k]是模M的2^k次單位根,其中alpha[0][]用於正變換,alpha[1][]用於逆變換
sll modulo = 1325137921;
char s[2][MAXN];
sll n[2][MAXN], res[MAXN];
int string_2_number(char x[], sll num[]){
x[0] = x[1] = "0"; //前面補0便於處理
int len = strlen(x+2), start = len 1;
int tot = (len &>&> 1) + 1; //採用百進位,4位一壓,數字串總長為tot
memset( num, 0, 4 * tot * sizeof(sll) ); //高位補0(最多需要補到原來的接近4倍)
for (int i = tot - 1; i&>=0; --i, start += 2)
num[i] = x[start] * 10 + x[start+1] - 528; // 實質為 num[i] = (x[start] - "0") * 10 + (x[start+1] - "0");
return tot;
}
inline sll power_mod (sll a,sll p,sll mod){
sll res = 1;
//a %= mod;
while (p){
if ( p1 )
res = res * a % mod;
a = a * a % mod;
p &>&>= 1;
}
return res;
}
void exGcd(sll a, sll b, sll gcd, sll x, sll y)
{
if(!b) {x = 1; y = 0; gcd = a;}
else {
exGcd(b, a % b, gcd, y, x);
y -= a/b * x;
}
}
inline void swap(sll x, sll y){sll tmp = x; x = y; y = tmp;}
//DIF-FFT: 輸入為自然順序,輸出為二進位倒序
void NTT_N2R(sll* a, int nBit){
int nLen = 1 &<&< nBit;
// butterfly operation
for (int layer = nBit; layer &>0; --layer)
{
int group = 1 &<&< layer, brother = group &>&> 1;
sll kernel = alpha[0][layer];
for (int k = 0; k &< nLen; k += group){ //遍歷每一組
sll w = 1;
for (int j = 0; j &< brother; ++j){ //對每一個元素,找到它的兄弟(和它進行蝴蝶操作的那個元素)
int cur = k + j, next = cur + brother;
sll u = a[cur], v = a[next];
a[cur] = u + v;
a[cur] = a[cur] &< modulo ? a[cur] : a[cur] - modulo;
a[next] = (u + modulo - v) * w % modulo; //中間運算結果最多達到2 * modulo^2,不會溢出(模更大時要當心!)
w = w * kernel % modulo;
}
}
}
}
//DIT-FFT: 輸入為二進位倒序,輸出為自然順序的FFT
void NTT_R2N(sll* a, int nBit){
int nLen = 1 &<&< nBit;
// butterfly operation
for (int layer = 1; layer &<= nBit; ++layer)
{
int group = 1 &<&< layer, brother = group &>&> 1;
sll kernel = alpha[1][layer], half = alpha[1][1];
for (int k = 0; k &< nLen; k += group){ //遍歷每一組
sll w = 1;
for (int j = 0; j &< brother; ++j){ //對每一個元素,找到它的兄弟(和它進行蝴蝶操作的那個元素)
int cur = k + j, next = cur + brother;
sll u = a[cur], t = w * a[next] % modulo;
a[cur] = u + t;
a[cur] = a[cur] &< modulo ? a[cur] : a[cur] - modulo;
a[next] = u - t;
a[next] = a[next] &< 0 ? a[next] + modulo : a[next];
w = w * kernel % modulo;
}
}
}
sll gcd = 0, invN = 0, invM = 0;
exGcd(nLen, modulo, gcd, invN, invM);
invN = (invN % modulo + modulo) % modulo;
for (int i=0; i&
if(high==-1)
printf("0
");
else{
printf("%d", a[high]);
for (int i=high-1; i&>=0; --i)
printf("%02d", a[i]);
printf("
");
}
}
int main(){
//alpha[0][]用於正變換,alpha[1][]用於逆變換
sll gcd, y;
exGcd(alpha[0][MAXK] = 1101238606, modulo, gcd, alpha[1][MAXK], y);
for (int i=16; i&>=0; --i){
alpha[0][i] = alpha[0][i+1] * alpha[0][i+1] % modulo;
alpha[1][i] = alpha[1][i+1] * alpha[1][i+1] % modulo;
}
while ( scanf("%s %s", s[0]+2, s[1]+2)==2 ){
int l1 = string_2_number(s[0], n[0]);
int l2 = string_2_number(s[1], n[1]);
int k = 0;
l1 = l1 &< l2 ? l2 : l1;
while ( (1&<&
大數乘法好像就是兩種演算法.冪運算倒是有幾個方法.
如果想系統了解推薦看看http://libtom.org裡面的東西。這個是開源的。裡面的一些演算法應該是最快演算法了。
大數乘法,GNU的GMP庫有非常好的實現,GMP的文檔中有對其中的演算法有描述,值得深入研究。
http://en.wikipedia.org/wiki/Multiplication_algorithm#Fast_multiplication_algorithms_for_large_inputs
http://en.wikipedia.org/wiki/Factorial#Computation請wiki:快速傅立葉變換與快速數論變換。
推薦閱讀: