逆元（除法取模）
我们需要解决的问题是求出 (b/a) mod p的值。
首先我们需要知道逆元这一玩意。
逆元的定义是这样的：
若a*c mod p=1（在以下的讨论中，由于写不出恒等式，所以均视为模p意义下的运算），则c
是a在mod p意义下的逆元。
这样在求(b/a) mod p的时候，就相当于求b*c mod p了。
下面有一些计算逆元的方法：
（记a的逆元为c，不考虑a=1的情况）
(a,p)=1
扩展欧几里得算法
若(a,p)=1，则(a*x-p*y) mod p=1的一个解x (mod p)就是a的逆元c。
显然，后面的p*y mod p是可以删去的，所以剩下的就是a*x mod p=1，所以x mod p就是a的逆元。
我们用扩展欧几里得就可以算出一个特解x。
不过有趣的是，用扩展欧几里得的时候可能是会得出负解x，所以这个时候就要把x变为正的。这个可以通过加p再模p的方法。
代码如下：
var
a,b,p,w,inverb:int64;//计算(a/b) mod p，inverb为b的逆元
x,y:int64;//扩展欧几里得的一组解
function multiply(a,b,p:int64):int64;
begin
multiply:=0;
while b<>0 do
begin
if b and 1<>0 then multiply:=(multiply+a) mod p;
a:=a shl 1 mod p;
b:=b shr 1;
end;
end;
function Gcd(a,b:int64):int64;//扩展欧几里得求一组特解
var
t:int64;
begin
if b=0 then
begin
x:=1;
y:=0;
exit(a);
end;
Gcd:=Gcd(b,a mod b);
t:=y;
y:=x-(a div b)*y;
x:=t;
end;
begin
read(a,b,p);//读入
Gcd(b,p);//计算x
x:=(x+(abs(x) div p+1)*p) mod p;//将特解x转为正数
inverb:=x mod p;//计算逆元
w:=multiply(a,inverb,p);//计算答案
writeln(w);//输出
end.
欧拉定理
若(a,p)=1，则a的逆元c mod p=a^(phi(p)-1) mod p。
这个也是很简单的：
a*a^(phi(p)-1) mod p=a^phi(p) mod p=1
代码如下：
type
inf=record
x:int64;
next:longint;
end;
var
a,b,p,inverb,w:int64;
top:longint;
Hash:array[0..100047] of longint;
Fin:array[1..100000] of inf;//小范围求欧拉函数时用通式求解（质因数分解用的是pollard_rho）
function multiply(a,b,p:int64):int64;
begin
multiply:=0;
while b<>0 do
begin
if b and 1<>0 then multiply:=(multiply+a) mod p;
a:=a shl 1 mod p;
b:=b shr 1;
end;
end;
function power(a,b,p:int64):int64;
var
s,w:int64;
begin
s:=1;
w:=a;
while b<>0 do
begin
if b and 1<>0 then s:=multiply(s,w,p);
w:=multiply(w,w,p);
b:=b shr 1;
end;
exit(s);
end;
function Miller(x:int64):boolean;
var
vi,vk:longint;
vj,vl,u,t:int64;
begin
if x=1 then exit(false);
if x=2 then exit(true);
u:=x-1;
t:=0;
while u and 1=0 do
begin
u:=u shr 1;
inc(t);
end;
for vi:=1 to 8 do
begin
vj:=2+random(x-2);
vj:=power(vj,u,x);
for vk:=1 to t do
begin
vl:=multiply(vj,vj,x);
if vl=1 then
if (vj<>1) and (vj<>x-1) then exit(false);
vj:=vl;
end;
if vj<>1 then exit(false);
end;
exit(true);
end;
function H(x:int64):longint;
begin
exit((x shl 1) mod 100047);
end;
procedure ins(x:int64);
var
v,u:longint;
begin
v:=H(x);
u:=Hash[v];
while u<>0 do
begin
if Fin[u].x=x then exit;
u:=Fin[u].next;
end;
inc(top);
Fin[top].x:=x;
Fin[top].next:=Hash[v];
Hash[v]:=top;
end;
function gcd(a,b:int64):int64;
var
c:int64;
begin
c:=a mod b;
while c<>0 do
begin
a:=b;
b:=c;
c:=a mod b;
end;
exit(b);
end;
function pollard(n,c:int64):int64;
var
x,y,d,s,t:int64;
begin
s:=1;
t:=2;
x:=1+random(n-2);
y:=x;
while true do
begin
x:=(multiply(x,x,n)+c) mod n;
d:=gcd(abs(x-y),n);
if (d<>1) and (d<>n) then exit(d);
if x=y then exit(n);
inc(s);
if s=t then
begin
t:=t shl 1;
y:=x;
end;
end;
end;
procedure decompose(x:int64);
var
y:int64;
begin
if Miller(x) then
begin
ins(x);
exit;
end;
y:=x;
while y>=x do
y:=pollard(x,2+random(x-2));
decompose(y);
decompose(x div y);
end;
function phi(x:int64):int64;//求欧拉函数
var
vi:longint;
begin
decompose(x);
for vi:=1 to top do
x:=x div Fin[vi].x*(Fin[vi].x-1);
exit(x);
end;
begin
read(a,b,p);//读入
inverb:=power(b,phi(p)-1,p);//计算逆元
w:=multiply(a,inverb,p);//计算答案
writeln(w);//输出
end.
p为质数
费马小定理
若p为质数，则a的逆元c mod p=a^(p-2) mod p
显然这也是很简单的：
由于a^(p-1) mod p=1，所以a*a^(p-2) mod p=1。
代码如下：
var
a,b,p,inverb,w:int64;
function multiply(a,b,p:int64):int64;
begin
multiply:=0;
while b<>0 do
begin
if b and 1<>0 then multiply:=(multiply+a) mod p;
a:=a shl 1 mod p;
b:=b shr 1;
end;
end;
function power(a,b,p:int64):int64;//写个快速幂就可以了，这还是很不错的选择。
var
s,w:int64;
begin
s:=1;
w:=a;
while b<>0 do
begin
if b and 1<>0 then s:=multiply(s,w,p);
w:=multiply(w,w,p);
b:=b shr 1;
end;
exit(s);
end;
begin
read(a,b,p);//读入
inverb:=power(b,p-2,p);//计算逆元
w:=multiply(a,inverb,p);//计算答案
writeln(w);//输出
end.
线性筛
首先是一个性质——逆元是一个完全积性函数。
这是什么意思呢？也就是说，令f(n)表示n的逆元，那么会有f(a*b)=f(a)*f(b)。
现在思路就比较简单了。一个朴素的思想就是先质因数分解x，之后求出x所有的质因子的逆元，相乘即可得到答案。
对于小数据范围（单独判断的问题）是可以的。如果我们要大范围地求出n以内的逆元，就可能需要线性筛了。
记得我们在求欧拉函数的时候是怎么做的吗？
欧拉函数是一个及积性函数（但不是完全积性函数）。也就是说，只有(a,b)=1时，才会有phi(a)*phi(b)=phi(a*b)。现在逆元是一个完全积性函数，自然也满足积性函数的性质。
所以参考线性筛求欧拉函数的做法来计算大规模的逆元。
详情转链接
任意条件
递推式
事实上f(n)（逆元的函数）是有一个递推式的。
令p=n*t+k，那么有f(n) mod p=n*(t*f(k))^2 mod p
嗯……这个证明就是抄别人的了。
p mod p=0=(n*t+k) mod p
即
n*t mod p=-k mod p
n*t*f(k) mod p=-1 mod p
(n*t*f(k)) mod p=1
f(n)=n*(t*f(k))^2 mod p
在模意义下n^(-1)就是f(n)（逆元的表示方法）
所以用一个递推式也是可以解出逆元的。
常处理的条件
求[1..n]的逆元
这个我们可以通过递推得到。类似于上面的部分。
假设我们已经得到的inv[1..i-1]，求inv[i]。
那么设intp=i*a+b，即i*a+b mod intp=0。两边同时除以i*b，就可以得到a/b+1/i mod intp=0。
也就是说，inv[i]=(intp div i)*(-inv[intp mod i])。注意转换为正数。
求1!,2!,...,n!的逆元
这个也是可以递推得到的。一般处理质数。
先求出阶乘fac[i]=fac[i-1]*i。
求出最后一项的逆元invfac[i]。
在逆向退回去，invfac[i]=invfac[i+1]*(i+1)。