math.tex

\section{math}

\subsection{prime/phi/mu}

\begin{lstlisting}
void pj() {
    int n = 10000000;
    pr = 0;
    mu[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (v[i] == 0) {
            v[i] = i;
            prime[pr++] = i;
            phi[i] = i - 1;
            mu[i] = -1;
        }
        for (int j = 0; j < pr && n / i >= prime[j]; j++) {
            v[i * prime[j]] = prime[j];
            if (v[i] == prime[j])
                mu[i * prime[j]] = 0;
            else
                mu[i * prime[j]] = -mu[i];
            phi[i*prime[j]] = phi[i] * (i%prime[j]?prime[j]-1:prime[j]);
            if (v[i] <= prime[j]) break;
        }
    }
}
\end{lstlisting}

\subsection{Congruence Equation}

\begin{lstlisting}
int ex_gcd(int a, int b, int& x, int& y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    int d = ex_gcd(b, a % b, x, y);
    int temp = x;
    x = y;
    y = temp - a / b * y;
    return d;
}
bool liEu(int a, int b, int c, int& x, int& y) {
    int d = ex_gcd(a, b, x, y);
    if (c % d != 0) return 0;
    int k = c / d;
    x *= k;
    y *= k;
    return 1;
}
\end{lstlisting}

\subsection{CRT}

\begin{lstlisting}
lt ai[maxn],bi[maxn];

lt mul(lt a,lt b,lt mod)
{
    lt res=0;
    while(b>0)
    {
        if(b&1) res=(res+a)%mod;
        a=(a+a)%mod;
        b>>=1;
    }
    return res;
}

lt exgcd(lt a,lt b,lt &x,lt &y)
{
    if(b==0){x=1;y=0;return a;}
    lt gcd=exgcd(b,a%b,x,y);
    lt tp=x;
    x=y; y=tp-a/b*y;
    return gcd;
}

lt excrt()
{
    lt x,y,k;
    lt M=bi[1],ans=ai[1];//第一个方程的解特判
    for(int i=2;i<=n;i++)
    {
        lt a=M,b=bi[i],c=(ai[i]-ans%b+b)%b;//ax≡c(mod b)
        lt gcd=exgcd(a,b,x,y),bg=b/gcd;
        if(c%gcd!=0) return -1; //判断是否无解，然而这题其实不用

        x=mul(x,c/gcd,bg);
        ans+=x*M;//更新前k个方程组的答案
        M*=bg;//M为前k个m的lcm
        ans=(ans%M+M)%M;
    }
    return (ans%M+M)%M;
}
\end{lstlisting}

\subsection{dyhshai}

\begin{lstlisting}
ll T, n, pri[maxn], cur, mu[maxn], sum_mu[maxn];
bool vis[maxn];
map<ll, ll> mp_mu;
ll S_mu(ll x) {
    if (x < maxn) return sum_mu[x];
    if (mp_mu[x]) return mp_mu[x];
    ll ret = 1ll;
    for (ll i = 2, j; i <= x; i = j + 1) {
        j = x / (x / i);
        ret -= S_mu(x / i) * (j - i + 1);
    }
    return mp_mu[x] = ret;
}
ll S_phi(ll x) {
    ll ret = 0ll;
    for (ll i = 1, j; i <= x; i = j + 1) {
        j = x / (x / i);
        ret += (S_mu(j) - S_mu(i - 1)) * (x / i) * (x / i);
    }
    return ((ret - 1) >> 1) + 1;
}
int main() {
    scanf("%lld", &T);
    mu[1] = 1;
    for (int i = 2; i < maxn; i++) {
    if (!vis[i]) {
        pri[++cur] = i;
        mu[i] = -1;
    }
    for (int j = 1; j <= cur && i * pri[j] < maxn; j++) {
        vis[i * pri[j]] = true;
        if (i % pri[j])
            mu[i * pri[j]] = -mu[i];
        else {
            mu[i * pri[j]] = 0;
            break;
        }
    }
    }
    for (int i = 1; i < maxn; i++) sum_mu[i] = sum_mu[i - 1] + mu[i];
    while (T--) {
        scanf("%lld", &n);
        printf("%lld %lld\n", S_phi(n), S_mu(n));
    }
    return 0;
}
\end{lstlisting}

\subsection{FFT}

\begin{lstlisting}
const int N = 4000005;
const double PI = acos ( -1. );
int r[N];
struct Complex {
    double r, i;
    Complex ( double _r = 0.0, double _i = 0.0 ) : r ( _r ), i ( _i ) {};
    inline void real ( const double &x ) {
        r = x;
    }
    inline double real() {
        return r;
    }
    inline Complex operator+ ( const Complex x ) const {
        return Complex ( r + x.r, i + x.i );
    }
    inline Complex operator- ( const Complex x ) const {
        return Complex ( r - x.r, i - x.i );
    }
    inline Complex operator* ( const Complex x ) const {
        return Complex ( r * x.r - i * x.i, x.r * i + r * x.i );
    }
    inline void operator/= ( const double x ) {
        r /= x;
    }
};
Complex a[N], b[N];
void change ( Complex y[], int len ) {
    for ( int i = 0; i < len; i++ ) {
        if ( i < r[i] ) {
            swap ( y[i], y[r[i]] );
        }
    }
}
void fft ( Complex y[], int len, int on ) {
    change ( y, len );
    for ( int m = 1; m < len; m <<= 1 ) {
        Complex wn ( cos ( PI / m ), sin ( on * PI / m ) );
        for ( int r = m << 1, j = 0; j < len; j += r ) {
            Complex w ( 1.0, 0 );
            for ( int k = 0; k < m; k++, w = w * wn ) {
                Complex u = y[j + k];
                Complex t = w * y[j + k + m];
                y[j + k] = u + t;
                y[j + k + m] = u - t;
            }
        }
    }
    if ( on == -1 )
        for ( int i = 0; i < len; i++ ) {
            y[i] /= len;
        }
}
int main() {
    int n, m;
    while ( ~scanf ( "%d %d", &n, &m ) ) {
        int l = 0, len = 1;
        while ( len <= n+m ) {
            len <<= 1;
            l++;
        }
        for ( int i = 0; i < len; i++ ) {
            r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
        }
        int tmp;
        for ( int i = 0; i <= n; i++ ) {
            scanf ( "%d", &tmp );
            a[i] = Complex ( tmp * 1.0, 0 );
        }
        for ( int i = n + 1; i < len; i++ ) {
            a[i] = Complex ( 0, 0 );
        }
        for ( int i = 0; i <= m; i++ ) {
            scanf ( "%d", &tmp );
            b[i] = Complex ( tmp * 1.0, 0 );
        }
        for ( int i = m + 1; i < len; i++ ) {
            b[i] = Complex ( 0, 0 );
        }
        fft ( a, len, 1 );
        fft ( b, len, 1 );
        for ( int i = 0; i < len; i++ ) {
            a[i] = a[i] * b[i];
        }
        fft ( a, len, -1 );
        for ( int i = 0; i <= n + m; i++ ) {
            printf ( "%lld%c", ( long long ) ( a[i].r + 0.5 ), i == n + m ? '\n' : ' ' );
        }
    }
}
\end{lstlisting}

\subsection{NTT}

\begin{lstlisting}
const int N = 4000005;
const double PI = acos ( -1. );
const int mod=998244353, G = 3, Gi = 332748118;//gmin=3;
long long quickpow(long long a,long long n){
    if(a==0)return 0;
    long long ans=1;
    while(n){
        if(n&1)ans=ans*a%mod;
        n>>=1;
        a=a*a%mod;
    }
    return ans;
}
int r[N];
long long a[N],b[N];
long long inv;
void exgcd(long long a,long long b,long long &x,long long &y){
    if(b){
        exgcd(b,a%b,y,x);
        y-=x*(a/b);
    }else{
        x=1;
        y=0;
        return ;
    }
}
long long gao(long long a){
    long long x,y;
    exgcd(a,mod,x,y);
    if(x<0)x+=mod;
    return x;
}
inline void change ( long long y[], int len ) {
    for ( int i = 0; i < len; i++ ) {
        if ( i < r[i] ) {
            swap ( y[i], y[r[i]] );
        }
    }
}
void fft ( long long y[], int len, int on ) {
    change ( y, len );
    for ( int m = 1; m < len; m <<= 1 ) {
        long long wn =quickpow(on==1?G:Gi,(mod-1)/(m<<1));
        for ( int r = m << 1, j = 0; j < len; j += r ) {
            long long  w=1;
            for ( int k = 0; k < m; k++, w = w * wn%mod ) {
                long long u = y[j + k];
                long long t = w * y[j + k + m]%mod;
                y[j + k] =(u + t)%mod;
                y[j + k + m] = (u - t+mod)%mod;
            }
        }
    }
    if ( on == -1 )
        for ( int i = 0; i < len; i++ ) {
            y[i]=y[i]*inv%mod;
        }
}
int main() {
    int n, m;
    while ( ~scanf ( "%d %d", &n, &m ) ) {
        int l = 0, len = 1;
        while ( len <= n+m ) {
            len <<= 1;
            l++;
        }
        inv=gao(len);
        for ( int i = 0; i < len; i++ ) {
            r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
        }
        for ( int i = 0; i <= n; i++ ) {
            scanf ( "%lld", &a[i] );
        }
        for ( int i = n + 1; i < len; i++ ) {
            a[i] =0;
        }
        for ( int i = 0; i <= m; i++ ) {
            scanf ( "%lld", &b[i] );
        }
        for ( int i = m + 1; i < len; i++ ) {
            b[i] = 0;
        }
        fft ( a, len, 1 );
        fft ( b, len, 1 );
        for ( int i = 0; i < len; i++ ) {
            a[i] = a[i] * b[i]%mod;
        }
        fft ( a, len, -1 );
        for ( int i = 0; i <= n + m; i++ ) {
            printf ( "%lld%c", a[i], i == n + m ? '\n' : ' ' );
        }
    }
}
\end{lstlisting}

\subsection{FWT}
\begin{lstlisting}
long long inv;
long long a1[150000], b1[150000];
long long a2[150000], b2[150000];
long long a3[150000], b3[150000];
void FWT1 ( long long n ) {
    for ( int d = 1; d < n; d <<= 1 )
        for ( int m = d << 1, i = 0; i < n; i += m )
            for ( int j = 0; j < d; j++ ) {
                long long  x1 = a1[i + j], y1 = a1[i + j + d];
                long long  x2 = a2[i + j], y2 = a2[i + j + d];
                long long  x3 = a3[i + j], y3 = a3[i + j + d];
                a1[i + j + d] =  ( x1 + y1 ) % mod ;
                a2[i + j]     =  ( x2 + y2 ) % mod ;
                a3[i + j]     =  ( x3 + y3 ) % mod ;
                a3[i + j + d] =  ( x3 - y3 + mod ) % mod ;
                //xor:a[i+j]=x+y,a[i+j+d]=(x-y+mod)%mod;
                //and:a[i+j]=x+y;
                //or:a[i+j+d]=x+y;
            }
}
void FWT2 ( long long n ) {
    for ( int d = 1; d < n; d <<= 1 )
        for ( int m = d << 1, i = 0; i < n; i += m )
            for ( int j = 0; j < d; j++ ) {
                long long  x1 = b1[i + j], y1 = b1[i + j + d];
                long long  x2 = b2[i + j], y2 = b2[i + j + d];
                long long  x3 = b3[i + j], y3 = b3[i + j + d];
                b1[i + j + d] =  ( x1 + y1 ) % mod ;
                b2[i + j]     =  ( x2 + y2 ) % mod ;
                b3[i + j]     =  ( x3 + y3 ) % mod ;
                b3[i + j + d] =  ( x3 - y3 + mod ) % mod ;
                //xor:a[i+j]=x+y,a[i+j+d]=(x-y+mod)%mod;
                //and:a[i+j]=x+y;
                //or:a[i+j+d]=x+y;
            }
}
void UFWT ( long long n ) {
    for ( long long d = 1; d < n; d <<= 1 )
        for ( long long  m = d << 1, i = 0; i < n; i += m )
            for ( long long  j = 0; j < d; j++ ) {
                long long  x1 = a1[i + j], y1 = a1[i + j + d];
                long long  x2 = a2[i + j], y2 = a2[i + j + d];
                long long  x3 = a3[i + j], y3 = a3[i + j + d];
                a1[i + j + d] =  ( y1 - x1 + mod) % mod ;
                a2[i + j]     =  ( x2 - y2 + mod ) % mod ;
                a3[i + j]     =  ( x3 + y3 ) * inv % mod ;
                a3[i + j + d] =  ( x3 - y3 + mod ) * inv % mod ;
                //xor:a[i+j]=(x+y)/2,a[i+j+d]=(x-y)/2;
                //and:a[i+j]=x-y;
                //or:a[i+j+d]=y-x;
            }
}
void solve ( long long n ) {
    FWT1 ( n );
    FWT2 ( n );
    for ( long long i = 0; i < n; i++ ) {
        a1[i] = a1[i] * b1[i] % mod;
        a2[i] = a2[i] * b2[i] % mod;
        a3[i] = a3[i] * b3[i] % mod;
    }
    UFWT ( n );
}
int main() {
    long long n;
    while ( ~scanf ( "%lld", &n ) ) {
        long long res = 1 << n;
        for ( long long i = 0; i < res; i++ ) {
            scanf ( "%lld", &a1[i] );
            a3[i] = a2[i] = a1[i];
        }
        for ( long long i = 0; i < res; i++ ) {
            scanf ( "%lld", &b1[i] );
            b3[i] = b2[i] = b1[i];
        }
        inv = mod - ( mod >> 1 );
        solve ( res );
        for ( long long i = 0; i < res; i++ ) {
            printf ( "%lld%c", a1[i], i == res - 1 ? '\n' : ' ' );
        }
        for ( long long i = 0; i < res; i++ ) {
            printf ( "%lld%c", a2[i], i == res - 1 ? '\n' : ' ' );
        }
        for ( long long i = 0; i < res; i++ ) {
            printf ( "%lld%c", a3[i], i == res - 1 ? '\n' : ' ' );
        }
    }
    return 0;
}
\end{lstlisting}

\subsection{liner basis}

\begin{lstlisting}
struct basis {
    static const int MAXL=62;
    long long a[MAXL+1];
    basis() {
        reset();
    }
    void reset() {
        memset(a,0,sizeof a);
    }
    void insert(long long x) {
        for(int i=MAXL; i>=0; --i) {
            if(!(x>>i)&1)continue;

            if(a[i])x^=a[i];
            else {
                for(int j=0; j<i; j++)if((x>>j)&1)x^=a[j];
                for(int j=i+1; j<=MAXL; j++)if((a[j]>>i)&1)a[j]^=x;
                a[i]=x;
                return ;
            }
        }
    }
    long long qmax() {
        long long ans=0;
        for(int i=0; i<=MAXL; i++)ans^=a[i];
        return ans;
    }
};
inline void insert(long long x) {
  for (int i = 55; i + 1; i--) {
    if (!(x >> i))  // x的第i位是0
      continue;
    if (!p[i]) {
      p[i] = x;
      break;
    }
    x ^= p[i];
  }
}
\end{lstlisting}

\subsection{high precision}

\begin{lstlisting}
#define MAXN 9999
// MAXN 是一位中最大的数字
#define MAXSIZE 10024
// MAXSIZE 是位数
#define DLEN 4
// DLEN 记录压几位
struct Big {
    int a[MAXSIZE], len;
    bool flag;  //标记符号'-'
    Big() {
        len = 1;
        memset(a, 0, sizeof a);
        flag = 0;
    }
    Big(const int);
    Big(const char*);
    Big(const Big&);
    Big& operator=(const Big&);  //注意这里operator有&，因为赋值有修改……
    //由于OI中要求效率
    //此处不使用泛型函数
    //故不重载
    // istream& operator>>(istream&,  BigNum&);   //重载输入运算符
    // ostream& operator<<(ostream&,  BigNum&);   //重载输出运算符
    Big operator+(const Big&) const;
    Big operator-(const Big&) const;
    Big operator*(const Big&)const;
    Big operator/(const int&) const;
    // TODO: Big / Big;
    Big operator^(const int&) const;
    // TODO: Big ^ Big;

    // TODO: Big 位运算;

    int operator%(const int&) const;
    // TODO: Big ^ Big;
    bool operator<(const Big&) const;
    bool operator<(const int& t) const;
    inline void print();
};
// README::不要随随便便把参数都变成引用，那样没办法传值
Big::Big(const int b) {
    int c, d = b;
    len = 0;
    // memset(a,0,sizeof a);
    CLR(a);
    while (d > MAXN) {
        c = d - (d / (MAXN + 1) * (MAXN + 1));
        d = d / (MAXN + 1);
        a[len++] = c;
    }
    a[len++] = d;
}
Big::Big(const char* s) {
    int t, k, index, l;
    CLR(a);
    l = strlen(s);
    len = l / DLEN;
    if (l % DLEN) ++len;
    index = 0;
    for (int i = l - 1; i >= 0; i -= DLEN) {
        t = 0;
        k = i - DLEN + 1;
        if (k < 0) k = 0;
        g(j, k, i) t = t * 10 + s[j] - '0';
        a[index++] = t;
    }
}
Big::Big(const Big& T) : len(T.len) {
    CLR(a);
    f(i, 0, len) a[i] = T.a[i];
    // TODO:重载此处？
}
Big& Big::operator=(const Big& T) {
    CLR(a);
    len = T.len;
    f(i, 0, len) a[i] = T.a[i];
    return *this;
}
Big Big::operator+(const Big& T) const {
    Big t(*this);
    int big = len;
    if (T.len > len) big = T.len;
    f(i, 0, big) {
        t.a[i] += T.a[i];
        if (t.a[i] > MAXN) {
            ++t.a[i + 1];
            t.a[i] -= MAXN + 1;
        }
    }
    if (t.a[big])
    t.len = big + 1;
    else
    t.len = big;
    return t;
}
Big Big::operator-(const Big& T) const {
    int big;
    bool ctf;
    Big t1, t2;
    if (*this < T) {
    t1 = T;
    t2 = *this;
    ctf = 1;
    } else {
    t1 = *this;
    t2 = T;
    ctf = 0;
    }
    big = t1.len;
    int j = 0;
    f(i, 0, big) {
    if (t1.a[i] < t2.a[i]) {
        j = i + 1;
        while (t1.a[j] == 0) ++j;
        --t1.a[j--];
        // WTF?
        while (j > i) t1.a[j--] += MAXN;
        t1.a[i] += MAXN + 1 - t2.a[i];
    } else
        t1.a[i] -= t2.a[i];
    }
    t1.len = big;
    while (t1.len > 1 && t1.a[t1.len - 1] == 0) {
    --t1.len;
    --big;
    }
    if (ctf) t1.a[big - 1] = -t1.a[big - 1];
    return t1;
}
Big Big::operator*(const Big& T) const {
    Big res;
    int up;
    int te, tee;
    f(i, 0, len) {
        up = 0;
        f(j, 0, T.len) {
            te = a[i] * T.a[j] + res.a[i + j] + up;
            if (te > MAXN) {
                tee = te - te / (MAXN + 1) * (MAXN + 1);
                up = te / (MAXN + 1);
                res.a[i + j] = tee;
            } else {
                up = 0;
                res.a[i + j] = te;
            }
        }
        if (up) res.a[i + T.len] = up;
    }
    res.len = len + T.len;
    while (res.len > 1 && res.a[res.len - 1] == 0) --res.len;
    return res;
}
Big Big::operator/(const int& b) const {
    Big res;
    int down = 0;
    gd(i, len - 1, 0) {
        res.a[i] = (a[i] + down * (MAXN + 1) / b);
        down = a[i] + down * (MAXN + 1) - res.a[i] * b;
    }
    res.len = len;
    while (res.len > 1 && res.a[res.len - 1] == 0) --res.len;
    return res;
}
int Big::operator%(const int& b) const {
    int d = 0;
    gd(i, len - 1, 0) d = (d * (MAXN + 1) % b + a[i]) % b;
    return d;
}
Big Big::operator^(const int& n) const {
    Big t(n), res(1);
    // TODO::快速幂这样写好丑= =//DONE:)
    int y = n;
    while (y) {
        if (y & 1) res = res * t;
        t = t * t;
        y >>= 1;
    }
    return res;
}
bool Big::operator<(const Big& T) const {
    int ln;
    if (len < T.len) return 233;
    if (len == T.len) {
        ln = len - 1;
        while (ln >= 0 && a[ln] == T.a[ln]) --ln;
        if (ln >= 0 && a[ln] < T.a[ln]) return 233;
        return 0;
    }
    return 0;
}
inline bool Big::operator<(const int& t) const {
    Big tee(t);
    return *this < tee;
}
inline void Big::print() {
    printf("%d", a[len - 1]);
    gd(i, len - 2, 0) { printf("%04d", a[i]); }
}

inline void print(Big s) {
    // s不要是引用，要不然你怎么print(a * b);
    int len = s.len;
    printf("%d", s.a[len - 1]);
    gd(i, len - 2, 0) { printf("%04d", s.a[i]); }
}
char s[100024];
\end{lstlisting}

\subsection{Fibonacci}

\begin{lstlisting}
pair<int, int> fib(int n) {
    if (n == 0) return {0, 1};
    auto p = fib(n >> 1);
    int c = p.first * (2 * p.second - p.first);
    int d = p.first * p.first + p.second * p.second;
    if (n & 1)
        return {d, c + d};
    else
        return {c, d};
}      
\end{lstlisting}

\subsection{BM}

\begin{lstlisting}
#define rep(i,a,n) for (int i=a;i<n;i++)
#define per(i,a,n) for (int i=n-1;i>=a;i--)
#define pb push_back
#define mp make_pair
#define all(x) (x).begin(),(x).end()
#define fi first
#define se second
#define SZ(x) ((int)(x).size())
typedef vector<int> VI;
typedef long long ll;
typedef pair<int,int> PII;
const ll mod=1000000007;
ll powmod(ll a,ll b) {ll res=1;a%=mod; assert(b>=0); for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;}
// head

int _,n;
namespace linear_seq {
    const int N=10010;
    ll res[N],base[N],_c[N],_md[N];

    vector<int> Md;
    void mul(ll *a,ll *b,int k) {
        rep(i,0,k+k) _c[i]=0;
        rep(i,0,k) if (a[i]) rep(j,0,k) _c[i+j]=(_c[i+j]+a[i]*b[j])%mod;
        for (int i=k+k-1;i>=k;i--) if (_c[i])
            rep(j,0,SZ(Md)) _c[i-k+Md[j]]=(_c[i-k+Md[j]]-_c[i]*_md[Md[j]])%mod;
        rep(i,0,k) a[i]=_c[i];
    }
    int solve(ll n,VI a,VI b) { // a 系数 b 初值 b[n+1]=a[0]*b[n]+...
//        printf("%d\n",SZ(b));
        ll ans=0,pnt=0;
        int k=SZ(a);
        assert(SZ(a)==SZ(b));
        rep(i,0,k) _md[k-1-i]=-a[i];_md[k]=1;
        Md.clear();
        rep(i,0,k) if (_md[i]!=0) Md.push_back(i);
        rep(i,0,k) res[i]=base[i]=0;
        res[0]=1;
        while ((1ll<<pnt)<=n) pnt++;
        for (int p=pnt;p>=0;p--) {
            mul(res,res,k);
            if ((n>>p)&1) {
                for (int i=k-1;i>=0;i--) res[i+1]=res[i];res[0]=0;
                rep(j,0,SZ(Md)) res[Md[j]]=(res[Md[j]]-res[k]*_md[Md[j]])%mod;
            }
        }
        rep(i,0,k) ans=(ans+res[i]*b[i])%mod;
        if (ans<0) ans+=mod;
        return ans;
    }
    VI BM(VI s) {
        VI C(1,1),B(1,1);
        int L=0,m=1,b=1;
        rep(n,0,SZ(s)) {
            ll d=0;
            rep(i,0,L+1) d=(d+(ll)C[i]*s[n-i])%mod;
            if (d==0) ++m;
            else if (2*L<=n) {
                VI T=C;
                ll c=mod-d*powmod(b,mod-2)%mod;
                while (SZ(C)<SZ(B)+m) C.pb(0);
                rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod;
                L=n+1-L; B=T; b=d; m=1;
            } else {
                ll c=mod-d*powmod(b,mod-2)%mod;
                while (SZ(C)<SZ(B)+m) C.pb(0);
                rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod;
                ++m;
            }
        }
        return C;
    }
    int gao(VI a,ll n) {
        VI c=BM(a);
        c.erase(c.begin());
        rep(i,0,SZ(c)) c[i]=(mod-c[i])%mod;
        return solve(n,c,VI(a.begin(),a.begin()+SZ(c)));
    }
};
\end{lstlisting}

\subsection{quadratic residue}

\begin{lstlisting}
ll a, w;
struct field {
    ll x, y;
    field(ll a = 0, ll b = 0) {
        x = a; y = b;
    }
};
field operator*(field a, field b) { return field(a.x*b.x%mod + a.y*b.y%mod*w%mod, a.x*b.y%mod + a.y*b.x%mod); }

ll ran() {
    static ll seed = 23333;
    return seed = (((((seed ^ 20030927) % mod + 330802) % mod * 9410) % mod - 19750115 + mod) % mod ^ 730903) % mod;
}

ll pows(ll a, ll b) {
    ll base = 1;
    while (b) {
        if (b & 1) 
            base = base * a%mod;
        a = a * a%mod;
        b /= 2;
    }
    return base;
}

field powfield(field a, ll b) {
    field base = field(1, 0);
    while (b) {
        if (b & 1) base = base * a;
        a = a * a; b /= 2;
    }
    return base;
}

ll legander(ll x) {
    int a = pows(x, (mod - 1) / 2);
    if (a + 1 == mod) return -1;
    return a;
}

ll surplus(ll x) {
    ll a;
    if (x == 0) return 0;
    if (legander(x) == -1) return -1;
    while (1) {
        a = ran() % mod;
        w = ((a*a - x) % mod + mod) % mod;
        if (legander(w) == -1) break;
    }
    field b = field(a, 1);
    b = powfield(b, (mod + 1) / 2);
    return b.x;
}

\end{lstlisting}