「套路学习」 O(1) 在线逆元

学习 / trick
Word count: 1.8k   |   Reading time≈ 9 min

不保证内容的正确性，本人代码交loj #161被卡了

Description

考虑这样一个问题，给定一个 $p$ ，$Q$ 次询问，每次询问一个数 $a$ 模 $p$ 的逆元，$Q \le 10^6, a <mod $。

由于 $a$ 的限制，我们无法直接预处理出所有的数的逆元，于是我们构造这么一个方程：

其中 $M$ 是我们设立的一个阈$(y\grave u)$值，假设我们可以快速找到这么一对满足条件的 $k, h$ 并且我们可以快速的预处理出 $M$ 以下的数的逆元，那么我们就可以解决这个问题。

我们把这个式子稍微改一下

这时，我们从数论书上抄了一个结论下来

tmp2

那么我们就可以从分母小于 $M$ 的分数中，找出一个 $\frac{q}{k}$ 使得 $|\frac{a}{p} - \frac{q}{k}| \le \frac{1}{k(M + 1)}\\$ ，然后再把这个式子乘回去

也就是 $h \le \frac{p}{M + 1}$ ，算上预处理的复杂度，这个部分的复杂度就是 $O(\frac{P}{M +1})$

现在第二部分就是如何找到一个与 $\frac{a}{p}$ 相邻的分数，我们暂且先考虑找到第一个大于的数，可以先暴力找出所有分母小于 $M$ 的所有分数排个序，然后把 $\frac ap$ 拿进去 lower_bound() 一下，复杂度就是 $O(M^2\log M^2 + q\log(M^2))$ ，考虑如何把 log 拿掉，我们发现分母小于 $M$ 的任意两个分数的差都大于 $\frac1{M^2}$ ，于是我们把每个分数乘上一个 $M^2$ 然后向下取整，每一个分数就会唯一对应一个 $[1, M^2]$ 里的整数。我们用一个 pair 存一下每一个值对应的分数，然后从后往前扫一遍这个数组，把空的位置设成这个位置后面的第一个值，就可以实现 $O(1)$ 查询了，小于的情况也是对称的，复杂度 $O(M^2+q)$

取 $M = p^{\frac 13}$ ，总复杂度 $O(p^{\frac 23} + q)$，实现了 O(1) 逆元，但事实上并没有比黑科技的 log 快，强还是松松强。

补充：还是别写这个垃圾了，太慢了

Code

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
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
// 日、不是说就慢一点么
#include <map>
#include <set>
#include <ctime>
#include <queue>
#include <stack>
#include <cmath>
#include <vector>
#include <bitset>
#include <cstdio>
#include <cctype>
#include <string>
#include <numeric>
#include <cstring>
#include <cassert>
#include <climits>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#include <functional>
using namespace std ;
namespace io {
    const int SIZE = (1 << 21) + 1;
    char ibuf[SIZE], *iS, *iT, obuf[SIZE], *oS = obuf, *oT = oS + SIZE - 1, c, qu[55]; int f, qr;
    // getchar
    #define gc() (iS == iT ? (iT = (iS = ibuf) + fread (ibuf, 1, SIZE, stdin), (iS == iT ? EOF : *iS ++)) : *iS ++)
    // print the remaining part
    inline void flush () {
        fwrite (obuf, 1, oS - obuf, stdout);
        oS = obuf;
    }
    // putchar
    inline void putc (char x) {
        *oS ++ = x;
        if (oS == oT) flush ();
    }
    // input a signed integer
    template <class I>
	inline void gi (I &x) {
        for (f = 1, c = gc(); c < '0' || c > '9'; c = gc()) if (c == '-') f = -1;
        for (x = 0; c <= '9' && c >= '0'; c = gc()) x = x * 10 + (c & 15); x *= f;
    }
    // print a signed integer
    template <class I>
	inline void print (I x) {
        if (!x) putc ('0'); if (x < 0) putc ('-'), x = -x;
        while (x) qu[++ qr] = x % 10 + '0',  x /= 10;
        while (qr) putc (qu[qr --]);
    }
    //no need to call flush at the end manually!
    struct Flusher_ {~Flusher_(){flush();}}io_flusher_;
}
using io :: gi;
using io :: putc;
using io :: print;

#define int long long
#define rep(i, a, b) for (ll i = (a); i <= (b); ++i)
#define per(i, a, b) for (ll i = (a); i >= (b); --i)
#define loop(it, v) for (auto it = v.begin(); it != v.end(); it++)
#define cont(i, x) for (ll i = head[x]; i; i = edge[i].nex)
#define clr(a) memset(a, 0, sizeof(a))
#define ass(a, cnt) memset(a, cnt, sizeof(a))
#define cop(a, b) memcpy(a, b, sizeof(a))
#define lowbit(x) (x & -x)
#define all(x) x.begin(), x.end()
#define SC(t, x) static_cast <t> (x)
#define ub upper_bound
#define lb lower_bound
#define pqueue priority_queue
#define mp make_pair
#define pb push_back
#define pof pop_front
#define pob pop_back
#define fi first
#define se second
#define y1 y1_
#define Pi acos(-1.0)
#define iv inline void
#define enter putchar('\n')
#define siz(x) ((ll)x.size())
#define file(x) freopen(x".in", "r", stdin),freopen(x".out", "w", stdout)
typedef double db ;
typedef long long ll ;
typedef unsigned long long ull ;
typedef pair <ll, ll> pii ;
typedef vector <ll> vi ;
typedef vector <pii> vii ;
typedef queue <ll> qi ;
typedef queue <pii> qii ;
typedef set <ll> si ;
typedef map <ll, ll> mii ;
typedef map <string, ll> msi ;
const ll maxn = 5e6 + 100 ;
const ll inf = 0x3f3f3f3f ;
const ll iinf = 1 << 30 ;
const ll linf = 2e18 ;
const ll mod = 1e9 + 7 ;
const double eps = 1e-9 ;
template <class T = ll> T chmin(T &a, T b) { return a = min(a, b);}
template <class T = ll> T chmax(T &a, T b) { return a = max(a, b);}
template <class T = ll> T read()
{
    T f = 1, a = 0;
    char ch = getchar() ;
    while (!isdigit(ch)) { if (ch == '-') f = -1 ; ch = getchar() ; }
    while (isdigit(ch)) { a =  (a << 3) + (a << 1) + ch - '0' ; ch = getchar() ; }
    return a * f ;
}

ll c = (ceil)(pow(mod, 1.0 / 3.0)) + 100, base = c * c * 2;

ll n;

ll inv[maxn], p[maxn], s[maxn];

pii pre[maxn], suf[maxn];

ll gcd(ll a, ll b)
{
    return !b ? a : gcd(b, a % b);
}

void init()
{
    inv[1] = 1;
    rep(i, 2, 5 * mod / c) inv[i] = (mod - mod / i) * inv[mod % i] % mod;
    pre[base] = suf[base] = mp(1, 1);
    pre[0] = suf[0] = mp(1, 0);
    rep(i, 2, c) rep(j, 1, i - 1)
    {
        ll now = base * j / i;
        if(pre[now].fi) continue;
        pre[now] = suf[now] = mp(i, j);
    }
    rep(i, 1, base) p[i] = (pre[i].fi != 0) ? i : p[i - 1];
    per(i, base, 1) s[i] = (suf[i].fi != 0) ? i : s[i + 1];
}

set <pair <db, pii>> nmsl;

ll get_inv(ll a)
{
    if(a <= mod / c) return inv[a];
    ll now = base * a / mod;
    // printf("now : %lld\n", now);
    db val = 1.0 * a / mod;
    nmsl.clear();
    ll tot = s[now];
    if(tot)
    {
        nmsl.insert(mp(fabs(1.0 * suf[tot].fi / suf[tot].se - val) * suf[tot].se, suf[tot]));
        tot = s[tot + 1];
        if(tot) nmsl.insert(mp(fabs(1.0 * suf[tot].fi / suf[tot].se - val) * suf[tot].se, suf[tot]));
    }
    tot = p[now - 1];
    nmsl.insert(mp(fabs(1.0 * suf[tot].fi / suf[tot].se - val) * suf[tot].se, suf[tot]));
    ll k = nmsl.begin() -> se.fi, h = nmsl.begin() -> se.se;
    // printf("%lld %lld %.10lf\n", k, h, 1.0 * h / k - val);
    // printf("%lld\n", a * k - h * mod);
    // cerr << siz(nmsl) << endl;
    // if(abs(a * k - h * mod) > 2 * mod / c) return puts("NMSL"), 0;
    ll tmp = (a * k >= h * mod) ? (inv[a * k - h * mod]) : (mod - inv[h * mod - a * k]);
    // cerr << "nmsl" << endl;
    return k * tmp % mod;
}

ll power(ll a, ll b)
{
    ll ret = 1;
    for(; b; b >>= 1, (a *= a) %= mod) if(b & 1) (ret *= a) %= mod;
    return ret;
}

ll a, ans;

ll fac[maxn];

signed main()
{
    // file("main");
    init();
    gi(n);
    fac[0] = 1;
    rep(i, 1, n) fac[i] = fac[i - 1] * 998244353 % mod;
    rep(i, 1, n)
    {
        // cerr << i << endl;
        // if(i % 10000 == 0) cerr << i << endl;
        gi(a);
        // cerr << a << endl;
        (ans += fac[n - i] * get_inv(a) % mod) %= mod;
        // printf("%lld\n", power(a, mod - 2));
    }
    print(ans);
    putc('\n');
    return 0;
}
  • Copyright： Copyright is owned by the author. For commercial reprints, please contact the author for authorization. For non-commercial reprints, please indicate the source.
  • Copyrights © 2020-2022 naitir
  • Visitors: | Views:

请我喝杯咖啡吧~

支付宝
微信