We use cookies to ensure you have the best browsing experience on our website. Please read our cookie policy for more information about how we use cookies.
importfunctoolsimportsysclassMathEngine():"""A collection of algorithms for Modular Arithmetics.References:[1]Granville,A."Binomial coefficients modulo prime powers".[2]Knuth,D."The Art of Computer Programming",vol1.[3]https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Polynomial_extended_Euclidean_algorithm"""def__init__(self):self.pe=Noneself._cache_factorial={}defbinomial(self,n,m):"""Find the binomial (n m) = n!/m!/(n-m)!."""ifm==0:ifn==0:return0else:return1b=self.factor_binomial(n,m)result=self.factors_to_integer(b)returnresultdefbinomial_prime_power(self,n,m,p):"""Find the max power of `p` which divides a binomial (n m)."""result=self._number_of_carries(n-m,m,p)returnresultdeffactor_binomial(self,n,m):"""Factor the binomial (n m) into primes.ThealgorithmisbasedonKummer'stheorem[1,page7].Result:{d1:p1,...,dn:pn},sothat(nm)=d1^p1*...*dn^pn."""result={}forpinself.primes(n+1):power=self.binomial_prime_power(n,m,p)ifpower>0:result[p]=powerreturnresultdeffactor_integer(self,n):"""Factor n into primes.Result:{d1:p1,...,dn:pn},sothatn=d1^p1*...*dn^pn."""result={}forpinself.primes(n+1):power=0whilen%p==0andn!=1:n=n// ppower+=1ifpower>0:result[p]=powerifn==1:breakreturnresultdeffactor_factorial(self,n):"""Factor the factorial of n into primes.ThealgorithmisbasedontheKummer'stheoremformulatedforbinomialcoefficients;forthecaseoffactorialssee[1,eq.(17)]or[2,eq.(1.2.5.8)].Result:{d1:p1,...,dn:pn}suchthatn!=d1^p1*...*dn^pn."""result={}forprimeinself.primes(n+1):k=nresult[prime]=0whilek>0:k=k// primeresult[prime]+=kreturnresultdeffactorial(self,n,p=0,mod=0):"""Find the product of all integers <= n, i.e. the factorial of n;ifp>0excludeintegersdivisiblebyp;ifmod>0findresultmodulomod."""_cache_key=(n,p,mod)if_cache_keyinself._cache_factorial:returnself._cache_factorial[_cache_key]result=1ifp==0:ifmod==0:forkinrange(1,n+1):result*=kelse:forkinrange(1,n+1):result=(result*k)%modelse:ifmod==0:forkinrange(1,n+1):ifk%p!=0:result*=kelse:forkinrange(1,n+1):ifk%p!=0:result=(result*k)%modself._cache_factorial[_cache_key]=resultreturnresultdeffactors_to_integer(self,n):result=1forprime,powerinn.items():result*=prime**powerreturnresultdefgcd_extended(self,a,b):"""Find gcd(a,b) and coefficients s,t of the Bezout's identity, such that`a*s+b*t=gcd(a,b)`.See[3]formoredetails.Result:(gcd(a,b),s,t)"""s,s_old=0,1t,t_old=1,0r,r_old=b,awhiler!=0:k=r_old// rr_old,r=r,r_old-k*rs_old,s=s,s_old-k*st_old,t=t,t_old-k*treturnr_old,s_old,t_olddefinteger_digits(self,n,p):"""Find digits of n in the p representation.Result:[a0,...,ak],suchthatn=a0+a1*p+...+ak*p^k,whereai<p."""result=[]whilen>0:result.append(n%p)n//= preturnresultdefinteger_digits_len(self,n,p):result=1whilen>0:result+=1n//= preturnresultdefmod_binomial(self,n,m,mod):"""Find the reminder from the division of the binomial (n m) by mod.Thealgorithmissimpleandslow,butlesspronetoerrorscomparingtothefasterandmorecomplexgeneralizedLucasalgorithm,formoredetailssee`mod_binomial_lucas`method."""binomial_factors=self.factor_binomial(n,m)result=1forprime,expinbinomial_factors.items():result=(result*pow(prime,exp,mod))%modreturnresultdefmod_binomial_lucas(self,n,m,p,q=1):"""Find the reminder form the division of the binomial (n m) by p^q,wheremodisaprimenumber.ThisroutineemploysLucas'Theorem(forq=1)anditsgeneralizationfor(q>1),see[1,Theorem1]."""e0=self.binomial_prime_power(n,m,p)qq=q-e0ifqq<0:return0p_e0,p_q,p_qq=p**e0,p**q,p**qqr=n-md=self.integer_digits_len(n,p)_ni,_mi,_ri=n,m,r_least_positive_residue_n=[]_least_positive_residue_m=[]_least_positive_residue_r=[]foriinrange(0,d):_least_positive_residue_n.append(_ni%p_qq)_least_positive_residue_m.append(_mi%p_qq)_least_positive_residue_r.append(_ri%p_qq)_ni//= p_mi//= p_ri//= p_factorial=[1]curr=1forxinrange(1,p_qq):ifx%p!=0:curr=(curr*x)%p_qq_factorial.append(curr)result,i=1,0aaa,bbb,ccc=1,1,1foriinrange(0,d):aa=_least_positive_residue_n[i]bb=_least_positive_residue_m[i]cc=_least_positive_residue_r[i]a=_factorial[aa]b=_factorial[bb]c=_factorial[cc]aaa=(aaa*a)%p_qqbbb=(bbb*b)%p_qqccc=(ccc*c)%p_qqi+=1result=aaa*self.mod_inverse(bbb*ccc%p_qq,p_qq)%p_qqifp!=2orqq<3:result*=(-1)**self._number_of_carries(m,r,p,qq-1)result=(result*p_e0)%p_qreturnresultdef_least_positive_residue(self,n,j,p,q):"""Helper for mod_binomial_lucas method"""result=(n//p**j) % p**qreturnresultdef_number_of_carries(self,m,r,p,j=0):"""Helper for mod_binomial_lucas method"""result,carry,jj=0,0,0while(m>0orr>0):if(m%p)+(r%p)+carry>=p:carry=1else:carry=0ifjj>=j:result+=carrym,r=m// p, r // pjj+=1returnresultdefmod_chinese(self,rs):n=functools.reduce(lambdaa,b:a*b,rs.keys(),1)result=0forni,aiinrs.items():nn=n// ni_,ri,si=self.gcd_extended(ni,nn)ei=si*nnresult+=ai*eiresult=result%nreturnresultdefmod_inverse(self,n,mod):r,s,t=self.gcd_extended(n,mod)ifr!=1:returnNonereturnsdefprimes(self,n):returnself.pe.primes(n)classSolver(object):def__init__(self):self.me=MathEngine()self.mod=142857defsolution(self,n,r):"""Find solution to the "nCr" problem."""mod_factors={3:3,11:1,13:1,37:1}mods={p**q:self.me.mod_binomial_lucas(n,r,p,q)forp,qinmod_factors.items()}result=self.me.mod_chinese(mods)returnresultif__name__=='__main__':T=int(sys.stdin.readline())solver=Solver()foriinrange(T):n,r=[int(num)fornuminsys.stdin.readline().split()]print(solver.solution(n,r))
nCr
You are viewing a single comment's thread. Return to all comments →
Python3 solution