2014-team3/code/nt

从 Trac 迁移的文章

这是从旧校内 Wiki 迁移的文章，可能存在一些样式问题，您可以向 memset0 反馈。

原文章内容如下：

{{{
#include<iostream>
#include<vector>
#include<algorithm>
#include<cassert>
using namespace std;
typedef unsigned long long num;
typedef long long int64;
num pm(num a,num n,num m){
    num r=1;
    for(;n;n>>=1,a=a*a%m)
        if(n&1)r=r*a%m;
    return r;
}
void ex(int64 a,int64 b,int64& x,int64& y){
    if(b==0)x=1,y=0;
    else {ex(b,a%b,y,x);y-=a/b*x;}
}
bool mr(num a,num n){ // probably prime
    int s=__builtin_ctzll(n-1),i;
    num d=n>>s,x=pm(a%n,d,n);
    for(i=0;i<s&&x!=1&&x!=n-1;++i)
        x=x*x%n;
    return x==n-1||i==0;
}
bool isprime(num n){
    if(n<=4)return n&2;  
    static const num sp[]={2,3,5,7,11,13,17,19,23,29,31,37};
    static const num s7[]={2,325,9375,28178,450775,9780504, 1795265022}; // at least 2^64
    for(int i=0;i<12;++i)
        if(n%sp[i]==0)return n==sp[i];
    if(n<4759123141ull) 
        return mr(2,n)&&mr(7,n) && mr(61,n);
    for(int i=0;i<7;++i)
        if(!mr(s7[i],n))return 0;
    return 1;
}
num Tonelli_Shanks(num n,num p){
//p is an odd prime  &&1<n<p && pm(n,p>>1,p)==1 
    if(p&2)return pm(n,p+1>>2,p);
    int s=__builtin_ctzll(p^1);
    num q = p>>s,z=2;
    for(;pm(z,p>>1,p)==1;++z);
    num c=pm(z,q,p),r=pm(n,q+1>>1,p),t=pm(n,q,p),tmp;
    for(int m=s,i;t!=1;){
        for(i=0,tmp=t;tmp!=1;++i)tmp=tmp*tmp%p;
        for(;i<--m;)c=c*c%p;
        r = r * c %p;
        c = c * c %p;
        t = t * c %p;
    }
    return r;
}
num Cipolla(num n,num p){
    num a=0;
    for(;pm((a*a+p-n)%p,p>>1,p)==1;++a);
    num I=(a*a+p-n)%p,rx=1,ry=0,tx=a,ty=1;
    auto mult=[=](num&a,num&b,num c,num d){
        num tmp=(b*d%p*I+a*c)%p;
        b=(a*d+b*c)%p;a=tmp;
    };
    for(num c=p+1>>1;c;c>>=1,mult(tx,ty,tx,ty))
        if(c&1)mult(rx,ry,tx,ty);
    return rx;
}

//2113929217,5
//1811939329,13
//2130706433,3
template<num P=2013265921,num g=31>
void ntt(num *a,size_t n,bool inv=false){
    num w=1,d=pm(g,(P-1)/n,P),t;
    size_t i,j,c,s;
    if(inv){
        for(i=1,j=n-1;i<j;swap(a[i++],a[j--]));
        for(t=pm(n,P-2,P),i=0;i<n;++i)a[i]=a[i]*t%P;
    }
    for(s=n>>1;s;s>>=w=1,d=d*d%P)
    for(c=0;c<s;++c,w=w*d%P)
    for(i=c;i<n;i+=s<<1){
        a[i|s]=(a[i]+P-(t=a[i|s]))*w%P;
        a[i]=(a[i]+t)%P;
    }
    for(i=1;i<n;++i){
        for(j=0,s=i,c=n>>1;c;c>>=1,s>>=1)j=j<<1|s&1;
        if(i<j)swap(a[i],a[j]);
    }
}
}}}

#include<iostream>
#include<vector>
#include<algorithm>
#include<cassert>
using namespace std;
typedef unsigned long long num;
typedef long long int64;
num pm(num a,num n,num m){
    num r=1;
    for(;n;n>>=1,a=a*a%m)
        if(n&1)r=r*a%m;
    return r;
}
void ex(int64 a,int64 b,int64& x,int64& y){
    if(b==0)x=1,y=0;
    else {ex(b,a%b,y,x);y-=a/b*x;}
}
bool mr(num a,num n){ // probably prime
    int s=__builtin_ctzll(n-1),i;
    num d=n>>s,x=pm(a%n,d,n);
    for(i=0;i<s&&x!=1&&x!=n-1;++i)
        x=x*x%n;
    return x==n-1||i==0;
}
bool isprime(num n){
    if(n<=4)return n&2;  
    static const num sp[]={2,3,5,7,11,13,17,19,23,29,31,37};
    static const num s7[]={2,325,9375,28178,450775,9780504, 1795265022}; // at least 2^64
    for(int i=0;i<12;++i)
        if(n%sp[i]==0)return n==sp[i];
    if(n<4759123141ull) 
        return mr(2,n)&&mr(7,n) && mr(61,n);
    for(int i=0;i<7;++i)
        if(!mr(s7[i],n))return 0;
    return 1;
}
num Tonelli_Shanks(num n,num p){
//p is an odd prime  &&1<n<p && pm(n,p>>1,p)==1 
    if(p&2)return pm(n,p+1>>2,p);
    int s=__builtin_ctzll(p^1);
    num q = p>>s,z=2;
    for(;pm(z,p>>1,p)==1;++z);
    num c=pm(z,q,p),r=pm(n,q+1>>1,p),t=pm(n,q,p),tmp;
    for(int m=s,i;t!=1;){
        for(i=0,tmp=t;tmp!=1;++i)tmp=tmp*tmp%p;
        for(;i<--m;)c=c*c%p;
        r = r * c %p;
        c = c * c %p;
        t = t * c %p;
    }
    return r;
}
num Cipolla(num n,num p){
    num a=0;
    for(;pm((a*a+p-n)%p,p>>1,p)==1;++a);
    num I=(a*a+p-n)%p,rx=1,ry=0,tx=a,ty=1;
    auto mult=[=](num&a,num&b,num c,num d){
        num tmp=(b*d%p*I+a*c)%p;
        b=(a*d+b*c)%p;a=tmp;
    };
    for(num c=p+1>>1;c;c>>=1,mult(tx,ty,tx,ty))
        if(c&1)mult(rx,ry,tx,ty);
    return rx;
}
//2113929217,5
//1811939329,13
//2130706433,3
template<num P=2013265921,num g=31>
void ntt(num *a,size_t n,bool inv=false){
    num w=1,d=pm(g,(P-1)/n,P),t;
    size_t i,j,c,s;
    if(inv){
        for(i=1,j=n-1;i<j;swap(a[i++],a[j--]));
        for(t=pm(n,P-2,P),i=0;i<n;++i)a[i]=a[i]*t%P;
    }
    for(s=n>>1;s;s>>=w=1,d=d*d%P)
    for(c=0;c<s;++c,w=w*d%P)
    for(i=c;i<n;i+=s<<1){
        a[i|s]=(a[i]+P-(t=a[i|s]))*w%P;
        a[i]=(a[i]+t)%P;
    }
    for(i=1;i<n;++i){
        for(j=0,s=i,c=n>>1;c;c>>=1,s>>=1)j=j<<1|s&1;
        if(i<j)swap(a[i],a[j]);
    }
}