bzoj 1409: Password 矩阵乘法+线性筛

题意

Rivest是密码学专家。近日他正在研究一种数列E = {E[1],E[2],……,E[n]}，
且E[1] = E[2] = p（p为一个质数），E[i] = E[i-2]*E[i-1] （若2< i<=n）。
例如{2,2,4,8,32,256,8192,……}就是p = 2的数列。在此基础上他又设计了一种加密算法，该算法可以通过一个密钥q (q < p)将一个正整数n加密成另外一个正整数d，计算公式为：d = E[n] mod q。现在Rivest想对一组数据进行加密，但他对程序设计不太感兴趣，请你帮助他设计一个数据加密程序。
0 < p n< 2^31 0 < q < p 0 < m <= 5000。

分析

注意到两项相乘其实就相当于指数相加。因为p,q互质所以可以用欧拉定理。那么只要把不大于sqrt(n)的素数都筛出来，然后求出phi(q)，就可以矩阵快速幂算指数了。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;

typedef long long LL;

const int N=50005;

int n,p,q,mo,prime[6005],tot;
bool not_prime[N];
struct Matrix
{
    LL a[3][3];

    void unit(int n)
    {
        memset(a,0,sizeof(a));
        for (int i=1;i<=n;i++) a[i][i]=1;
    }
}a;

int read()
{
    int x=0,f=1;char ch=getchar();
    while (ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while (ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}

void get_prime(int n)
{
    for (int i=2;i<=n;i++)
    {
        if (!not_prime[i]) prime[++tot]=i;
        for (int j=1;j<=tot&&i*prime[j]<=n;j++)
        {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0) break;
        }
    }
}

int ksm(int x,int y,int mo)
{
    int ans=1%mo;
    while (y)
    {
        if (y&1) ans=(LL)ans*x%mo;
        x=(LL)x*x%mo;y>>=1;
    }
    return ans;
}

void mul(Matrix &c,Matrix a,Matrix b)
{
    memset(c.a,0,sizeof(c.a));
    for (int i=1;i<=2;i++)
        for (int k=1;k<=2;k++)
            for (int j=1;j<=2;j++)
                (c.a[i][j]+=a.a[i][k]*b.a[k][j])%=mo;
}

Matrix Matrix_ksm(Matrix x,int y)
{
    Matrix ans;ans.unit(2);
    while (y)
    {
        if (y&1) mul(ans,ans,x);
        mul(x,x,x);y>>=1;
    }
    return ans;
}

int main()
{
    int T=read();p=read();
    get_prime(50000);
    while (T--)
    {
        n=read();q=read();mo=q;int tmp=q;
        if (n<=2)
        {
            printf("%d\n",p%q);
            continue;
        }
        for (int i=1;i<=tot&&prime[i]*prime[i]<=tmp;i++)
            if (tmp%prime[i]==0)
            {
                mo=mo/prime[i]*(prime[i]-1);
                while (tmp%prime[i]==0) tmp/=prime[i];
            }
        if (tmp>1) mo=mo/tmp*(tmp-1);
        a.a[1][1]=0;a.a[1][2]=a.a[2][1]=a.a[2][2]=1;
        a=Matrix_ksm(a,n-2);
        printf("%d\n",ksm(p,(a.a[1][2]+a.a[2][2])%mo,q));
    }
    return 0;
}