【递推+矩阵快速幂+模运算+欧拉降幂】HDU - 5895 Mathematician QSC

HDU - 5895

解题思路

太难了搞了我整整一天,数论咋这么难

参考题解

1.累加法化简 g ( n ) g(n) g(n)

img

在看看g(n)的定义 g(n) = f(1)^2 + f(2)^2 + … + f(n)^2;

因为

img

所以:

img

多写几个就可以看出累加法:

img

所以:

img

2.矩阵快速幂求 f ( n ) f(n) f(n)


[ f ( 1 ) , f ( 0 ) ] × A n = [ f ( n + 1 ) , f ( n ) ] [f(1),f(0)] \times A^{n} = [f(n+1),f(n)] [f(1),f(0)]×An=[f(n+1),f(n)]
其中 f ( 1 ) = 1 , f ( 0 ) = 0 , f(1)=1, f(0)=0, f(1)=1,f(0)=0,
A = [ 2 1 1 0 ] A= \left[ \begin{matrix} 2 & 1\\ 1 & 0 \end{matrix} \right] A=[2110]
设结果矩阵为 r e s res res,则 f ( n + 1 ) = r e s [ 0 ] [ 0 ] , f ( n ) = r e s [ 0 ] [ 1 ] f(n+1)=res[0][0],f(n)=res[0][1] f(n+1)=res[0][0],f(n)=res[0][1]

3.模运算的处理。

现在的式子化为: x g ( n ) % ( s + 1 ) = x f ( n ) × f ( n + 1 ) 2 % ( s + 1 ) x^{g(n)}\%(s+1)=x^{\frac {f(n) \times f(n+1)}{2}}\%(s+1) xg(n)%(s+1)=x2f(n)×f(n+1)%(s+1)

涉及幂运算的取模,运用欧拉降幂公式的第三种情况
a b % p = { a b % ϕ ( p ) % p g c d ( a , p ) = 1 a b % p g c d ( a , p ) ≠ 1 , b < ϕ ( p ) a b % ϕ ( p ) + ϕ ( p ) % p g c d ( a , p ) ≠ 1 , b ≥ ϕ ( p ) a^b\%p=\left\{ \begin{aligned} &a^{b\%\phi(p)}\%p& \quad gcd(a,p)=1\\ &a^b\%p& \quad gcd(a,p) \neq1,b<\phi(p)\\ &a^{b\%\phi(p)+\phi(p)}\%p& \quad gcd(a,p) \neq1,b\ge\phi(p)\\ \end{aligned} \right . ab%p=ab%ϕ(p)%pab%pab%ϕ(p)+ϕ(p)%pgcd(a,p)=1gcd(a,p)=1,b<ϕ(p)gcd(a,p)=1,bϕ(p)
化为: x f ( n ) × f ( n + 1 ) 2 % ϕ ( s + 1 ) + ϕ ( s + 1 ) % ( s + 1 ) x^{\frac {f(n) \times f(n+1)}{2}\%\phi(s+1)+\phi(s+1)}\%(s+1) x2f(n)×f(n+1)%ϕ(s+1)+ϕ(s+1)%(s+1)

单独求一个欧拉函数的板子:

LL phi(LL n) {
    LL res = n;
    for(int i = 2; i * i <= n; i++) {
        if (n % i == 0) {
            res = res - res / i;
            while (n % i == 0) n /= i;
        }
    }
    if (n > 1) res = res / n * (n - 1);
    return res;
}

不知道啥是欧拉函数戳这里,很详细~

单独拎出 x x x 的指数部分的前半部分看,要求 f ( n ) × f ( n + 1 ) 2 % ϕ ( s + 1 ) \frac {f(n) \times f(n+1)}{2}\%\phi(s+1) 2f(n)×f(n+1)%ϕ(s+1)

因为 2 2 2 ϕ ( s + 1 ) \phi(s+1) ϕ(s+1) 不一定互质,因此不能用乘法逆元求解。

正确做法是应用这个公式: a b % m o d = a % m o d × b b \frac ab\%mod=\frac {a\%mod\times b}{b} ba%mod=ba%mod×b

  • 证明:
    a b m o d    k = d a b = k x + d a = k b x + b d a m o d    k b = b d a m o d    k b b = d \frac ab\mod k = d \\ \frac ab=kx+d\\ a=kbx+bd\\ a\mod kb=bd\\ \frac {a \mod kb}{b}=d bamodk=dba=kx+da=kbx+bdamodkb=bdbamodkb=d

因此 f ( n ) × f ( n + 1 ) 2 % ϕ ( s + 1 ) = f ( n ) × f ( n + 1 ) % 2 × ϕ ( s + 1 ) 2 \frac {f(n) \times f(n+1)}{2}\%\phi(s+1)= \frac {f(n) \times f(n+1) \%2 \times\phi(s+1)}{2} 2f(n)×f(n+1)%ϕ(s+1)=2f(n)×f(n+1)%2×ϕ(s+1),可劲儿求吧。

注意事项

f ( n ) f(n) f(n) 的时候运用矩阵快速幂取模,注意此时的模不是 s + 1 s+1 s+1,而是 2 × ϕ ( s + 1 ) 2 \times\phi(s+1) 2×ϕ(s+1)

参考代码

#include<stdio.h>
#include<iostream>
#include<vector>
#include<cstring>
#include<cstdio>
#include<climits>
#include<cmath>
#include<algorithm>
#include<queue>
#include<deque>
#include<map>
#include<set>
#include<stack>
//#define LOCAL  //提交时一定注释
#define VI vector<int>
#define io ios::sync_with_stdio(false); cin.tie(0); cout.tie(0)
using namespace std;
typedef long long LL;
typedef double db;
const int inf = 0x3f3f3f3f;
const LL INF = 1e18;
typedef vector<VI> mat;

inline int readint() {int x; scanf("%d", &x); return x;}

LL mod;

mat multi(mat &A, mat &B) {
    mat C(A.size(), VI(B[0].size()));
    for(int i = 0; i < A.size(); i++) {
        for(int j = 0; j < B.size(); j++) {
            for(int k = 0; k < B.size(); k++) {
                C[i][j] = (C[i][j] + 1LL * A[i][k] * B[k][j] % mod) % mod;
            }
        }
    }
    return C;
}

mat fm(mat A, LL n) {
    mat res(A.size(), VI(A.size()));
    for(int i = 0; i < A.size(); i++)
        res[i][i] = 1;
    while (n) {
        if (n & 1) res = multi(res, A);
        A = multi(A, A);
        n >>= 1;
    }
    return res;
}

LL phi(LL n) {
    LL res = n;
    for(int i = 2; i * i <= n; i++) {
        if (n % i == 0) {
            res = res - res / i;
            while (n % i == 0) n /= i;
        }
    }
    if (n > 1) res = res / n * (n - 1);
    return res;
}

LL qm(LL a, LL b, LL c) {
    LL res = 1;
    while (b) {
        if (b & 1) res = res * a % c;
        a = a * a % c;
        b >>= 1;
    }
    return res;
}

int main() {
#ifdef LOCAL
    freopen("input.txt", "r", stdin);
//   freopen("output.txt", "w", stdout);
#endif
    int t = readint();
    mat A(2, VI(2)), B(1, VI(2));
    A[0][0] = 2;
    A[0][1] = A[1][0] = 1;
    B[0][0] = 1, B[0][1] = 0;
    LL n, y, x, s;
    mat res, A_n;
    while (t--) {
        scanf("%lld%lld%lld%lld", &n, &y, &x, &s);
        mod = 2 * phi(s + 1);
        A_n = fm(A, n * y);  //求A^n
        res = multi(B, A_n);
        LL mi = (res[0][0] % mod * res[0][1] % mod) % mod;
        mi = mi / 2 + mod / 2;
        printf("%lld\n", qm(x, mi, s + 1));
    }
    return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值