数学(论)里的一些定理（莫比乌斯反演，傅立叶变换，数论变换...）

莫比乌斯反演

莫比乌斯反演在数论中占有重要的地位，许多情况下能大大简化运算。那么我们先来认识莫比乌斯反演公式。

 

定理：和是定义在非负整数集合上的两个函数，并且满足条件，那么我们得到结论

 

     

 

在上面的公式中有一个函数，它的定义如下：

 

    （1）若，那么

    （2）若，均为互异素数，那么

    （3）其它情况下

 

 

对于函数，它有如下的常见性质：

 

    （1）对任意正整数有

  

                            

 

        （2）对任意正整数有

 

         

 

线性筛选求莫比乌斯反演函数代码。

void Init()
{
    memset(vis,0,sizeof(vis));
    mu[1] = 1;
    cnt = 0;
    for(int i=2; i<N; i++)
    {
        if(!vis[i])
        {
            prime[cnt++] = i;
            mu[i] = -1;
        }
        for(int j=0; j<cnt&&i*prime[j]<N; j++)
        {
            vis[i*prime[j]] = 1;
            if(i%prime[j]) mu[i*prime[j]] = -mu[i];
            else
            {
                mu[i*prime[j]] = 0;
                break;
            }
        }
    }
}

有了上面的知识，现在我们来证明莫比乌斯反演定理。

证明

证明完毕！

嗯，有了莫比乌斯反演，很多问题都可以简化了，接下来我们来看看莫比乌斯反演在数论中如何简化运算的。

题目：http://bz.cdqzoi.com/JudgeOnline/problem.php?id=2818

题意：给一个正整数，其中，求使得为质数的的个数，。

分析：对于本题，因为是使得为质数，所以必然要枚举小于等于的质数，那么对于每一个质数，只

     需要求在区间中，满足有序对互质的对数。

     也就是说，现在问题转化为：在区间中，存在多少个有序对使得互质，这个问题就简单啦，因为

     是有序对，不妨设，那么我们如果枚举每一个，小于有多少个与互素，这正是欧拉函数。所以

     我们可以递推法求欧拉函数，将得到的答案乘以2即可，但是这里乘以2后还有漏计算了的，那么有哪些呢？

     是且为素数的情况，再加上就行了。

代码：

#include <iostream>
#include <string.h>
#include <stdio.h>
#include <bitset>

using namespace std;
typedef long long LL;
const int N = 10000010;

bitset<N> prime;
LL phi[N];
LL f[N];
int p[N];
int k;

void isprime()
{
    k = 0;
    prime.set();
    for(int i=2; i<N; i++)
    {
        if(prime[i])
        {
            p[k++] = i;
            for(int j=i+i; j<N; j+=i)
                prime[j] = false;
        }
    }
}

void Init()
{
    for(int i=1; i<N; i++)  phi[i] = i;
    for(int i=2; i<N; i+=2) phi[i] >>= 1;
    for(int i=3; i<N; i+=2)
    {
        if(phi[i] == i)
        {
            for(int j=i; j<N; j+=i)
                phi[j] = phi[j] - phi[j] / i;
        }
    }
    f[1] = 0;
    for(int i=2;i<N;i++)
        f[i] = f[i-1] + (phi[i]<<1);
}

LL Solve(int n)
{
    LL ans = 0;
    for(int i=0; i<k&&p[i]<=n; i++)
        ans += 1 + f[n/p[i]];
    return ans;
}

int main()
{
    Init();
    isprime();
    int n;
    scanf("%d",&n);
    printf("%I64d
",Solve(n));
    return 0;
}

上题不算太难，普通的欧拉函数就可以搞定，接下来我们来看看它的升级版。

题意：给定两个数和，其中，，求为质数的有多少对？其中和的范

     围是。

 

分析：本题与上题不同的是和不一定相同。在这里我们用莫比乌斯反演来解决，文章开头也说了它能大大简化

     运算。我们知道莫比乌斯反演的一般描述为：

 

     

 

     其实它还有另一种描述，本题也是用到这种。那就是：

 

     

 

     好了，到了这里，我们开始进入正题。。。

 

     对于本题，我们设

 

     为满足且和的的对数

     为满足且和的的对数

 

     那么，很显然，反演后得到

 

     因为题目要求是为质数，那么我们枚举每一个质数，然后得到

 

     

 

     如果直接这样做肯定TLE，那么我们必须优化。

 

     我们设，那么继续得到。

 

     到了这里，可以看出如果我们可以先预处理出所有的对应的的值，那么本题就解决了。

 

     我们设，注意这里为素数，。

 

     那么，我们枚举每一个，得到，现在分情况讨论：

 

     （1）如果整除，那么得到

 

       

 

     （2）如果不整除，那么得到

 

       

代码：

#include <iostream>
#include <string.h>
#include <stdio.h>

using namespace std;
typedef long long LL;
const int N = 10000005;

bool vis[N];
int p[N];
int cnt;
int g[N],u[N],sum[N];

void Init()
{
    memset(vis,0,sizeof(vis));
    u[1] = 1;
    cnt = 0;
    for(int i=2;i<N;i++)
    {
        if(!vis[i])
        {
            p[cnt++] = i;
            u[i] = -1;
            g[i] = 1;
        }
        for(int j=0;j<cnt&&i*p[j]<N;j++)
        {
            vis[i*p[j]] = 1;
            if(i%p[j])
            {
                u[i*p[j]] = -u[i];
                g[i*p[j]] = u[i] - g[i];
            }
            else
            {
                u[i*p[j]] = 0;
                g[i*p[j]] = u[i];
                break;
            }
        }
    }
    sum[0] = 0;
    for(int i=1;i<N;i++)
        sum[i] = sum[i-1] + g[i];
}

int main()
{
    Init();
    int T;
    scanf("%d",&T);
    while(T--)
    {
        LL n,m;
        cin>>n>>m;
        if(n > m) swap(n,m);
        LL ans = 0;
        for(int i=1,last;i<=n;i=last+1)
        {
            last = min(n/(n/i),m/(m/i));
            ans += (n/i)*(m/i)*(sum[last]-sum[i-1]);
        }
        cout<<ans<<endl;
    }
    return 0;
}

多项式乘法运算初级版——快速傅里叶变换

快速傅里叶变换在信息学竞赛中主要用于求卷积，或者说多项式乘法。我们知道，多项式乘法的普通算法时间复杂度

是，通过快速傅里叶变换可以使时间降为，那么接下来会详细介绍快速傅里叶变换的原理。

 

首先来介绍多项式的两种表示方法，即系数表示法和点值表示法。从某种意义上说，这两种方法是等价的。先设

    

（1）系数表示法

 

    对于一个次数界为的多项式来说，其系数表示法就是一个由系数组成的向量，很

    明显，这样的多项式乘法运算的时间复杂度为。

 

（2）点值表示法

 

    对于一个次数界为的多项式来说，其点值是个点值对所形成的集合

 

    

 

    其中各不相同，并且当时，有。可以看出一个多项式可以有多种不同的点值

    表示法，而通过这个不同的点值对可以表示一个唯一的多项式。而通过点值表示法来计算多项式的乘法，时间

    复杂度为。

 

    从原则上来说，计算多项式的点值是简单易行的，因为我们只需要先选取个相异的点，然后通过秦九韶算法可

    以在时间内求出所有的，实际上如果我们的选得巧妙的话，就可以加速这一过程，使其运行时间变

    为。

 

    根据多项式的系数表示法求其点值表示法的过程称为求值，而根据点值表示法求其系数表示法的过程称为插值。

 

    对于求卷积或者说多项式乘法运算问题，先是通过傅里叶变换对系数表示法的多项式进行求值运算，这一步的时

    间复杂度为，然后在时间内进行点值相乘，再进行插值运算。

 

那么，接下来就是我们今天的重点了，如何高效地对一个多项式进行求值运算，即将多项式的表示法变为点值表示法。

 

如果选取单位复根作为求值点，则可以通过对系数向量进行离散傅里叶变换（DFT），得到相应的点值表示。同样地

也可以通过对点值对进行逆DFT运算，获得相应的系数向量。DFT和逆DFT的时间复杂度均为。

 

一. 求DFT

 

    选取次单位复根作为来求点值是比较巧妙的做法。

    次单位复根是满足的复数，次单位复根恰好有个，它们是，，为

    了解释这一式子，利用复数幂的定义，值称为主次单位根，所有其

    它次单位复根都是的次幂。

 

    个次单位复根在乘法运算下形成一个群，该群的结构与加法群模相同。

 

    接下来认识几个关于次单位复根的重要性质。

   

    （1）相消引理

 

        对于任何整数，有

 

    （2）折半引理

 

        如果且为偶数，则

 

    （3）求和引理

 

        对任意整数和不能被整除的非零整数，有

 

          

 

     回顾一下，我们希望计算次数界为的多项式

 

     

 

     在处的值，假定是2的幂，因为给定的次数界总可以增大，如果需要，总可以添加值为零

     的新的高阶系数。假定已知的系数形式为，对，定义结果

     如下

                 

     向量是系数向量的离散傅里叶变换，写作。

     通过使用一种称为快速傅里叶变换（FFT）的方法，就可以在时间内计算出，而直接

     计算的方法所需时间为，FFT主要是利用单位复根的特殊性质。FFT方法运用了分治策略，它用

     中偶数下标的系数与奇数下标的系数，分别定义了两个新的次数界为的多项式和

     

     则进一步有

 

     这样在处的值得问题就转换为求次数界为的多项式和在点

     处的值。由于在奇偶分类时导致顺序发生变化，所以需要先通过Rader算法进行

     倒位序，在FFT中最重要的一个操作是蝴蝶操作，通过蝴蝶操作可以将前半部分和后半部分的值求出。

 

题目：http://acm.hdu.edu.cn/showproblem.php?pid=1402

题意：大数乘法，需要用FFT实现。

代码：

#include <iostream>
#include <string.h>
#include <stdio.h>
#include <math.h>

using namespace std;
const int N = 500005;
const double PI = acos(-1.0);

struct Virt
{
    double r, i;

    Virt(double r = 0.0,double i = 0.0)
    {
        this->r = r;
        this->i = i;
    }

    Virt operator + (const Virt &x)
    {
        return Virt(r + x.r, i + x.i);
    }

    Virt operator - (const Virt &x)
    {
        return Virt(r - x.r, i - x.i);
    }

    Virt operator * (const Virt &x)
    {
        return Virt(r * x.r - i * x.i, i * x.r + r * x.i);
    }
};

//雷德算法--倒位序
void Rader(Virt F[], int len)
{
    int j = len >> 1;
    for(int i=1; i<len-1; i++)
    {
        if(i < j) swap(F[i], F[j]);
        int k = len >> 1;
        while(j >= k)
        {
            j -= k;
            k >>= 1;
        }
        if(j < k) j += k;
    }
}

//FFT实现
void FFT(Virt F[], int len, int on)
{
    Rader(F, len);
    for(int h=2; h<=len; h<<=1) //分治后计算长度为h的DFT
    {
        Virt wn(cos(-on*2*PI/h), sin(-on*2*PI/h));  //单位复根e^(2*PI/m)用欧拉公式展开
        for(int j=0; j<len; j+=h)
        {
            Virt w(1,0);            //旋转因子
            for(int k=j; k<j+h/2; k++)
            {
                Virt u = F[k];
                Virt t = w * F[k + h / 2];
                F[k] = u + t;     //蝴蝶合并操作
                F[k + h / 2] = u - t;
                w = w * wn;      //更新旋转因子
            }
        }
    }
    if(on == -1)
        for(int i=0; i<len; i++)
            F[i].r /= len;
}

//求卷积
void Conv(Virt a[],Virt b[],int len)
{
    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);
}

char str1[N],str2[N];
Virt va[N],vb[N];
int result[N];
int len;

void Init(char str1[],char str2[])
{
    int len1 = strlen(str1);
    int len2 = strlen(str2);
    len = 1;
    while(len < 2*len1 || len < 2*len2) len <<= 1;

    int i;
    for(i=0; i<len1; i++)
    {
        va[i].r = str1[len1-i-1] - '0';
        va[i].i = 0.0;
    }
    while(i < len)
    {
        va[i].r = va[i].i = 0.0;
        i++;
    }
    for(i=0; i<len2; i++)
    {
        vb[i].r = str2[len2-i-1] - '0';
        vb[i].i = 0.0;
    }
    while(i < len)
    {
        vb[i].r = vb[i].i = 0.0;
        i++;
    }
}

void Work()
{
    Conv(va,vb,len);
    for(int i=0; i<len; i++)
        result[i] = va[i].r+0.5;
}

void Export()
{
    for(int i=0; i<len; i++)
    {
        result[i+1] += result[i]/10;
        result[i] %= 10;
    }
    int high = 0;
    for(int i=len-1; i>=0; i--)
    {
        if(result[i])
        {
            high = i;
            break;
        }
    }
    for(int i=high; i>=0; i--)
        printf("%d",result[i]);
    puts("");
}

int main()
{
    while(~scanf("%s%s",str1,str2))
    {
        Init(str1,str2);
        Work();
        Export();
    }
    return 0;
}

题目：http://acm.hdu.edu.cn/showproblem.php?pid=4609

题意：给定n条长度已知的边，求能组成多少个三角形。

分析：用一个num数组来记录次数，比如num[i]表示长度为i的边有num[i]条。然后对num[]求卷积，除去本身重

     复的和对称的，然后再整理一下就好了。

代码：

#include <iostream>
#include <string.h>
#include <algorithm>
#include <stdio.h>
#include <math.h>

using namespace std;
typedef long long LL;

const int N = 400005;
const double PI = acos(-1.0);

struct Virt
{
    double r,i;

    Virt(double r = 0.0,double i = 0.0)
    {
        this->r = r;
        this->i = i;
    }

    Virt operator + (const Virt &x)
    {
        return Virt(r+x.r,i+x.i);
    }

    Virt operator - (const Virt &x)
    {
        return Virt(r-x.r,i-x.i);
    }

    Virt operator * (const Virt &x)
    {
        return Virt(r*x.r-i*x.i,i*x.r+r*x.i);
    }
};

//雷德算法--倒位序
void Rader(Virt F[],int len)
{
    int j = len >> 1;
    for(int i=1; i<len-1; i++)
    {
        if(i < j) swap(F[i], F[j]);
        int k = len >> 1;
        while(j >= k)
        {
            j -= k;
            k >>= 1;
        }
        if(j < k) j += k;
    }
}

//FFT实现
void FFT(Virt F[],int len,int on)
{
    Rader(F,len);
    for(int h=2; h<=len; h<<=1) //分治后计算长度为h的DFT
    {
        Virt wn(cos(-on*2*PI/h),sin(-on*2*PI/h));  //单位复根e^(2*PI/m)用欧拉公式展开
        for(int j=0; j<len; j+=h)
        {
            Virt w(1,0);            //旋转因子
            for(int k=j; k<j+h/2; k++)
            {
                Virt u = F[k];
                Virt t = w*F[k+h/2];
                F[k] = u+t;      //蝴蝶合并操作
                F[k+h/2] = u-t;
                w = w*wn;      //更新旋转因子
            }
        }
    }
    if(on == -1)
        for(int i=0; i<len; i++)
            F[i].r /= len;
}

//求卷积
void Conv(Virt F[],int len)
{
    FFT(F,len,1);
    for(int i=0; i<len; i++)
        F[i] = F[i]*F[i];
    FFT(F,len,-1);
}

int a[N];
Virt F[N];
LL num[N],sum[N];
int len,n;

void Init()
{
    memset(num,0,sizeof(num));
    scanf("%d",&n);
    for(int i=0; i<n; i++)
    {
        scanf("%d",&a[i]);
        num[a[i]]++;
    }
    sort(a, a + n);
    int len1 = a[n-1] + 1;
    len  = 1;
    while(len < len1*2) len <<= 1;
    for(int i=0; i<len1; i++)
        F[i] = Virt(num[i],0);
    for(int i=len1; i<len; i++)
        F[i] = Virt(0,0);
}

void Work()
{
    Conv(F,len);
    for(int i=0; i<len; i++)
        num[i] = (LL)(F[i].r+0.5);
    len = a[n-1]*2;
    for(int i=0; i<n; i++)
        num[a[i]+a[i]]--;
    for(int i=1; i<=len; i++)
        num[i] >>= 1;
    sum[0] = 0;
    for(int i=1; i<=len; i++)
        sum[i] = sum[i-1] + num[i];
    LL cnt = 0;
    for(int i=0; i<n; i++)
    {
        cnt+=sum[len]-sum[a[i]];
        //减掉一个取大，一个取小的
        cnt-=(LL)(n-1-i)*i;
        //减掉一个取本身，另外一个取其它
        cnt-=(n-1);
        //减掉大于它的取两个的组合
        cnt-=(LL)(n-1-i)*(n-i-2)/2;
    }
    LL tot = (LL)n*(n-1)*(n-2)/6;
    printf("%.7lf
",(double)cnt/tot);
}

int main()
{
    int T;
    scanf("%d",&T);
    while(T--)
    {
        Init();
        Work();
    }
    return 0;
}

多项式乘法运算终极版——NTT（快速数论变换）

在上一篇文章中 http://blog.csdn.net/acdreamers/article/details/39005227 介绍了用快速傅里叶变

换来求多项式的乘法。可以发现它是利用了单位复根的特殊性质，大大减少了运算，但是这种做法是对复数系数的矩阵

加以处理，每个复数系数的实部和虚部是一个正弦及余弦函数，因此大部分系数都是浮点数，我们必须做复数及浮点数

的计算，计算量会比较大，而且浮点数的计算可能会导致误差增大。

 

今天，我将来介绍另一种计算多项式乘法的算法，叫做快速数论变换（NTT），在离散正交变换的理论中，已经证明在

复数域内，具有循环卷积特性的唯一变换是DFT，所以在复数域中不存在具有循环卷积性质的更简单的离散正交变换。

因此提出了以数论为基础的具有循环卷积性质的快速数论变换。

 

回忆复数向量，其离散傅里叶变换公式如下

 

   

 

离散傅里叶逆变换公式为

 

   

 

今天的快速数论变换（NTT）是在上进行的，在快速傅里叶变换（FFT）中，通过次单位复根来运算的，即满

足的，而对于快速数论变换来说，则是可以将看成是的等价，这里是模素数

的原根（由于是素数，那么原根一定存在）。即

 

        

 

所以综上，我们得到数论变换的公式如下

 

    

 

而数论变换的逆变换公式为

 

    

 

这样就把复数对应到一个整数，之后一切都是在系统内考虑。

 

上述数论变换（NTT）公式中，要求是素数且必须是的因子。由于经常是2的方幂，所以可以构造形

如的素数。通常来说可以选择为费马素数，这样的变换叫做费马数数论变换。

 

这里我们选择，，这样得到模的原根值为。

 

题目：http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1028

分析：题目意思就是大数相乘，此处用快速数论变换（NTT）实现。

代码：

#include <iostream>
#include <string.h>
#include <stdio.h>

using namespace std;
typedef long long LL;

const int N = 1 << 18;
const int P = (479 << 21) + 1;
const int G = 3;
const int NUM = 20;

LL  wn[NUM];
LL  a[N], b[N];
char A[N], B[N];

LL quick_mod(LL a, LL b, LL m)
{
    LL ans = 1;
    a %= m;
    while(b)
    {
        if(b & 1)
        {
            ans = ans * a % m;
            b--;
        }
        b >>= 1;
        a = a * a % m;
    }
    return ans;
}

void GetWn()
{
    for(int i=0; i<NUM; i++)
    {
        int t = 1 << i;
        wn[i] = quick_mod(G, (P - 1) / t, P);
    }
}

void Prepare(char A[], char B[], LL a[], LL b[], int &len)
{
    len = 1;
    int len_A = strlen(A);
    int len_B = strlen(B);
    while(len <= 2 * len_A || len <= 2 * len_B) len <<= 1;
    for(int i=0; i<len_A; i++)
        A[len - 1 - i] = A[len_A - 1 - i];
    for(int i=0; i<len - len_A; i++)
        A[i] = '0';
    for(int i=0; i<len_B; i++)
        B[len - 1 - i] = B[len_B - 1 - i];
    for(int i=0; i<len - len_B; i++)
        B[i] = '0';
    for(int i=0; i<len; i++)
        a[len - 1 - i] = A[i] - '0';
    for(int i=0; i<len; i++)
        b[len - 1 - i] = B[i] - '0';
}

void Rader(LL a[], int len)
{
    int j = len >> 1;
    for(int i=1; i<len-1; i++)
    {
        if(i < j) swap(a[i], a[j]);
        int k = len >> 1;
        while(j >= k)
        {
            j -= k;
            k >>= 1;
        }
        if(j < k) j += k;
    }
}

void NTT(LL a[], int len, int on)
{
    Rader(a, len);
    int id = 0;
    for(int h = 2; h <= len; h <<= 1)
    {
        id++;
        for(int j = 0; j < len; j += h)
        {
            LL w = 1;
            for(int k = j; k < j + h / 2; k++)
            {
                LL u = a[k] % P;
                LL t = w * (a[k + h / 2] % P) % P;
                a[k] = (u + t) % P;
                a[k + h / 2] = ((u - t) % P + P) % P;
                w = w * wn[id] % P;
            }
        }
    }
    if(on == -1)
    {
        for(int i = 1; i < len / 2; i++)
            swap(a[i], a[len - i]);
        LL Inv = quick_mod(len, P - 2, P);
        for(int i = 0; i < len; i++)
            a[i] = a[i] % P * Inv % P;
    }
}

void Conv(LL a[], LL b[], int n)
{
    NTT(a, n, 1);
    NTT(b, n, 1);
    for(int i = 0; i < n; i++)
        a[i] = a[i] * b[i] % P;
    NTT(a, n, -1);
}

void Transfer(LL a[], int n)
{
    int t = 0;
    for(int i = 0; i < n; i++)
    {
        a[i] += t;
        if(a[i] > 9)
        {
            t = a[i] / 10;
            a[i] %= 10;
        }
        else t = 0;
    }
}

void Print(LL a[], int n)
{
    bool flag = 1;
    for(int i = n - 1; i >= 0; i--)
    {
        if(a[i] != 0 && flag)
        {
            printf("%d", a[i]);
            flag = 0;
        }
        else if(!flag)
            printf("%d", a[i]);
    }
    puts("");
}

int main()
{
    GetWn();
    while(scanf("%s%s", A, B)!=EOF)
    {
        int len;
        Prepare(A, B, a, b, len);
        Conv(a, b, len);
        Transfer(a, len);
        Print(a, len);
    }
    return 0;
}