[SDOI2010]古时猪文

www.myexceptions.net  网友分享于：2013-02-20  浏览：17次
[SDOI2010]古代猪文

```#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <cctype>
#include <string>
#include <map>
#include <set>
#include <queue>
#include <algorithm>
#include <fstream>
#include <iostream>
using namespace std;

#ifdef WIN32
#define fmt64 "%I64d"
#else
#define fmt64 "%lld"
#endif
#define PI M_PI
#define oo 0x13131313
#define iter iterator
#define PB push_back
#define PO pop_back
#define MP make_pair
#define fst first
#define snd second
#define cstr(a) (a).c_str()

#define FOR(i, j, k) for(i = (j); i <= (k); ++i)
#define ROF(i, j, k) for(i = (j); i >= (k); --i)
#define FER(e, d, u) for(e = d[u]; e; e = e->n)
#define FRE(i, a) for(i = (a).begin(); i != (a).end(); ++i)

typedef unsigned int uint;
typedef long long int64;
typedef unsigned long long uint64;
typedef long double real;

template<class T> inline bool minim(T &a, const T &b) {return b < a ? a = b, 1 : 0;}
template<class T> inline bool maxim(T &a, const T &b) {return b > a ? a = b, 1 : 0;}
template<class T> inline T sqr(const T &a) {return a * a;}

#define maxn 100005
#define maxm 35700

int pt, pm[maxn / 5];
bool np[maxn + 1];
int ft, fac[20], cnt[20];
int dt, d[maxn];

void get_prime()
{
int i, j, k;
FOR(i, 2, maxn) {
if (!np[i]) pm[++pt] = i;
FOR(j, 1, pt) {
if ((k = i * pm[j]) > maxn) break;
np[k] = 1;
if (!(i % pm[j])) break;
}
}
}

void divide(int k)
{
int i, j; j = ft = 0;
FOR(i, 1, pt) {
if (sqr(pm[i]) > k) break;
if (!(k % pm[i])) {
fac[++ft] = pm[i], cnt[ft] = 1;
for (k /= pm[i]; !(k % pm[i]); k /= pm[i]) ++cnt[ft];
}
}
if (k != 1) fac[++ft] = k, cnt[ft] = 1;
}

void dfs(int t, int rest, int now)
{
if (!rest) {d[++dt] = now; return;}
if (t > ft) return;
dfs(t + 1, rest, now);
for (int i = 1; i <= cnt[t]; ++i) {
if (i > rest) break;
dfs(t + 1, rest - i, now *= fac[t]);
}
}

const int mm = 2 * 3 * 4679 * 35617;
const int m[4] = {2, 3, 4679, 35617};
const int M[4] = {499955829, 333303886, 213702, 28074};
int inv[4];

int f[4][maxm + 1], v[4][maxm + 1];

int ex_gcd(int a, int b, int &x, int &y)
{
if (!b) return x = 1, y = 0, a;
int ret = ex_gcd(b, a % b, x, y), t;
return t = x, x = y, y = t - (a / b) * y, ret;
}

int inverse(int a, int b)
{
int x, y;
ex_gcd(a, b, x, y);
if ((x %= b) < 0) x += b;
return x;
}

void prepare()
{
int i, j, k, l;
FOR(l, 0, 3) {
inv[l] = inverse(M[l], m[l]);
f[l][0] = f[l][1] = v[l][0] = v[l][1] = 1;
}
FOR(l, 0, 3) FOR(i, 2, m[l] - 1) {
f[l][i] = (f[l][i - 1] * i) % m[l];
if (!np[i]) v[l][i] = inverse(i, m[l]);
FOR(j, 1, pt) {
if ((k = pm[j] * i) > m[l]) break;
(v[l][k] = v[l][i] * v[l][pm[j]]) %= m[l];
if (!(i % pm[j])) break;
}
}
FOR(l, 0, 3) FOR(i, 2, m[l] - 1)
(v[l][i] *= v[l][i - 1]) %= m[l];
}

int C(int a, int b, int l)
{
return a < b ? 0 : (int64(f[l][a] * v[l][b]) * v[l][a - b]) % m[l];
}

int lucas(int a, int b, int l)
{
int64 ret = 1;
for (; b; a /= m[l], b /= m[l])
(ret *= C(a % m[l], b % m[l], l)) %= m[l];
return ret;
}

int main()
{
freopen("ancient.in", "r", stdin);
freopen("ancient.out", "w", stdout);

int i, j; int64 k; j = k = 0;

int n; int64 G; cin >> n >> G;
if (G == mm + 1) puts("0"), exit(0);
get_prime(), divide(n);
FOR(i, 1, ft) j += cnt[i];
FOR(i, 0, j) dfs(1, i, 1);

prepare();
FOR(i, 1, dt) {
int64 sum = 0;
FOR(j, 0, 3) {
int64 t = lucas(n, d[i], j);
(sum += inv[j] * M[j] * t) %= mm;
}
(k += sum) %= mm;
}

int64 ans = 1;
for (i = k; i; (G *= G) %= mm + 1, i >>= 1)
if (i & 1) (ans *= G) %= mm + 1;
cout << ans << endl;

return 0;
}
```

1楼a25201233天前 10:18