实数gcd，大数快速乘与组合数取模

实数gcd

在计算几何的题目中，有时会需要你算两个角度的最大公因数(会有误差)，具体操作和整数的gcd类似，但要注意输入可能会有0。

模板如下：

double fgcd(double a,double b)
{
if(fabs(a)<eps) return b;
if(fabs(b)<eps) return a;
return fgcd(b,fmod(a,b));
}//fmod为math库的小数取余函数

大数快速乘

适用于那些中间变量会越界的乘法，快速乘和快速幂的想法基本类似，相当于降一级运算。

模板如下：

LL pow_mul(LL a,LL n,LL m)
{
if(n == 0) return 1;
LL x = pow_mul(a,n>>1,m);
LL ans = (x+x)%m;
if(n&1) ans = (ans+a)%m;
return ans;
}

组合数取模：逆元+lucas定理+CRT+extend_欧几里得算法

题目：2015长春网络赛hdoj5446

地址：http://acm.hdu.edu.cn/showproblem.php?pid=5446

代码：

#include <iostream>
#pragma comment(linker, "/STACK:1024000000,1024000000") 
#include <stdio.h>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <string>
#include <string.h>
#include <sstream>
#include <cctype>
#include <climits>
#include <set>
#include <map>
#include <deque>
#include <queue>
#include <vector>
#include <iterator>
#include <algorithm>
#include <stack>
#include <functional>
/*int类型最大值INT_MAX，short最大值为SHORT_MAX
long long最大值为LONG_LONG_MAX*/
//cout << "OK" << endl;
#define _clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int INF = INT_MAX;
const double eps = 1e-8;
const double EULER = 0.577215664901532860;
const double PI = 3.1415926535897932384626;
const double E = 2.71828182845904523536028;
typedef long long LL;

int dir[4][2]={{0,-1},{0,1},{1,0},{-1,0}};
LL p[20],ans[20] = {0},n,m,k,t; 
LL pow_mod(LL a,LL n,LL m)
{
	if(n == 0) return 1;
	LL x = pow_mod(a,n>>1,m);
	LL ans = x*x%m;
	if(n&1) ans = ans*a%m;
	return ans;
}
LL f(LL m,LL n,LL p)
{
	if(m > n) return 0;
	LL ans = 1;
	for(int i=1; i<=m; i++)
    {
        LL a = (n + i - m) % p;
        LL b = i % p;
        ans=ans*(a * pow_mod(b, p-2,p) % p) % p;
        //cout << a << " " << pow_mod(b, p-2,p) << endl;
    }
    return ans;
}
LL Lucas(LL m, LL n,LL p)
{
    if(m == 0) return 1;
    return f(m % p, n % p,p) * Lucas(m / p, n / p,p) % p;
}
LL extend_gcd(LL a,LL b,LL &x,LL &y)
{
	if(!b)
	{
		x = 1;y = 0;
		return a;
	}
	else
	{
		LL r = extend_gcd(b,a%b,y,x);
		y-=x*(a/b);
		return r;
	}
}
LL multi(LL a,LL b,LL p)   //大数 （a×b）%p
{
    LL ret=0;
    while(b)
    {
        if(b&1)
            ret=(ret+a)%p;
            a=(a+a)%p;
            b>>=1;
    }
    return ret;

}
LL CRT()
{
	LL M = 1,ret = 0;
	for(int i = 0;i<k;i++)M*=p[i];
	for(int i = 0;i<k;i++)
	{
		LL x,y;
		LL tm = M/p[i];
		extend_gcd(tm,p[i],x,y);
		//cout << tm << " " << x << endl;
		ret = (ret+multi(multi(x,tm,M),ans[i],M))%M;
	}
	return (ret+M)%M;
}
int main()
{
    //freopen("sample.in", "r", stdin);
	//freopen("sample.out", "w", stdout);
	
	scanf("%I64d",&t);
	while(t--)
	{
		_clr(ans,0);
		scanf("%I64d%I64d%I64d",&n,&m,&k);
		for(int i = 0;i<k;i++)
		{
			scanf("%I64d",&p[i]);
			ans[i]= Lucas(m,n,p[i]);
			//cout << ans[i] << endl;
		}
		printf("%I64d\n",CRT());
	}

    //fclose(stdin);
    //fclose(stdout); 
    return 0;
}