zoukankan      html  css  js  c++  java
  • 1402大数相乘 FFT算法学习!

    网上找的,慢慢研究下

    View Code
      1 #include <cmath>  
    2 #include <complex>
    3 #include <cstring>
    4 using namespace std;
    5
    6 const double PI = acos(-1);
    7 typedef complex<double> cp;
    8 typedef long long int64;
    9
    10 const int N = 1 << 16;
    11 int64 a[N], b[N], c[N << 1];
    12
    13 void bit_reverse_copy(cp a[], int n, cp b[])
    14 {
    15 int i, j, k, u, m;
    16 for (u = 1, m = 0; u < n; u <<= 1, ++m);
    17 for (i = 0; i < n; ++i)
    18 {
    19 j = i; k = 0;
    20 for (u = 0; u < m; ++u, j >>= 1)
    21 k = (k << 1) | (j & 1);
    22 b[k] = a[i];
    23 }
    24 }
    25
    26 void FFT(cp _x[], int n, bool flag)
    27 {
    28 static cp x[N << 1];
    29 bit_reverse_copy(_x, n, x);
    30 int i, j, k, kk, p, m;
    31 for (i = 1, m = 0; i < n; i <<= 1, ++m);
    32 double alpha = 2 * PI;
    33 if (flag) alpha = -alpha;
    34 for (i = 0, k = 2; i < m; ++i, k <<= 1)
    35 {
    36 cp wn = cp(cos(alpha / k), sin(alpha / k));
    37 for (j = 0; j < n; j += k)
    38 {
    39 cp w = 1, u, t;
    40 kk = k >> 1;
    41 for (p = 0; p < kk; ++p)
    42 {
    43 t = w * x[j + p + kk];
    44 u = x[j + p];
    45 x[j + p] = u + t;
    46 x[j + p + kk] = u - t;
    47 w *= wn;
    48 }
    49 }
    50 }
    51 memcpy(_x, x, sizeof(cp) * n);
    52 }
    53
    54 void polynomial_multiply(int64 a[], int na, int64 b[], int nb, int64 c[], int &nc)
    55 {
    56 int i, n;
    57 i = max(na, nb) << 1;
    58 for (n = 1; n < i; n <<= 1);
    59 static cp x[N << 1], y[N << 1];
    60 for (i = 0; i < na; ++i) x[i] = a[i];
    61 for (; i < n; ++i) x[i] = 0;
    62 FFT(x, n, 0);
    63 for (i = 0; i < nb; ++i) y[i] = b[i];
    64 for (; i < n; ++i) y[i] = 0;
    65 FFT(y, n, 0);
    66 for (i = 0; i < n; ++i) x[i] *= y[i];
    67 FFT(x, n, 1);
    68 for (i = 0; i < n; ++i)
    69 {
    70 c[i] =(int64)(x[i].real() / n + 0.5);
    71 }
    72 for (nc = na + nb - 1; nc > 1 && !c[nc - 1]; --nc);
    73 }
    74
    75 const int LEN = 5, MOD = 100000;
    76 void convert(char *s, int64 a[], int &n)
    77 {
    78 int len = strlen(s), i, j, k;
    79 for (n = 0, i = len - LEN; i >= 0; i -= LEN)
    80 {
    81 for (j = k = 0; j < LEN; ++j)
    82 k = k * 10 + (s[i + j] - '0');
    83 a[n++] = k;
    84 }
    85 i += LEN;
    86 if (i)
    87 {
    88 for (j = k = 0; j < i; ++j)
    89 k = k * 10 + (s[j] - '0');
    90 a[n++] = k;
    91 }
    92 }
    93
    94 void print(int64 a[], int n)
    95 {
    96 printf("%I64d", a[--n]);
    97 while (n) printf("%05I64d", a[--n]);
    98 puts("");
    99 }
    100
    101 char buf[N + 10];
    102
    103 int main()
    104 {
    105 int na, nb, nc;
    106
    107 while (scanf("%s", buf) != EOF)
    108 {
    109 bool sign = false;
    110 if (buf[0] == '-')
    111 {
    112 sign = !sign;
    113 convert(buf + 1, a, na);
    114 }
    115 else convert(buf, a, na);
    116 scanf("%s", buf);
    117 if (buf[0] == '-')
    118 {
    119 sign = !sign;
    120 convert(buf + 1, b, nb);
    121 }
    122 else convert(buf, b, nb);
    123 polynomial_multiply(a, na, b, nb, c, nc);
    124 int64 t1, t2;
    125 t1 = 0;
    126 for (int i = 0; i < nc; ++i)
    127 {
    128 t2 = t1 + c[i];
    129 c[i] = t2 % MOD;
    130 t1 = t2 / MOD;
    131 }
    132 for (; t1; t1 /= MOD) c[nc++] = t1 % MOD;
    133 if (sign) putchar('-');
    134 print(c, nc);
    135 }
    136 return 0;
    137 }
  • 相关阅读:
    OSI安全体系结构
    PHP 二维数组根据相同的值进行合并
    Java实现 LeetCode 17 电话号码的字母组合
    Java实现 LeetCode 16 最接近的三数之和
    Java实现 LeetCode 16 最接近的三数之和
    Java实现 LeetCode 16 最接近的三数之和
    Java实现 LeetCode 15 三数之和
    Java实现 LeetCode 15 三数之和
    Java实现 LeetCode 15 三数之和
    Java实现 LeetCode 14 最长公共前缀
  • 原文地址:https://www.cnblogs.com/bersaty/p/2281186.html
Copyright © 2011-2022 走看看