8次FFT50分求助(WA)

P4245 【模板】任意模数多项式乘法

wtyqwq @ 2020-03-23 19:55:43

\text{double} 都改成 \text{long} \text{double} 了,\text{int} 全部改为了 \text{long} \text{long} 都没有用。

#pragma GCC optimize(3)
#pragma GCC optimize("Ofast")
#include <cstdio>
#include <cmath>
#include <cctype>
#include <iomanip>
#define R register
#define N (1 << 17)
#define int long long
#define ll int
#define double long double
#define PI 3.141592653589793238462643383279502884
using namespace std;
char Buf[1 << 24], *S = Buf, *T = Buf;
//#define getchar() (S == T && (T = (S = Buf) + fread(Buf, 1, 1 << 24, stdin), S == T) ? EOF : *S++)
inline int input() {
    R int x(0), f(0);
    R char c(getchar());
    while (!isdigit(c)) f |= (c == '-'), c = getchar();
    while (isdigit(c)) x = (x << 1) + (x << 3) + (c ^ 48), c = getchar();
    return f ? -x : x;
}
struct complex {
    double x, y;
    inline complex(double x = 0, double y = 0) : x(x), y(y) {}
    inline complex operator+(const complex b) const { return complex(x + b.x, y + b.y); }
    inline complex operator-(const complex b) const { return complex(x - b.x, y - b.y); }
    inline complex operator*(const complex b) const { return complex(x * b.x - y * b.y, x * b.y + y * b.x); }
};
int n, m, t, lim, k, mod;
int rev[N << 2], ans[N << 1];
complex a[N << 2], b[N << 2], c[N << 2], d[N << 2], X[N << 2];
inline void FFT(complex *f, int flag) {
    for (int i = 0; i < lim; ++i)
        if (i < rev[i])
            swap(f[i], f[rev[i]]);
    for (R int l = 2; l <= lim; l <<= 1) {
        int len = l >> 1;
        R complex OMG(cos(PI / len), sin(PI / len));
        OMG.y *= flag;
        for (R int k = 0; k < lim; k += l) {
            R complex buf(1, 0);
            for (R int i = k; i < k + len; ++i) {
                complex tmp = buf * f[len + i];
                f[len + i] = f[i] - tmp, f[i] = f[i] + tmp, buf = buf * OMG;
            }
        }
    }
}
inline double Round(double x) {
    if (x == 0)
        return x;
    x > 0 ? x += 0.5 : x -= 0.5;
    return x;
}
inline void solve(complex *a, complex *b, int res) {
    for (R int i = 0; i < lim; ++i) X[i] = a[i] * b[i];
    FFT(X, -1);
    //for (R int i = 0; i < lim; ++i) printf("%.5lf\n", X[i]);
    for (R int i = 0; i <= t; ++i) ans[i] = (ll)((ll)ans[i] + (ll)Round(X[i].x / lim) % (ll)mod * (ll)res % mod) % mod;
}
inline void MTT(complex *a, complex *b, complex *c, complex *d) {
    FFT(a, 1), FFT(b, 1), FFT(c, 1), FFT(d, 1);
    //for (int i = 0; i < lim; ++i) printf("%.4lf\n", a[i]);
    //for (int i = 0; i < lim; ++i) printf("%.4lf\n", b[i]);
    //for (int i = 0; i < lim; ++i) printf("%.4lf\n", c[i]);
    //for (int i = 0; i < lim; ++i) printf("%.4lf\n", d[i]);
    solve(a, c, (1ll << 30) % mod), solve(a, d, (1 << 15) % mod);
    solve(c, b, (1ll << 15) % mod), solve(b, d, 1);
}
signed main() {
    n = input(), m = input(), mod = input();
    for (int i = 0; i <= n; ++i) {
        int x = input();
        a[i].x = (double)(x >> 15ll);
        c[i].x = (double)(x & ((1ll << 15ll) - 1ll));
        //printf("%.4lf %.4lf\n", a[i].x, c[i].x);
    }
    for (int i = 0; i <= m; ++i) {
        int x = input();
        b[i].x = (double)(x >> 15ll);
        d[i].x = (double)(x & ((1ll << 15) - 1ll));
    }
    for (t = n + m, lim = 1ll; lim <= t;) lim <<= 1, ++k;
    for (int i = 0ll; i < lim; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (k - 1));
    MTT(a, c, b, d);
    for (int i = 0ll; i <= t; ++i) printf("%lld ", ans[i]);
    return 0;
}

by Smile_Cindy @ 2020-03-23 20:15:39

@Oler_Accepted 我7次DFT也没掉精度,看一下吧

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const long double PI=acos(-1);
const int MAX_N=1<<18;
const int base=1<<15;
int n,m,Mod; 
struct Complex
{
    long double x,y;
    Complex(){}
    Complex(long double _x,long double _y)
    {
        x=_x;
        y=_y;
    }
};
Complex operator + (Complex w1,Complex w2)
{
    return Complex(w1.x+w2.x,w1.y+w2.y);
}
Complex operator - (Complex w1,Complex w2)
{
    return Complex(w1.x-w2.x,w1.y-w2.y);
}
Complex operator * (Complex w1,Complex w2)
{
    return Complex(w1.x*w2.x-w1.y*w2.y,w1.x*w2.y+w1.y*w2.x);
}
int rev[MAX_N];
void FFT(Complex A[],int n,int flag)
{
    for(int i=0;i<n;i++)if(rev[i]>i)swap(A[i],A[rev[i]]);
    for(int i=1;i<n;i<<=1)
    {
        Complex w1=Complex(cos(2*PI*flag/(i<<1)),sin(2*PI*flag/(i<<1)));
        for(int j=0;j<n;j+=(i<<1))
        {
            Complex w2=Complex(1,0);
            for(int k=0;k<i;k++)
            {
                Complex t1=A[j+k],t2=A[i+j+k]*w2;
                A[j+k]=t1+t2;
                A[i+j+k]=t1-t2;
                w2=w2*w1;
            }
        }
    }
    if(flag==-1)for(int i=0;i<n;i++)A[i].x/=n;
}
int len=1,L=0;
void init()
{
    while(len<=n+m)
    {
        len<<=1;
        L++;
    }
    for(int i=1;i<len;i++)rev[i]=rev[i>>1]>>1|(i&1)<<(L-1);
}
Complex A1[MAX_N],A2[MAX_N],B1[MAX_N],B2[MAX_N];
Complex A[MAX_N],B[MAX_N],C[MAX_N];
int F[MAX_N],G[MAX_N];
int ans[MAX_N];
signed main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout.tie(0);
    cin>>n>>m>>Mod;
    for(int i=0;i<=n;i++)cin>>F[i];
    for(int i=0;i<=m;i++)cin>>G[i];
    init();
    for(int i=0;i<=n;i++)A1[i].x=F[i]/base,B1[i].x=F[i]%base;
    for(int i=0;i<=m;i++)A2[i].x=G[i]/base,B2[i].x=G[i]%base;
    FFT(A1,len,1);
    FFT(A2,len,1);
    FFT(B1,len,1);
    FFT(B2,len,1);
    for(int i=0;i<len;i++)A[i]=A1[i]*A2[i];
    for(int i=0;i<len;i++)B[i]=A1[i]*B2[i]+A2[i]*B1[i];
    for(int i=0;i<len;i++)C[i]=B1[i]*B2[i];
    FFT(A,len,-1);
    FFT(B,len,-1);
    FFT(C,len,-1);
    for(int i=0;i<len;i++)
    {
        ans[i]=(LL)(A[i].x+0.5)%Mod;
        ans[i]=(1LL*ans[i]*base+(LL)(B[i].x+0.5))%Mod;
        ans[i]=(1LL*ans[i]*base+(LL)(C[i].x+0.5))%Mod;
    }
    for(int i=0;i<=n+m;i++)cout<<ans[i]<<" \n"[i==n+m];
    return 0;
}

by wtyqwq @ 2020-03-23 20:16:14

金勾,对不起


by wtyqwq @ 2020-03-23 20:16:52

@Alpha Thanks


by 行吟啸九州 @ 2020-03-23 20:24:52

@刘兆洲 那个做法,MYY自己的代码都掉精,更过不去


by wtyqwq @ 2020-03-23 21:01:38

@Alpha 竟然是 \pi 的问题。。。必须用 acos(-1)。。。


by Smile_Cindy @ 2020-03-23 21:02:47

@Oler_Accepted 一般人都用acos(-1),因为这样的话不论精度多少都能用,否则PI至少输20位。


by wtyqwq @ 2020-03-23 21:04:24

@Alpha 嗯,好的


by iostream @ 2020-03-23 21:08:59

@Oler_Accepted 其实可以用 double 的要不然可能会太慢,我的记录。


上一页 |