• 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
  • 相关阅读:
    理解HTTP幂等性
    企业技术树
    数据库MySQL-Oracle-DB2-SQLServer分页查询
    Redis安装教程
    Redis VS Memcached
    Redis简介
    Redis系列文章导读
    坐标轴
    图例
    画网格
  • 原文地址:https://www.cnblogs.com/Randolph87/p/3669526.html
Copyright © 2020-2023  润新知