3123 高精度练习之超大整数乘法

题目描述 Description

    给出两个正整数A和B,计算A*B的值。保证A和B的位数不超过100000位。

输入描述 Input Description

    读入两个用空格隔开的正整数

输出描述 Output Description

    输出A*B的值

样例输入 Sample Input

    4 9

样例输出 Sample Output

    36

数据范围及提示 Data Size & Hint

    两个正整数的位数不超过100000位

什么都不说了,裸FFT

看了一天的算导,勉强看懂了怎么O(nlogn)求dft,求dft^(-1)就实在看不懂了,唉,只好记下结论

首先是理解系数表达和点值表达,因为点值表达乘法是O(n)的,所以我们要通过系数表达算出点值表达,相乘后再算出点值表达

这个是分治求的dft奇偶分成两个序列,分别求出这两个序列的dft然后通过下面的代码就可以O(n)合并这个两个dft了

1 for i:=0 to n>>(t+1)-1 do
2   begin
3     p:=i<<(t+1)+s;
4     wt:=w[i<<t]*a[p+1<<t];
5     tt[i]:=a[p]+wt;
6     tt[i+n>>(t+1)]:=a[p]-wt;
7   end;

先把单位根的几个性质搞懂,然后就基本上可以理解这个分治了(最好是自己手算一下......我就懒得说什么了,自己都不是很理解,代码还是抄别人的)

zky神犇手动搬运算导FFT

pascal代码(我从这里抄的)

  1 const
  2         s=4;
  3 type
  4         cp=record
  5           x,y:double;
  6         end;
  7         arr=array[0..1 shl 17]of cp;
  8 var
  9         c:array[0..100000]of int64;
 10         a,b,w,tt:arr;
 11         n,i,p:longint;
 12         st:ansistring;
 13         wt:cp;
 14 
 15 operator *(var a,b:cp)c:cp;
 16 begin
 17         c.x:=a.x*b.x-a.y*b.y;
 18         c.y:=a.x*b.y+a.y*b.x;
 19 end;
 20 
 21 operator +(var a,b:cp)c:cp;
 22 begin
 23         c.x:=a.x+b.x;
 24         c.y:=a.y+b.y;
 25 end;
 26 
 27 operator -(var a,b:cp)c:cp;
 28 begin
 29         c.x:=a.x-b.x;
 30         c.y:=a.y-b.y;
 31 end;
 32 
 33 procedure dft(var a:arr;s,t:longint);
 34 begin
 35         if n>>t=1 then exit;
 36         dft(a,s,t+1);
 37         dft(a,s+1<<t,t+1);
 38         for i:=0 to n>>(t+1)-1 do
 39           begin
 40             p:=i<<(t+1)+s;
 41             wt:=w[i<<t]*a[p+1<<t];
 42             tt[i]:=a[p]+wt;
 43             tt[i+n>>(t+1)]:=a[p]-wt;
 44           end;
 45         for i:=0 to n>>t-1 do
 46           a[i<<t+s]:=tt[i];
 47 end;
 48 
 49 procedure get(var a:arr);
 50 var
 51         i,l,ll:longint;
 52         c:char;
 53 begin
 54         read(c);
 55         st:='';
 56         while (ord(c)>=ord('0')) and (ord(c)<=ord('9')) do
 57           begin
 58             st:=st+c;
 59             read(c);
 60           end;
 61         while length(st)mod s<>0 do
 62           insert('0',st,0);
 63         ll:=length(st)div s;
 64         for l:=1 to ll do
 65           val(copy(st,l*s-s+1,s),a[ll-l].x);
 66         while n>>1<l do
 67           n:=n<<1;
 68 end;
 69 
 70 begin
 71         n:=1;
 72         get(a);
 73         get(b);
 74         for i:=0 to n-1 do
 75           w[i].x:=cos(pi*2*i/n);
 76         for i:=0 to n-1 do
 77           w[i].y:=sin(pi*2*i/n);
 78         dft(a,0,0);
 79         dft(b,0,0);
 80         for i:=0 to n-1 do
 81           a[i]:=a[i]*b[i];
 82         for i:=0 to n-1 do
 83           w[i].y:=-w[i].y;
 84         dft(a,0,0);
 85         for i:=0 to n-1 do
 86           begin
 87             c[i]:=c[i]+round(a[i].x/n);
 88             c[i+1]:=c[i] div 10000;
 89             c[i]:=c[i] mod 10000;
 90           end;
 91         while (c[i]=0) and (i>0) do
 92           dec(i);
 93         for p:=i downto 0 do
 94           begin
 95             str(c[p],st);
 96             while (i<>p) and (length(st)<s) do
 97               st:='0'+st;
 98             write(st);
 99           end;
100 end.
View Code

又写了一遍233

 1 const
 2     maxn=400400;
 3 type
 4     cp=record
 5         x,y:double;
 6     end;
 7     arr=array[0..maxn]of cp;
 8 var
 9     a,b,w,tt:arr;
10     n:longint;
11     s:ansistring;
12     c:array[0..maxn]of longint;
13 
14 operator -(a,b:cp)c:cp;
15 begin
16     c.x:=a.x-b.x;
17     c.y:=a.y-b.y;
18 end;
19 
20 operator +(a,b:cp)c:cp;
21 begin
22     c.x:=a.x+b.x;
23     c.y:=a.y+b.y;
24 end;
25 
26 operator *(a,b:cp)c:cp;
27 begin
28     c.x:=a.x*b.x-a.y*b.y;
29     c.y:=a.x*b.y+a.y*b.x;
30 end;
31 
32 procedure get(var a:arr);
33 var
34     c:char;
35     i,len:longint;
36 begin
37     s:='';
38     read(c);
39     while c in['0'..'9'] do
40         begin
41             s:=s+c;
42             read(c);
43         end;
44     len:=length(s);
45     for i:=1 to len do
46         a[len-i].x:=ord(s[i])-ord('0');
47     while n<len<<1 do
48         n:=n<<1;
49 end;
50 
51 procedure dft(var a:arr;s,t:longint);
52 var
53     i,p:longint;
54     wt:cp;
55 begin
56     if n>>t=1 then exit;
57     dft(a,s+1<<t,t+1);dft(a,s,t+1);
58     for i:=0 to n>>t>>1-1 do
59         begin
60             p:=i<<t<<1+s;
61             wt:=w[i<<t]*a[p+1<<t];
62             tt[i]:=a[p]+wt;
63             tt[i+n>>t>>1]:=a[p]-wt;
64         end;
65     for i:=0 to n>>t-1 do
66         a[s+i<<t]:=tt[i];
67 end;
68 
69 procedure main;
70 var
71     i:longint;
72 begin
73     n:=1;get(a);get(b);
74     for i:=0 to n-1 do w[i].x:=cos(pi*2*i/n);
75     for i:=0 to n-1 do w[i].y:=sin(pi*2*i/n);
76     dft(a,0,0);dft(b,0,0);
77     for i:=0 to n-1 do a[i]:=a[i]*b[i];
78     for i:=0 to n-1 do w[i].y:=-w[i].y;
79     dft(a,0,0);
80     for i:=0 to n-1 do c[i]:=round(a[i].x/n);
81     for i:=0 to n-2 do
82         begin
83             inc(c[i+1],c[i]div 10);
84             c[i]:=c[i]mod 10;
85         end;
86     i:=n-1;
87     while (c[i]=0) and (i>0) do dec(i);
88     for i:=i downto 0 do write(c[i]);
89 end;
90 
91 begin
92     main;
93 end.
View Code
原文地址:https://www.cnblogs.com/Randolph87/p/3669526.html