1
0
正常的算法: 一》C(m, n) % p (p为素数, m, n <= 10^8) 算法伪代码:
1) 预处理 f[x] = (x!) mod p. (0 <= x < p)
2) 求m!中把p因子去掉,余下的数乘起来 mod p.
函数 G(m)求解 2)问题: 伪代码:
G(m) 1) if (0 == m) return 1; 2) m = k*p + t; 3) G[m] = f[p-1] ^ k * f[t] * G[m / k] mod p;
(其中 f[p-1]^k mod p 二分做)
3) C(m, n) = m! / (n! * (m-n)!); A = m! mod p; B = n!*(m-n)! mod p; C = A / B mod p; 在这里有俩种算法: 1) 求B的逆元; ( O(n^2) 枚举, 或者.. ) 2) 扩展欧几里德求 BC = A(mod p) (= 是同余的意思) 最小的 C 就是答案; 这里的 A, B 转换成 对 2)问题 求解.
#include <iostream>
#include <map>
#include <set>
#include <string>
#include <algorithm>
#include <math.h>
#include <queue>
using namespace std;
template <class T> void out(T x, int n){ for(int i = 0; i < n; ++i) cout << x[i] << ' '; cout << endl; }
template <class T> void out(T x, int n, int m){ for(int i = 0; i < n; ++i) out(x[i], m); cout << endl; }
#define OUT(x) (cout << #x << " = " << x << endl)
#define FOR(i, a, b) for(int i = (int)a; i < (int)b; ++i)
#define REP(i, b) FOR(i, 0, b)
#define FORD(i, a, b) for(int i = (int)a; i >= (int)b; --i)
#define N 10007
#define MOD(x, y) ((x) % (y) + y) % (y)
long long f[N+1];
long long ext_gcd(long long a, long long b, long long &x, long long &y) {
if (0 == b)
{
x = 1; y = 0;
return a;
}
else
{
long long tx, ty, d;
d = ext_gcd(b, MOD(a, b), tx, ty);
x = ty;
y = tx - a / b * ty;
return d;
}
}
long long modular_linear(long long a, long long b, long long n) {
long long tx, ty, d, x0, i;
d = ext_gcd(a, n, tx, ty);
if (0 == b % d)
{
x0 = MOD(tx * b / d, n);
return x0;
}
else
{
return 1;
}
}
void init(){
f[0] = 1;
FOR (i, 1, N)
{
f[i] = (f[i-1] * i) % N;
}
}
// m! 有多少p因子
long long calP(long long m, long long p){
long long ret = 0;
while (m)
{
ret += m / p;
m /= p;
}
return ret;
}
long long square(long long x){ return x * x; }
long long bigmod(long long b, long long p, long long m) {
if (0 == p)
return 1;
else if (0 == p % 2)
return square(bigmod(b, p/2, m)) % m;
else
return ((b % m) * bigmod(b, p-1, m)) % m;
}
// 计算m!-mod p
long long calM(long long m, long long p){
if (0 == m) return 1;
long long ret = 0, k, t;
k = m / p;
t = m % p;
ret = (bigmod(f[p-1], k, N) * f[t]) % p;
ret *= calM(m / p, p);
return ret % p;
}
long long calCmnModP(long long m, long long n, long long p){
long long A, B, u, v;
u = calP(m, p);
v = calP(n, p) + calP(m-n, p);
if (u > v) return 0;
else
{
A = calM(m, p);
B = (calM(n, p) * calM(m-n, p)) % p;
return modular_linear(B, A, p);
}
}
int main(){
int t, n, m;
scanf("%d", &t);
init();
while (t--)
{
scanf("%d %d", &m, &n);
printf("%d\n", calCmnModP(m, n, N));
}
return 0;
}
二》<h2 class="entry_title">大整数求组合数取余(Lucas定理)</h2><div class="archive_info"><span class="date">2013年02月25日</span><span class="category"> ⁄ 综合</span>⁄ 共 915字 ⁄ 字号 <span class="font"><a>小</a> <a>中</a> <a>大</a></span><span class="comment"> ⁄ <span>评论关闭</span></span><span class="edit"></span></div><div class="ad_r"></div><p>【卢卡斯(Lucas)定理】</p><p>Lucas定理用来求C(a,b)mod p的值,其中p为素数。</p><p>数学表达式为:</p><p>Lucas(a,b,q)=C(a%q,b%q)*Lucas(a/p,b/p,p);</p><p>Lucas(a,0,q)=0;</p><p>通过这个定理就可以很方便的把大数的组合转化成小数。但其中还是要求C(a%q,b%q)%p,所以这里引入逆元来求。</p><p>【定义】若整数a,b,p, 满足a·b≡1(mod p).则称a 为b 模p 的乘法逆元, 即a=b- 1mod p.其中, p 是模数。</p><p>应用到组合数中来就是:</p><p> a!/[b!*(a-b)!] % p == a! * [b!*(a-b)!]-1 %p</p><p>【逆元求法】:</p><p>应用费马小定理,ap-1=1 mod p ,即 a*ap-2=1 mod p</p><p>也就是说 ap-2就是a的逆元。</p><p>当然这里求出来的逆元是在取模p的逆元,对我们最终目标没有影响。这也是比较方便而且比较好的方法。</p><p>弄个模板来。</p><p>
</p><p>
</p><p><pre name="code" class="cpp">#include <cstdio>
#include <iostream>
#include <cmath>
#include <cstring>
#include <algorithm>
using namespace std;
#define maxn 100010
typedef long long LL;
LL m,n,p;
LL Pow(LL a,LL b,LL mod)
{
LL ans=1;
while(b)
{
if(b&1)
{
b--;
ans=(ans*a)%mod;
}
else
{
b/=2;
a=(a*a)%mod;
}
}
return ans;
}
LL C(LL n,LL m)
{
if(n<m)
return 0;
LL ans=1;
for(int i=1;i<=m;i++)
{
ans=ans*(((n-m+i)%p)*Pow(i,p-2,p)%p)%p;
}
return ans;
}
LL Lucas(LL n,LL m)
{
if(m==0)
return 1;
return (Lucas(n/p,m/p)*C(n%p,m%p))%p;
}
int main()
{
int t;
scanf("%d",&t);
while(t--)
{
scanf("%lld%lld%lld",&n,&m,&p);
printf("%lld\n",Lucas(n,m));
}
return 0;
}
三》<h2 id="t_618194800100l5rv" class="titName SG_txta">大数取余</h2>
<div id="sina_keyword_ad_area2" class="articalContent "><p>有一类题目会因为求出的结果太大而只要求输出对某个数m取余后的结果,而且这个m是比较小的数,比如不超过32位整数…
而这类大数都是可以由较小的数经过某些运算得到的…
于是我整理了一下对付几种运算的方法…包括四则运算,指数,组合数,塔函数的应对方法…</p><p>那么就开始吧…慢慢来…
首先是最常识的加减法:</p><pre>add_mod(a,b,m){
return ((a%m)+(b%m))%m;
}
别小看加法哦…很多用dp解的题目中靠着加法可是能达到很大的数呢… minus_mod(a,b,m){
return (a-b+m)%m;
}
减法…会遇到吗? 接着是依然很简单的乘法: multiply_mod(a,b,m){
return ((a%m)*(b%m))%m;
}
这是当m*m不会溢出时可以用的,同时也是通常的情况… 但是如果m*m连long long都会溢出的话…就需要把一个数一位位拆开来做了… multiply_mod(a,b,m){
if(b==0)return 0;
return (((b&1)?a:0)+(multiply_mod(a,b>>1,m)<<1)%m)%m;
}
然后是除法,但有点限制: (a/b)%m 特殊条件:m和b互质 前提:a能被b整除 这个有点特殊,意为虽然不知道a是多少,但是已知c,而且c=a%m,用c和b来求(a/b)%m的方法 虽然需要m和b互质,但是不互质的话也是可以做的,因为a也一定是gcd(b,m)的倍数,具体可以看看这里 需要用到扩展欧几里德来求… 至于扩展欧几里德是什么…去Google一下吧… extend_euclid(a,b,&x0,&y0){
if(b==0){
x0=1;
y0=0;
return a;
}
r=extend_euclid(b,a%b);
t=x0;
x0=y0;
y0=t-a/b*y0;
return r;
}
divide_mod(a,b,m){
extend_euclid(b,m,x0,y0);
return (a*(((x0%m)+m)%m))%m;
}
这个在求组合数的时候可能用到… 不过似乎很少遇到需要用除法取余的情况呢… 然后是更大却很简单的幂运算: (a^b)%m 这是初学递归或者二分时就会遇上的一个很简单的方法,和之前的乘法差不多 power_mod(a,b,m){
if(b==0)return 1;
if(b&1)return (a*power_mod((a*a)%m,b>>1,m))%m;
else return power_mod((a*a)%m,b>>1);
}
其中乘法取余会运算中溢出的话可以改成之前那个multiply_mod() 恩…开始有趣了…下面是组合数: C(a,b)%m 特殊条件:m是质数 如果b或者a-b比较小的话,可以用之前计算(a/b)%m的方式来把组合数公式展开来计算 不过当b和a-b与m相比很大时,有更好的方法: a,b在m进制下表示为 a=(ak,…,a0),b=(bk,…,b0)0=<ai,bi<m 于是会有这样的性质: C(a,b)=C(ak,bk)*...*C(a0,b0) (mod m)
最后是难以想象的大的数…塔函数: (a↑↑b)%m 这里可以看到其一些性质 比如ProjectEuler 282 Ackermann 函数就是非常恶心的大数 第一层是很小的常数,第二层是n个a相加,也就是乘法,第三层是n个a相乘,也就是幂,第四层是n个a叠着做幂即塔函数,第k+1层是n个a做第k层运算… 用小数字居然也能表示出如此之大的数…佩服啊… 在上一篇中也有提到: n>=phi(m)时,a^n=a^(n%phi(m)+phi(m)) (mod m)
其中phi()是欧拉函数 由于phi(x)<x,所以就算是对a↑↑b这样大的数也总会在足够短的时间内收敛,然后计算出(modm)的值 这个不止在塔函数中,也可以用在各种指数异常大的情况下 特别的,在b>m的情况下,(a↑↑b)%m的值将是定值 恩…就到这里了…
PS:逆元素-http://baike.baidu.com/view/3140397.htm?from_id=11054145&type=syn&fromtitle=%E9%80%86%E5%85%83&fr=aladdin
非正常的思路: (找到了余数的规律 = = 。。)
代码还过了。。。 卢卡思、费马、欧几里德怕不是要死的活过来!??? #include<stdio.h>
int main()
{
int i,n,p,t;
while (scanf("%d", &n) != EOF)
{
for (i = 0; i < n; i++)
{
scanf("%d%d", &t, &p);
printf("%d\n", (t / p) % p);
}
}
return 0;
}
|