[Bzoj 3129] [Sdoi2013] equation "tolerance + combination number modulo + Chinese remainder theorem"

Source: Internet
Author: User
Tags modulus

Title Link: BZOJ-3129

Problem analysis

The idea of using the partition method, if there are no restrictions, then the program number is C (m-1, n-1).

If there is a limitation of Xi >= Ai, then we can subtract the M from Ai-1, which is equivalent to a fixed portion of this part to Xi, which translates into unrestricted conditions.

What if there are some restrictions on the XI <= Ai? It is not possible to ask directly, but notice that there are no more than 8 limits, and we can use the principle of repulsion.

Consider the allowance: consider which restrictions are violated, that is, what restrictions are XI <= ai but XI > ai, which translates to the constraints of the XI >= AI.

Then we can find out the answer in the time of 2^8 * T (combinatorial number).

How do you find the number of combinations? It is not possible to directly preprocess the inverse of factorial, because modulus is not always prime.

We want to divide the modulus into a form such as Pi^ai, so that they are between 22 coprime, we can find out the answer respectively, and then use the Chinese remainder theorem combined together.

Chinese remainder theorem: if there are n equations x = XI (mod mi), M = m1 * m2 *. * mn, in the sense of MoD M, there is a unique solution to the equation set.

x = Sigma (MI * INV (MI) * xi)% M, where mi = M/mi, INV (MI) is the inverse of mi in mod mi sense.

So our question is, how to find C (n, m)% (P^a).

This requires the use of "combined number modulo", which is specifically used to solve this problem.

Using a method similar to fast factorial, the factorial of the bar up and down in the combined number is split into the form E * p^f, then the e is calculated directly, and the F bar is subtracted and then calculated.

How to x! and split it into E * p^f?

Assuming that we want the MoD number to be p^a, then we need to preprocess the prefix product of the remainder of the multiples of p in [1, P^a-1] (similar to factorial less than multiples of p).

Then we know [1, X] contains p in the number of x/p, we will extract 1 p in these numbers, then we get p^ (x/p), and then this x/p number becomes [1, x/p], you can recursively go down.

The rest of the sections can be segmented to be divided into [1, P^a-1], [P^a + 1, p^a + p^a-1] ..... In this way, the product of each paragraph is the same (mod p^a sense), the direct fast power can be.

Finally, there will be a section [1, x% (P^a)], which is also a direct preprocessing value.

So this problem is done (call ~).

Also note that when writing code, I ask the inverse to use Euler's theorem but it is wrong.

Euler's theorem: A^phi (b) = 1 (mod b) Condition: GCD (A, b) = 1

Note that it is A^phi (b) and not a^ (b-1) ! When B is not prime, kneel!

Code
#include <iostream> #include <cstdlib> #include <cstring> #include <algorithm> #include < Cmath> #include <cstdio>using namespace std;typedef long long ll;typedef double lf;const int maxp = 10201 +, Ma XN1 = 8 + 5;int T, p, N, N1, N2, M, Top, Ans;int a[maxn1]; LL Temp; LL Fac[10][maxp], pr[10], pi[10], pa[10], phi_pi[10], mi[10], inv_mi[10], xi[10]; ll Pow (ll A, ll B, ll Mod) {ll ret, f;ret = 1; f = a;while (b) {if (b & 1) {ret *= f;ret%= Mod;} b >>= 1;f *= f;f%= Mod;} return ret;} void Prepare () {int x, sqrtx;x = p; SQRTX = (int) sqrt ((LF) x); Top = 0;for (int i = 2; I <= sqrtx; ++i) {if (x% i! = 0) Continue; Pr[++top] = i; Pa[top] = 0; Pi[top] = 1;while (x% i = = 0) {++pa[top]; Pi[top] *= i;x/= i;} Phi_pi[top] = Pi[top]/pr[top] * (Pr[top]-1);} if (x > 1) {pr[++top] = x; Pa[top] = 1; Pi[top] = x; Phi_pi[top] = pi[top]-1;} for (int i = 1; I <= Top; ++i) {mi[i] = p/pi[i];inv_mi[i] = Pow (Mi[i], Phi_pi[i]-1, pi[i]); Fac[i][0] = 1;for (int j = 1; J &LT Pi[i]; ++J) {if (j% pr[i]! = 0) Fac[i][j] = fac[i][j-1] * J% pi[i];else fac[i][j] = fac[i][j-1];}}} struct Es{ll e, f;}; ES Calc (int x, int k) {ES ret, tc;if (x < pr[k]) {ret.e = FAC[K][X];RET.F = 0;return ret;} RET.F = X/PR[K];TC = Calc (X/pr[k], k); Ret.f + = TC.F;RET.E = TC.E * fac[k][x% pi[k]]% PI[K];RET.E = RET.E * Pow (fac[k ][PI[K]-1], x/pi[k], pi[k])% pi[k];return ret;} ll C (int x, int y, int k) {ll ret;int pf;es Ex, Ey, Exy; Ex = Calc (x, k); Ey = Calc (y, k); Exy = Calc (X-y, k); ret = EX.E * Pow (EY.E, Phi_pi[k]-1, pi[k])% pi[k] * Pow (EXY.E, Phi_pi[k]-1, pi[k])% pi[k];p f = Ex.f-ey.f-exy.f;if (PF >= pa[k]) ret = 0;ELSE ret = ret * POW (Pr[k], pf, pi[k])% pi[k];return ret; int C (int x, int y) {if (x = = y) return 1;if (x < y) return 0;if (y = = 0) return 1;  LL ret = 0;for (int i = 1; I <= top; ++i) xi[i] = C (x, y, i); for (int i = 1; I <= Top; ++i) {ret + = xi[i] * mi[i]% P * Inv_mi[i]% P;ret%= p;} return (int) ret; void DFS (int x, int Cnt, int Sum) {if (x = = N1) {if (Cnt & 1) Temp-= C (m-sum-1, n-1); else Temp + = C (m-sum-1, n-1); temp = (temp% p + p)% P;return; }dfs (x + 1, CNT, sum);D FS (x + 1, cnt + 1, sum + a[x + 1]);} int Solve () {temp = 0;dfs (1, 0, 0);D FS (1, 1, a[1]); return (int) temp;} int main () {scanf ("%d%d", &t, &p); Prepare (); for (int case = 1; Case <= T; ++case) {scanf ("%d%d%d%d", &n, &n1, &n2, &m); for (int i = 1; I <= N1; ++i) scanf ("%d", &a[i]); int Nu m;for (int i = 1; I <= n2; ++i) {scanf ("%d", &num); m-= Num-1;} if (N1 > 0) ans = Solve (); Else ans = C (m-1, n-1);p rintf ("%d\n", ans);} return 0;}

  

[Bzoj 3129] [Sdoi2013] equation "tolerance + combination number modulo + Chinese remainder theorem"

Contact Us

The content source of this page is from Internet, which doesn't represent Alibaba Cloud's opinion; products and services mentioned on that page don't have any relationship with Alibaba Cloud. If the content of the page makes you feel confusing, please write us an email, we will handle the problem within 5 days after receiving your email.

If you find any instances of plagiarism from the community, please send an email to: info-contact@alibabacloud.com and provide relevant evidence. A staff member will contact you within 5 working days.

A Free Trial That Lets You Build Big!

Start building with 50+ products and up to 12 months usage for Elastic Compute Service

  • Sales Support

    1 on 1 presale consultation

  • After-Sales Support

    24/7 Technical Support 6 Free Tickets per Quarter Faster Response

  • Alibaba Cloud offers highly flexible support services tailored to meet your exact needs.