总共四个文件。
/*
文件名
Bignum.h
*/
#include <iostream>
#include <cstdio>
#include <string>
#include <cstring>
#include <vector>
#include "FFT.h"
#include "Vectorcal.h"
using namespace std;
struct Bignum {
int s[5010];
bool flag;
int len;
Bignum()
{
memset(s, 0, sizeof(s));
flag = 0;
len = 0;
}
void init(string num)
{
memset(s, 0, sizeof(s));
flag = len = 0;
int i = 0;
if(num[0]=='-')
{
flag = 1;
i++;
}
for(int p = 0;num[i]!='\0';p++,i++)
{
s[p] = num[i]-'0';
len++;
}
}
void print()
{
if(flag==1)
printf("-");
for(int i=0;i<len;i++)
printf("%d", s[i]);
printf("\n");
}
void rev()
{
int i = 0, j = len-1;
while(i<j)
{
swap(s[i], s[j]);
i++; j--;
}
}
bool zero()
{
if(len==1&&s[0]==0)
return 1;
return 0;
}
//赋值符号
Bignum operator = (const Bignum b)
{
flag = b.flag;
len = b.len;
for(int i=0;i<len;i++)
{
s[i] = b.s[i];
}
return *this;
}
//重载比较符号
bool operator == (const Bignum b)
{
if(flag!=b.flag||len!=b.len)
return false;
for(int i=0;i<len;i++)
{
if(s[i]!=b.s[i])
return false;
}
return true;
}
bool operator < (const Bignum b)
{
bool res = 0;
if(len<b.len)
res = 1;
else if(len>b.len)
res = 0;
else
{
for(int i=0;i<len;i++)
{
if(s[i]<b.s[i])
{
res = 1;
break;
}
else if(s[i]>b.s[i])
{
res = 0;
break;
}
}
}
if(flag==0&&b.flag==0)
return res;
else if(flag==1&&b.flag==1)
return !res;
else if(flag==1&&b.flag==0)
return true;
else if(flag==0&&b.flag==1)
return false;
return false;
}
bool operator > (const Bignum b)
{
bool res = 0;
if(len>b.len)
res = 1;
else if(len<b.len)
res = 0;
else
{
for(int i=0;i<len;i++)
{
if(s[i]>b.s[i])
{
res = 1;
break;
}
else if(s[i]<b.s[i])
{
res = 0;
break;
}
}
}
if(flag==0&&b.flag==0)
return res;
else if(flag==1&&b.flag==1)
return !res;
else if(flag==1&&b.flag==0)
return false;
else if(flag==0&&b.flag==1)
return true;
return false;
}
int CMP(Bignum a, Bignum b)
{
if(a<b)
return -1;
else if(a==b)
return 0;
else if(a>b)
return 1;
}
//重载四则运算
Bignum operator + (const Bignum num)
{
Bignum c;
Bignum a = *this, b = num;
if(a.flag==b.flag)
{
c.flag = a.flag;
VectorAdd(a.s, a.len, b.s, b.len, c.s, c.len);
}
else
{
if(VectorCmp(a.s, a.len, b.s, b.len)) //a的绝对值更大
{
c.flag = a.flag;
VectorSub(a.s, a.len, b.s, b.len, c.s, c.len);
}
else
{
c.flag = b.flag;
VectorSub(b.s, b.len, a.s, a.len, c.s, c.len);
}
}
return c;
}
Bignum operator - (const Bignum num)
{
Bignum c;
Bignum a = *this, b = num;
if(a.flag!=b.flag)
{
c.flag = a.flag;
VectorAdd(a.s, a.len, b.s, b.len, c.s, c.len);
}
else
{
if(VectorCmp(a.s, a.len, b.s, b.len))
{
VectorSub(a.s, a.len, b.s, b.len, c.s, c.len);
c.flag = a.flag;
}
else
{
VectorSub(b.s, b.len, a.s, a.len, c.s, c.len);
c.flag = !a.flag;
}
}
return c;
}
Bignum operator * (const Bignum num)
{
Bignum c;
Bignum a = *this, b = num;
if(a.flag==b.flag)
c.flag = 0;
else
c.flag = 1;
/* a.rev();
b.rev();
for(int i=0;i<a.len;i++)
{
int x = 0;
for(int j=0;j<b.len;j++)
{
c.s[i+j] = a.s[i]*b.s[j]+x+c.s[i+j];
x = c.s[i+j]/10;
c.s[i+j] %= 10;
}
c.s[i+b.len] = x;
}
c.len = a.len+b.len;
while(c.s[c.len-1]==0&&c.len>1)
c.len--;
c.rev();
return c;*/
VectorMul(a.s, a.len, b.s, b.len, c.s, c.len);
return c;
}
Bignum operator / (const Bignum num)
{
Bignum c;
Bignum a = *this, b = num;
if(a.zero())
{
c.init("0");
return c;
}
if(a.flag==b.flag)
c.flag = 0;
else
c.flag = 1;
a.flag = b.flag = 0;
c.len = a.len-b.len+1;
if(c.len<=0)
c.len = 1;
/*
for(int i=c.len-1;i>=0;i--)
{
Bignum temp;
temp = b;
temp.len += i;
while(a>temp||a==temp)
{
a = a-temp;
c.s[i]++;
}
}
while(c.s[c.len-1]==0&&c.len>1)
c.len--;
c.rev();*/
if(a<b)
{
c.len = 1;
c.s[0] = 0;
return c;
}
else if(a==b)
{
c.len = 1;
c.s[0] = 1;
return c;
}
VectorDiv(a.s, a.len, b.s, b.len, c.s, c.len);
return c;
}
Bignum operator % (const Bignum num)
{
Bignum c;
Bignum a = *this, b = num;
Bignum temp = (a/b);
temp = temp*b;
c = a-temp;
return c;
}
};
Bignum qpow(Bignum a, Bignum b, Bignum mod)
{
Bignum div, res, t;
div.init("2");
res.init("1");
t = a;
while(!b.zero())
{
if(b.s[b.len-1]%2==1)
{
res = res*t;
res = res%mod;
}
t = t*t;
t = t%mod;
b = b/div;
}
return res;
}
Bignum GCD(Bignum a, Bignum b)
{
Bignum z;
z.init("0");
if(a<b)
{
Bignum t = a; a = b; b = t;
}
if(a%b==z)
return b;
return GCD(b, a%b);
}
Bignum EX_GCD(Bignum m, Bignum n, Bignum &x, Bignum &y)
{
Bignum one, zero;
one.init("1");
zero.init("0");
if(n==zero)
{
x = one; y = zero;
return m;
}
Bignum a, a1, b, b1, c, d, q, r, t;
a1.init("1"); b.init("1");
a.init("0"); b1.init("0");
c = m; d = n;
q = c/d; r = c%d;
while(!(r==zero))
{
c = d;
d = r;
t = a1;
a1 = a;
a = t - q * a;
t = b1;
b1 = b;
b = t - q * b;
q = c/d;
r = c%d;
}
x = a; y = b;
return d;
}
/*
Bignum EX_GCD(Bignum a, Bignum b, Bignum &x, Bignum &y)
{
Bignum z, o;
z.init("0");
o.init("1");
if(b==z)
{
x = o; y = z;
return a;
}
Bignum r = EX_GCD(b, a%b, x, y);
Bignum temp = y;
y = x-((a/b)*y);
x = temp;
return r;
}
*/
bool CheckPrime(Bignum a, Bignum n)
{
Bignum one, two, zero;
one.init("1");
two.init("2");
zero.init("0");
Bignum tem = n - one;
int j = 0;
while(tem % two == zero)
{
tem = tem / two;
j++;
}
Bignum x = qpow(a, tem, n);
if(x == one || x == n - one)
return true;
while(j--)
{
x = (x*x)%n;
if(x == n - one)
return true;
}
return false;
}
bool Miller_Robin(Bignum num) //检验num是不是素数 primelen
{
int primelen = num.len;
Bignum one, two, zero, a, b, temp;
one.init("1");
two.init("2");
zero.init("0");
for(int i = 1; i <= 15; i++) //做15次随机检验
{
a.len = primelen-1;
for(int i=0;i<primelen-1;i++)
{
a.s[i] = ((int)rand())%10;
if(i==0)
a.s[i] = ((int)rand())%9+1;
}
if(!CheckPrime(a, num))
return false;
}
return true;
}
Bignum GetPr(int len) //得到一个指定位数的大素数
{
Bignum zero, one, two, res;
one.init("1");
two.init("2");
zero.init("0");
for(int i=0;i<len;i++)
{
res.s[i] = ((int)rand())%10;
if(i==0&&res.s[i]==0)
res.s[i] = ((int)rand())%9+1;
}
res.flag = 0;
res.len = len;
if(res%two==zero) //保证res是奇数
res = res+one;
int done = 0;
while(1)
{
if(Miller_Robin(res))
return res;
res = res + two;
}
}
/*
文件名:
Vectorcal.h
*/
#include <cstdio>
#include <algorithm>
#include <iostream>
#include <string>
#include "VectorDiv.h"
using namespace std;
int FFT_done;
void Reverse(int a[], int len)
{
int i = 0, j = len-1;
while(i<j)
{
swap(a[i], a[j]);
i++; j--;
}
}
void VectorCpy(int a[], int b[], int len)
{
for(int i=0;i<len;i++)
a[i] = b[i];
}
//向量无符号加法
void VectorAdd(int a[], int lena, int b[], int lenb, int c[], int &lenc)
{
Reverse(a, lena);
Reverse(b, lenb);
int i = 0;
int x = 0;
while(i<lena||i<lenb)
{
c[i] = a[i]+b[i]+x;
x = c[i]/10;
c[i] %= 10;
lenc++;
i++;
}
if(x)
{
lenc++;
c[lenc-1] = x;
}
Reverse(c, lenc);
}
//向量无符号减法
void VectorSub(int a[], int lena, int b[], int lenb, int c[], int &lenc)
{
Reverse(a, lena);
Reverse(b, lenb);
int i = 0, flag = 0;
int x = 0;
while(i<lena||i<lenb)
{
if(a[i]<b[i])
{
a[i] += 10;
a[i+1]--;
}
c[i] = a[i]-b[i];
lenc++;
i++;
}
while(c[lenc-1]==0&&lenc>0)
lenc--;
if(lenc==0)
{
c[0] = 0;
lenc = 1;
flag = 0;
}
Reverse(c, lenc);
}
//向量无符号比较大小
int VectorCmp(int a[], int lena, int b[], int lenb)
{
if(lena>lenb)
return 1;
if(lena<lenb)
return 0;
for(int i=0;i<lena;i++)
{
if(a[i]>b[i])
return 1;
else if(a[i]<b[i])
return 0;
}
return 1;
}
void VectorDiv(int a1[], int lena, int b1[], int lenb, int c1[], int &lenc)
{
for(int i=0;i<=max(lena, lenb)*64;i++)
{
a[i] = b[i] = 0;
t3[i] = 0;
bi[i] = 0;
x[i] = y[i] = 0;
}
int n, m, d;
if(!FFT_done)
{
FFT_done = 1;
FFT_Init();
}
std::string str1, str2;
for(int i=0;i<lena;i++)
str1.push_back(a1[i]+'0');
for(int i=0;i<lenb;i++)
str2.push_back(b1[i]+'0');
n = get(a, str1), m = get(b, str2);
d = divide(n,m);
lack(n, m, d);
lenc = d;
for(int i=d-1,p=0;i>=0;i--,p++)
c1[p] = x[i];
}
/*
文件名:
FFT.h
*/
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <cstring>
#define MAXN 5010
#define PI (acos(-1.0))
using namespace std;
struct t_Comp {
double r, i;
t_Comp(double _r = 0.0, double _i = 0.0)
{
r = _r; i = _i;
}
t_Comp operator + (const t_Comp &x) const
{
return t_Comp(r+x.r, i+x.i);
}
t_Comp operator - (const t_Comp &x) const
{
return t_Comp(r-x.r, i-x.i);
}
t_Comp operator * (const t_Comp &x) const
{
return t_Comp(r*x.r-i*x.i, r*x.i+i*x.r);
}
};
char str1[MAXN], str2[MAXN];
t_Comp x1[MAXN], x2[MAXN];
int sum[MAXN];
int rev[MAXN];
void change(t_Comp y[], int len)
{
for(int i=0;i<len;i++)
{
rev[i] = rev[i>>1]>>1;
if(i&1)
rev[i] |= len >> 1;
}
for(int i=0;i<len;i++)
{
if(i<rev[i])
swap(y[i], y[rev[i]]);
}
return;
}
void FFT(t_Comp y[], int len, int flag)
{
change(y, len);
for(int h=2;h<=len;h<<=1)
{
t_Comp wn(cos(2*PI/h), sin(flag*2*PI/h));
for(int j=0;j<len;j+=h)
{
t_Comp w(1, 0);
for(int k=j;k<j+h/2;k++)
{
t_Comp u = y[k];
t_Comp t = w*y[k+h/2];
y[k] = u+t;
y[k+h/2] = u-t;
w = w*wn;
}
}
}
if(flag==-1)
{
for(int i=0;i<len;i++)
y[i].r /= len;
}
}
void VectorMul(int a[], int lena, int b[], int lenb, int c[], int &lenc)
{
int len = 1;
while((len<lena*2)||(len<lenb*2)) //len只能是2的幂次
len <<= 1;
for(int i=0;i<lena;i++)
x1[i] = t_Comp(a[lena-1-i], 0);
for(int i=lena;i<len;i++) //少的位补0
x1[i] = t_Comp(0, 0);
for(int i=0;i<lenb;i++)
x2[i] = t_Comp(b[lenb-1-i], 0);
for(int i=lenb;i<len;i++)
x2[i] = t_Comp(0, 0);
FFT(x1, len, 1); //做DFT
FFT(x2, len, 1);
for(int i=0;i<len;i++)
x1[i] = x1[i]*x2[i];
FFT(x1, len, -1); //做IDFT
for(int i=0;i<len;i++)
c[i] = int(x1[i].r+0.5);;
for(int i=0;i<len;i++) //处理进位
{
c[i+1] += c[i]/10;
c[i] %= 10;
}
len = lena+lenb-1;
while(c[len]==0&&len>0) //去除前导0
len--;
int i = 0, j = len;
while(i<j)
{
swap(c[i], c[j]);
i++; j--;
}
lenc = len+1;
}
/*
int main()
{
int a[210] = {0}, b[210] = {0}, c[210] = {0};
int lenr = 0;
scanf("%s%s", str1, str2);
int len1 = strlen(str1), len2 = strlen(str2);
for(int i=0;i<len1;i++)
a[i] = str1[i]-'0';
for(int i=0;i<len2;i++)
b[i] = str2[i]-'0';
VectorMul(a, len1, b, len2, c, lenr);
for(int i=0;i<lenr;i++)
printf("%d", c[i]);
printf("\n%d\n", lenr);
system("pause");
return 0;
}*/
/*
文件名: VectorDiv.h
*/
#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#define inline __inline__ __attribute__((always_inline))
#define L 524288
//#define L 10010
typedef long long LL;
int REV[L<<1];
struct complex
{
double a,b;
inline void operator+=(const complex &z){a += z.a, b += z.b;}
inline complex operator-(const complex &z){return {a-z.a,b-z.b};}
inline complex operator*(const complex &z){return {a*z.a-b*z.b,a*z.b+b*z.a};}
}W[2][L];
void FFT_Init()
{
int i, j;
const int l = L>>1;
const double pi = acos(-1);
for(i=j=1;j<L;j<<=1)
for(;i<j<<1;i++)
REV[i<<1]=REV[i],REV[i<<1|1]=REV[i]+j;
for(i=0;i<l;i++)
W[0][i+l] = {cos(pi*i/l), sin(pi*i/l)};
for(i=l-1;i;i--)
W[0][i] = W[0][i<<1];
memcpy(W[1],W[0],sizeof W[1]);
for(i=1;i<L;i++)
W[1][i].b *= -1;
}
void FFT(complex *f,int len,int sign)
{
int i = 1,j,k; complex *p,*q,*w,t; sign = (sign < 0);
for(int *r=REV+len+1;i<len;i++,r++)if(i<(k=*r))t=f[i],f[i]=f[k],f[k]=t;
for(i=1;i<len;i<<=1)for(j=0;j<len;j+=i,j+=i)
for(q=(p=f+j)+i,w=W[sign]+i,k=0;k<i;k++,p++,q++)t=*q**w++,*q=*p-t,*p+=t;
if(sign){double p=1./len;for(i=0;i<len;i++,f++)f->a *= p, f->b *= p;}
}
int get(char *s, std::string str)
{
int c, l = 0;
do
{
*s++ = char(str[l]^48);
l++;
}
while(str[l] >= '0' && str[l] <= '9');
return l;
}
char a[L],b[L],bi[L],x[L],y[L];
complex t1[L],t2[L];
LL t3[L+1];
void Newton(int len) //计算b的倒数，有效位数为len，结果存进bi里
{
bool g = 1;
int l = 16, l2 = 32, l4 = 64, i;
long double d = 0, e = 1;
for(i=0;i<20;i++)
d = d+e*b[i], e *= 0.1;
d = 10./d;
if(d < 10)
{
for(i=0;i<=l;i++)
{
bi[i] = d, d = (d-bi[i])*10;
}
bi[l-1] += (bi[l] > 4);
}
else bi[0] = 10;
while(l < len)
{
p:
memset(t1, 0, sizeof(complex)*l4);
memset(t2, 0, sizeof(complex)*l4);
for(i=0;i<l2;i++)t1[i].a = b[i];
for(i=0;i<l;i++)t2[i].a = bi[i];
FFT(t1,l4,1), FFT(t2,l4,1);
for(i=0;i<l4;i++)
{
t1[i] = t1[i]*t2[i];
t1[i].a = 20-t1[i].a, t1[i].b = -t1[i].b;
t2[i] = t1[i]*t2[i];
}
FFT(t2,l4,-1); t3[l4] = 0;
for(i=l4-1;i>=0;i--)
{
t3[i] = LL(floor(t2[i].a+0.5))+t3[i+1]/10, t3[i+1] %= 10;
if(t3[i+1] < 0)t3[i+1] += 10, t3[i]--;
}
if(t3[0] > 9)
{
bi[0] = t3[0]/10, t3[0] %= 10;
for(i=1;i<l2;i++)bi[i] = t3[i-1];
bi[l2-1] += (t3[l2-1] > 4);
}
else
{
for(i=0;i<l2;i++)bi[i] = t3[i];
bi[l2-1] += (t3[l2] > 4);
}
l <<= 1, l2 <<= 1, l4 <<= 1;
}
if(g)
{
g = 0, l >>= 1, l2 >>= 1, l4 >>= 1;
goto p;
}
}
int divide(int n, int m) //计算a/b，整数部分存进x里
{
int p = n-m+16,l,i;
Newton(p);
for(l=1;l<n+p;l<<=1);
memset(t1,0,sizeof(complex)*l);
memset(t2,0,sizeof(complex)*l);
for(i=0;i<n;i++)t1[i].a = a[i];
for(i=0;i<p;i++)t2[i].a = bi[i];
FFT(t1,l,1), FFT(t2,l,1);
for(i=0;i<l;i++)t1[i] = t1[i]*t2[i];
FFT(t1,l,-1); t3[l] = 0;
for(i=l-1;i>=0;i--)t3[i] = LL(t1[i].a+0.5)+t3[i+1]/10, t3[i+1] %= 10;
if(t3[0] > 9)
{
x[0] = t3[0]/10, t3[0] %= 10, l = n-m+1;
for(i=0;i<n-m;i++)x[i+1] = t3[i];
}
else for(l=n-m,i=0;i<n-m;i++)x[i] = t3[i];
return l;
}
void lack(int n,int m,int &d) //微调
{
int l,i,j; char t; LL tl = 0;
for(i=0,j=n-1;i<j;i++,j--)
{
t = a[i], a[i] = a[j], a[j] = t;
}
for(i=0,j=m-1;i<j;i++,j--)
{
t = b[i], b[i] = b[j], b[j] = t;
}
for(i=0,j=d-1;i<j;i++,j--)
{
t = x[i], x[i] = x[j], x[j] = t;
}
for(l=1;l<n;l<<=1);
memset(t1,0,sizeof(complex)*l);
memset(t2,0,sizeof(complex)*l);
for(i=0;i<m;i++)
t1[i].a = b[i];
for(i=0;i<d;i++)
t2[i].a = x[i];
t2[0].a += 1.;
FFT(t1, l, 1);
FFT(t2, l, 1);
for(i=0;i<l;i++)
t1[i] = t1[i]*t2[i];
FFT(t1, l, -1);
for(i=0;i<=n;i++)
{
tl += LL(t1[i].a+0.5), y[i] = tl%10, tl /= 10;
}
for(i=n;i>=0;i--)
{
{
if(y[i] > a[i])
return;
else if(y[i] < a[i])
break;
}
}
for(x[0]++,i=0;x[i]>9;i++)
{
x[i+1] += x[i]/10;
x[i] %= 10;
}
if(x[d])
d++;
}
/*
int main()
{
int n, m, d;
FFT_Init();
std::string str1, str2;
std::cin>>str1>>str2;
n = get(a, str1), m = get(b, str2);
d = divide(n,m);
lack(n,m,d);
//NTT::getout(x,d);
system("pause");
return 0;
}
*/
版权声明：本文为Berserker____原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。