P3307 [SDOI2013]项链

发表于 | 阅读数

  • 珠子=正三菱柱,三个面上面的数字$x$,必须满足$x\in [1,a]$，且珠子上面的数字的最大公约数要恰好为1。三棱柱可以选择翻转.

  • 相邻的两个珠子必须不同。

  • 两串项链如果能够经过旋转变成一样的，那么这两串项链就是相同的

$n<=10^{14}，a<=10^7，T<=10$

好题！

首先算珠子种类，$s3$表示$3$个$\gcd=1$的个数，同理$s2,s1$。

$Burnside引理$

然后之后$Polya$里的$m=res$

这里矩阵快速幂就太大了。

  • 如果$i-1$与$1$不同则$f(i)=f(i-1)(m-2)$

  • 如果$i-1$与$1$相同则$i-2$与$1$不同。 $f(i)=f(i-2)(m-1)$

$x^2-(m-1)x-(m-2)=(x+1)(x-(m-1))=0$

$f(i)=a(-1)^i+b(m-1)^i,f(1)=0,f(2)=m(m-1)$

$f(1)$由于$1$与$1$相邻。所以为$0$。

然后$dfs$一遍搜因数并且一遍处理欧拉函数。

BZOJ上的数据经过加强， n 可能是 $10^9+7$ 的倍数（没有逆元），而最后要除以 $n$ ，所以需要一些特殊的技巧。

  • $MOD=(1e9+7)^2$
  • 此时要用快速乘
  • 答案等于在模 $10^9+7$ 意义下 $\frac{sum}{1e9+7}$ 除以 $\frac{n}{1e9+7}$ 的结果。
  • 需要$exgcd$求$6$在$(1e9+7)^2$的逆元。
代码 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153

#include <bits/stdc++.h>
using namespace std;
typedef unsigned long long ll;
#define pii pair<ll, ll>
#define mk make_pair
const ll N = 1e7 + 10;
const ll MOD = 1000000007ll;
const ll Inv61 = 166666668ll;
const ll Inv62 = 833333345000000041ll;

ll mod;
ll Inv6;
ll O;
inline ll mult(ll x, ll y)
{
    if (!O)
        return x * y % mod;
    return (x * y - (ll)((long double)x / mod * y) * mod + mod) % mod;
}

int pcnt, pr[N], npr[N];
int mu[N], usum[N];

void Prime_init(ll n)
{
    npr[1] = 1;
    mu[1] = 1;
    for (ll i = 1; i <= n; i++)
    {
        if (!npr[i])
            pr[++pcnt] = i, mu[i] = -1;
        for (ll j = 1; j <= pcnt && pr[j] * i <= n; j++)
        {
            npr[pr[j] * i] = 1;

            if (i % pr[j] == 0)
            {
                mu[pr[j] * i] = 0;
                break;
            }
            else
            {
                mu[pr[j] * i] = mu[i] * (-1);
            }
        }
    }
    for (ll i = 1; i <= n; i++)
    {
        usum[i] = usum[i - 1] + mu[i];
    }
}

ll k;
ll p[100], e[100];
ll cnt;
ll ans;
ll qpow(ll a, ll x)
{
    ll res = 1;
    while (x)
    {
        if (x & 1)
            res = mult(res, a);
        a = mult(a, a);
        x >>= 1;
    }
    return res;
}
ll calc(ll x)
{
    return (qpow(k - 1, x) + mult((x & 1 ? mod - 1 : 1), (k - 1))) % mod;
}

void dfs(ll x, ll now, ll phi, ll n)
{
    if (x > cnt)
    {
        ans = (ans + mult(phi % mod, calc(n / now))) % mod;
        //cout << phi << " " << calc(n / now) << endl;
        return;
    }
    dfs(x + 1, now, phi, n);
    for (ll i = 1; i <= e[x]; i++)
    {
        phi *= (p[x] - (i == 1));
        now *= p[x];
        dfs(x + 1, now, phi, n);
    }
}
void solve(ll n)
{
    cnt = 0;
    ll tmp = n;
    for (ll i = 2; i * i <= n; i++)
    {
        if (n % i == 0)
        {
            ll o = 0;
            while (n % i == 0)
            {
                n /= i;
                o++;
            }
            p[++cnt] = i;
            e[cnt] = o;
        }
    }
    if (n > 1)
    {
        p[++cnt] = n;
        e[cnt] = 1;
    }
    dfs(1, 1, 1, tmp);
}
void init(ll n)
{
    k = 2;
    for (ll i = 1, j; i <= n; i = j + 1)
    {
        j = n / (n / i);
        k = (k + mult(mult(mult(n / i, n / i), n / i + 3), (usum[j] - usum[i - 1] + mod) % mod)) % mod;
    }
    k = mult(Inv6, k);
}

int main()
{
    ll T;
    cin >> T;
    Prime_init(1e7);
    while (T--)
    {
        ll n, a;
        cin >> n >> a;
        ans = 0;

        O = n % MOD ? 0 : 1;
        if (O)
            mod = MOD * MOD, Inv6 = Inv62;
        else
            mod = MOD, Inv6 = Inv61;

        init(a);
        solve(n);

        if (O)
            ans = mult(ans / MOD, qpow(n / MOD, MOD - 2)) % MOD;
        else
            ans = mult(ans, qpow(n % MOD, MOD - 2));
        cout << ans << endl;
    }
}