2015 ICL, Finals, Div. 1 Ceizenpok’s formula && Gym - 100633J

最新推荐文章于 2019-09-04 23:36:18 发布

原创最新推荐文章于 2019-09-04 23:36:18 发布 · 382 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#Gym

数论专栏收录该内容

7 篇文章

订阅专栏

上来先感谢大佬的文章：https://blog.youkuaiyun.com/clove_unique/article/details/54571216

题目链接：http://codeforces.com/gym/100633/problem/J

题意：

求
$C k n % p$ $C_n^k\%p$ 的值，其中 $1\le n \le 10^{18}$ , $0\le k \le n$ , $2 \le p \le 1000000$ .

分析：

若p为质数，可以用lucas定理直接得，但题目中p不一定为质数，则需要用拓展lucas定理+中国剩余定理。

结论：

拓展lucas定理：

可求方程：

a n s \equiv c (m o d p k 1 1)

$ans\equiv c(mod\, p_1^{k_1})$

中国剩余定理（又名孙子定理）

可求

⎧ ⎩ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ a n s \equiv c 1 (m o d m 1) a n s \equiv c 2 (m o d m 2) . . . . a n s \equiv c q (m o d m q)

$\begin{cases} ans\equiv c_1(mod \, m_1)\\ ans\equiv c_2(mod\, m_2)\\ ....\\ ans\equiv c_q(mod\, m_q) \end{cases}$
的最小正整数解，解为

(\sum i = 1 q * c i * M m i * i n v (M m i, m i)) % M

$(\sum_{i=1}^q*c_i*\frac{M}{m_i}*inv(\frac{M}{m_i},m_i))\%M$ 其中

M = \prod 1 q m i

$M=\prod_{1}^{q} m_i$

题解：

对于p，将它表示成

p = p k 1 1 * p k 2 2 . . . p k q q

$p=p_1^{k_1}*p_2^{k_2}...p_q^{k_q}$ 的形式，对于每个

pi p i $p_i$ ，用lucas定理求解，再用中国剩余定理求和。

唯一的问题在于在求

n ! ( n - m ) ! m ! % p k i i

$\frac{n!}{(n-m)!m!}\,\%p_i^{k_i}$ 的时候,由于分母与

pkii p i k i $p_i^{k_i}$ 不一定互质，没法求逆元，所以需要先把分子分母中

pi p i $p_i$ 全部除去，最后再乘回来。
代码实现：

long long x=0,ans;
for(long long i=n;i;i/=pi) x+=i/pi;
for(long long i=k;i;i/=pi) x-=i/pi;
for(long long i=n-k;i;i/=pi) x-=i/pi;
ans=a*mod_inverser(b,pk)%pk*mod_inverser(c,pk)*pow_mod(pi,x,pk)%pk;

解释一下 $\frac{i}{p_i}$ 的含义：
以求 $n!$ 为例：假设 $n=100,p_i=2$
$100!=100*99*98*...*1$
$=(99*97*95*...*1)*2^{50}*(1*2*3*...*50)$
$=(99*97*95*...*1)*2^{50}*(1*3*...*49)*2^{25}*(1*2*...*25)$
……
所以，一直除 $p_i$ 可将 $n!$ 表示为

p x i * a

$p_i^{x}*a$ a为其余的项的积。通过这个方法可求得

n!(n−m)!m!%pkii n ! ( n − m ) ! m ! % p i k i $\frac{n!}{(n-m)!m!}\,\%p_i^{k_i}$

代码：

#include <iostream>
#include <bits/stdc++.h>
using namespace std;

long long exgcd(long long a,long long b,long long &x,long long &y)
{
    long long d=a;
    if(b!=0)
    {
        d=exgcd(b,a%b,y,x);
        y-=(a/b)*x;
    }else{
        x=1;
        y=0;
    }
    return d;
}
long long pow_mod(long long a,long long b,long long p)
{
    long long res=1;
    while(b)
    {
        if(b&1) res=res*a%p;
        a=a*a%p;
        b>>=1;
    }
    return res;
}
long long mod_inverser(long long a,long long m)
{
    long long x,y;
    exgcd(a,m,x,y);
    return (m+x%m)%m;
}
long long cal(long long n,long long pi,long long pk)
{
    if(n==0) return 1;
    long long ans=1;
    if(n/pk)
    {
         for(long long i=2;i<=pk;i++)
         {
             if(i%pi)
                ans=ans*i%pk;
         }
         ans=pow_mod(ans,n/pk,pk);
    }
    for(long long i=2;i<=n%pk;i++)
    {
        if(i%pi)
            ans=ans*i%pk;
    }
    return ans*cal(n/pi,pi,pk)%pk;
}
long long C(long long n,long long k,long long p,long long pi,long long pk)
{
    if(k>n) return 0;
    long long a=cal(n,pi,pk),b=cal(k,pi,pk),c=cal(n-k,pi,pk);
    long long x=0,ans;
    for(long long i=n;i;i/=pi) x+=i/pi;
    for(long long i=k;i;i/=pi) x-=i/pi;
    for(long long i=n-k;i;i/=pi) x-=i/pi;
    ans=a*mod_inverser(b,pk)%pk*mod_inverser(c,pk)*pow_mod(pi,x,pk)%pk;
    return ans*(p/pk)%p*mod_inverser(p/pk,pk)%p;
}

int main()
{
    long long n,k,p;
    cin>>n>>k>>p;
    long long ans=0;
    for(long long x=p,i=2;i<=p;i++)
    {
        if(x%i==0)
        {
            long long pk=1;
            while(x%i==0)
            {
                pk*=i;
                x/=i;
            }
            ans=(ans+C(n,k,p,i,pk))%p;
        }
    }
    cout<<ans<<endl;
    return 0;
}