absi2011's Blog & Daily Life.

全新的开始       我要省选翻盘       I wanna AK in 高考\化学       自此,生无可恋
[破碎的状态] [-21] 100644 H
[破碎的状态] [-19] PA 2015 sia

[破碎的状态] [-21] 100633 J

absi2011 posted @ Jul 01, 2016 04:55:47 PM in 刷题记录 with tags 小高考 CF Gym 数学题 瞎搞 , 528 阅读

题意:

求[tex]C_n^k\ mod\ m[/tex]

[tex]n<=10^{18},0<=k<=n,m<=10^6[/tex]

解法:

众所周知..

[tex]C_n^k = \frac{n!}{k! (n-k)!}[/tex]

我们根据中国剩余定理..

首先..拆成若干个[tex]p^k[/tex]的形式

对于每个[tex]p^k[/tex],我们只要对于n,k和n-k知道(以下均以n为例):

1,[tex] n! = p^{a} * q[/tex]中a的值

2,[tex]q\ mod\ p^k[/tex]的值

这时候我们只要把这个a和q'和剩下两个值进行下操作,实际上是

[tex]a_{sum} = a_1 - a_2 - a_3[/tex]

[tex]q_{sum} = \frac {q_1}{q_2 q_3} = q_1 * q_2^{\phi(p^k)-1} * q_3^{\phi(p^k)-1}[/tex]

那么..

a我们可以通过不断除的办法来求..而q..

我们可以发现....

对于每个数x,[/tex]gcd(x,p)=1[/tex],则直接乘以[tex]x%p[/tex]

对于每个数[tex]x = p * a[/tex],我们要乘以[tex]x/p^2 % p[/tex]

对于每个数[tex]x = p^2 * a[/tex],我们要乘以[tex]x/p^2 % p[/tex]

......

我们可以把所有gcd(x,p)=1的给找出来,乘以一下,因为值太多,我们可以统计有多少个x % p = 1, x % p = 2 .....

这个数字是n/p或者n/p+1...直接乘以下就好了..

那么把剩下的数字/p,就转化回了原问题

考虑到我不会中国剩余定理的快速求的姿势..我暴力了一发.....写的不太好看

综上..

代码如下:

#include<set>
#include<map>
#include<list>
#include<queue>
#include<stack>
#include<string>
#include<time.h>
#include<math.h>
#include<memory>
#include<vector>
#include<bitset>
#include<fstream>
#include<stdio.h>
#include<utility>
#include<sstream>
#include<string.h>
#include<iostream>
#include<stdlib.h>
#include<algorithm>
using namespace std;
int modo;
int phi_modo;
long long power(long long x,long long y)
{
    if (y==0) return 1ll;
    long long t=power(x,y/2);
    t=t*t%modo;
    if (y%2==1) t=t*x%modo;
    return t;
}
int facts[1000005];
struct result
{
    long long modo_ans;
    long long cnt;
    result ()
    {
        modo_ans=1;
        cnt=0;
    }
    friend result operator * (const result &a,const result &b)
    {
        result c;
        c.modo_ans=(a.modo_ans*b.modo_ans)%modo;
        c.cnt=a.cnt+b.cnt;
        return c;
    }
    friend result operator / (const result &a,const result &b)
    {
        result c;
        c.modo_ans=(a.modo_ans*power(b.modo_ans,phi_modo-1))%modo;
        c.cnt=a.cnt-b.cnt;
        return c;
    }
};
result fact(long long x,int p)
{
    if (x==0)
    {
        return result();
    }
    if (x<modo)
    {
        result a;
        a.modo_ans=1;
        a.modo_ans=facts[x];
        a.cnt=x/p;
        return a*fact(x/p,p);
    }
    else
    {
        result a;
        a.modo_ans=1;
        int t=x%modo;
        a.modo_ans=facts[t];
        long long temp=facts[modo];
        a.cnt=x/p;
        a.modo_ans=a.modo_ans*power(temp,x/modo)%modo;
        return a*fact(x/p,p);
    }
}
int main()
{
    #ifdef absi2011
    freopen("input.txt","r",stdin);
    freopen("output.txt","w",stdout);
    #endif
    long long n,m,k;
    cin>>n>>m>>k;
    int i;
    int ans=1;
    int modo_ans=0;
    for (i=2;i<=k;i++)
    {
        if (k%i==0)
        {
            modo=1;
            int sum=0;
            for (;k%i==0;)
            {
                sum++;
                k/=i;
                modo*=i;
            }
            phi_modo=modo/i*(i-1);
            int j;
            facts[0]=1;
            for (j=1;j<=modo;j++)
            {
                facts[j]=facts[j-1];
                if (j%i!=0) facts[j]=facts[j-1]*(long long)j%modo;
            }
            result a=fact(n,i)/fact(m,i)/fact(n-m,i);
            if (a.cnt>=sum)
            {
                a.cnt=0;
                a.modo_ans=0;
            }
            else
            {
                a.modo_ans=a.modo_ans*power(i,a.cnt)%modo;
            }
            for (j=0;j<modo;j++)
            {
                if ((ans*(long long)j+modo_ans)%modo==a.modo_ans) break;
            }
            modo_ans=(ans*(long long)j+modo_ans)%(ans*modo);
            ans*=modo;
        }
    }
    printf("%d\n",modo_ans);
    return 0;
}

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter