class Modulo_Polynomial(): __slots__=("Poly", "max_degree") def __init__(self, Poly=[], max_degree=2*10**5): """ 多項式の定義 P: 係数のリスト max_degree ※Mod: 法はグローバル変数から指定 """ if Poly: self.Poly=[p%Mod for p in Poly[:max_degree]] else: self.Poly=[0] self.max_degree=max_degree def __str__(self): return str(self.Poly) def __repr__(self): return self.__str__() def __iter__(self): yield from self.Poly #= def __eq__(self,other): from itertools import zip_longest return all([a==b for a,b in zip_longest(self.Poly,other.Poly,fillvalue=0)]) #+,- def __pos__(self): return self def __neg__(self): return self.scale(-1) #items def __getitem__(self, index): if isinstance(index, slice): return Modulo_Polynomial(self.Poly[index], self.max_degree) else: if index<0: raise IndexError("index is negative (index: {})".format(index)) elif index>=len(self.Poly): return 0 else: return self.Poly[index] def __setitem__(self, index, value): if index<0: raise IndexError("index is negative (index: {})".format(index)) elif index>=self.max_degree: return if index>=len(self.Poly): self.Poly+=[0]*(index-len(self.Poly)+1) self.Poly[index]=value%Mod #Boole def __bool__(self): return any(self.Poly) #簡略化 def reduce(self): """ 高次の 0 を切り捨て """ P=self.Poly for d in range(len(P)-1,-1,-1): if P[d]: break self.resize(d+1) return #シフト def __lshift__(self,other): if other<0: return self>>(-other) if other>self.max_degree: return Modulo_Polynomial([0],self.max_degree) G=[0]*other+self.Poly return Modulo_Polynomial(G,self.max_degree) def __rshift__(self,other): if other<0: return self<<(-other) return Modulo_Polynomial(self.Poly[other:],self.max_degree) #次数 def degree(self): P=self.Poly for d in range(len(self.Poly)-1,-1,-1): if P[d]: return d return -float("inf") #加法 def __add__(self,other): P=self; Q=other if Q.__class__==Modulo_Polynomial: N=min(P.max_degree,Q.max_degree) A=P.Poly; B=Q.Poly else: N=P.max_degree A=P.Poly; B=Q return Modulo_Polynomial(Calc.Add(A,B),N) def __radd__(self,other): return self+other #減法 def __sub__(self,other): P=self; Q=other if Q.__class__==Modulo_Polynomial: N=min(P.max_degree,Q.max_degree) A=P.Poly; B=Q.Poly else: N=P.max_degree A=P.Poly; B=Q return Modulo_Polynomial(Calc.Sub(A,B),N) def __rsub__(self,other): return (-self)+other #乗法 def __mul__(self,other): P=self Q=other if Q.__class__==Modulo_Polynomial: a=b=0 for x in P.Poly: if x: a+=1 for y in Q.Poly: if y: b+=1 if a>b: P,Q=Q,P P.reduce();Q.reduce() U,V=P.Poly,Q.Poly M=min(P.max_degree,Q.max_degree) if a<2*P.max_degree.bit_length(): B=[0]*(len(U)+len(V)-1) for i in range(len(U)): if U[i]: for j in range(len(V)): B[i+j]+=U[i]*V[j] if B[i+j]>Mod: B[i+j]-=Mod else: B=Calc.Convolution(U,V)[:M] B=Modulo_Polynomial(B,M) B.reduce() return B else: return self.scale(other) def __rmul__(self,other): return self.scale(other) #除法 def __floordiv__(self,other): if not other: raise ZeroDivisionError if isinstance(other,int): return self/other self.reduce() other.reduce() return Modulo_Polynomial(Calc.Floor_Div(self.Poly, other.Poly), max(self.max_degree, other.max_degree)) def __rfloordiv__(self,other): if not self: raise ZeroDivisionError if isinstance(other,int): return Modulo_Polynomial([],self.max_degree) #剰余 def __mod__(self,other): if not other: return ZeroDivisionError self.reduce(); other.reduce() r=Modulo_Polynomial(Calc.Mod(self.Poly, other.Poly), min(self.max_degree, other.max_degree)) r.reduce() return r def __rmod__(self,other): if not self: raise ZeroDivisionError r=other-(other//self)*self r.reduce() return r def __divmod__(self,other): q=self//other r=self-q*other; r.reduce() return (q,r) #累乗 def __pow__(self,other): if other.__class__==int: n=other m=abs(n) Q=self A=Modulo_Polynomial([1],self.max_degree) while m>0: if m&1: A*=Q m>>=1 Q*=Q if n>=0: return A else: return A.inverse() else: P=Log(self) return Exp(P*other) #逆元 def inverse(self, deg=None): assert self.Poly[0], "定数項が0" if deg==None: deg=self.max_degree return Modulo_Polynomial(Calc.Inverse(self.Poly, deg), self.max_degree) #除法 def __truediv__(self,other): if isinstance(other, Modulo_Polynomial): if Calc.is_sparse(other.Poly): d,f=Calc.coefficients_list(other.Poly) K=len(d) H=[0]*self.max_degree alpha=pow(other[0], Mod-2, Mod) H[0]=alpha*self[0]%Mod for i in range(1, self.max_degree): c=0 for j in range(1, K): if d[j]<=i: c+=f[j]*H[i-d[j]]%Mod else: break c%=Mod H[i]=alpha*(self[i]-c)%Mod H=Modulo_Polynomial(H, min(self.max_degree, other.max_degree)) return H else: return self*other.inverse() else: return pow(other,Mod-2,Mod)*self def __rtruediv__(self,other): return other*self.inverse() #スカラー倍 def scale(self, s): return Modulo_Polynomial(Calc.Times(self.Poly,s),self.max_degree) #最高次の係数 def leading_coefficient(self): for x in self.Poly[::-1]: if x: return x return 0 def censor(self, N=-1, Return=False): """ N 次以上の係数をカット """ if N==-1: N=self.max_degree N=min(N, self.max_degree) if Return: return Modulo_Polynomial(self.Poly[:N],self.max_degree) else: self.Poly=self.Poly[:N] def resize(self, N, Return=False): """ 強制的に Poly の配列の長さを N にする. """ N=min(N, self.max_degree) P=self if Return: if len(P.Poly)>N: E=P.Poly[:N] else: E=P.Poly+[0]*(N-len(P.Poly)) return Modulo_Polynomial(E,P.max_degree) else: if len(P.Poly)>N: del P.Poly[N:] else: P.Poly+=[0]*(N-len(P.Poly)) #代入 def substitution(self, a): """ a を (形式的に) 代入した値を求める. a: int """ y=0 t=1 for p in self.Poly: y=(y+p*t)%Mod t=(t*a)%Mod return y #================================================= class Calculator: def __init__(self): self.primitive=self.__primitive_root() self.__build_up() def __primitive_root(self): p=Mod if p==2: return 1 if p==998244353: return 3 if p==10**9+7: return 5 if p==163577857: return 23 if p==167772161: return 3 if p==469762049: return 3 fac=[] q=2 v=p-1 while v>=q*q: e=0 while v%q==0: e+=1 v//=q if e>0: fac.append(q) q+=1 if v>1: fac.append(v) g=2 while g>rank2, Mod) iroot[-1]=pow(root[-1], Mod-2, Mod) for i in range(rank2)[::-1]: root[i]=root[i+1]*root[i+1]%Mod iroot[i]=iroot[i+1]*iroot[i+1]%Mod prod=iprod=1 for i in range(rank2-1): rate2[i]=root[i+2]*prod%Mod irate2[i]=iroot[i+2]*prod%Mod prod*=iroot[i+2]; prod%=Mod iprod*=root[i+2]; iprod%=Mod prod=iprod = 1 for i in range(rank2-2): rate3[i]=root[i + 3]*prod%Mod irate3[i]=iroot[i + 3]*iprod%Mod prod*=iroot[i + 3]; prod%=Mod iprod*=root[i + 3]; iprod%=Mod self.root=root; self.iroot=iroot self.rate2=rate2; self.irate2=irate2 self.rate3=rate3; self.irate3=irate3 def Add(self, A, B): """ 必要ならば末尾に元を追加して, [A[i]+B[i]] を求める. """ if type(A)==int: A=[A] if type(B)==int: B=[B] m=min(len(A), len(B)) C=[(A[i]+B[i])%Mod for i in range(m)] C.extend(A[m:]) C.extend(B[m:]) return C def Sub(self, A, B): """ 必要ならば末尾に元を追加して, [A[i]-B[i]] を求める. """ if type(A)==int: A=[A] if type(B)==int: B=[B] m=min(len(A), len(B)) C=[0]*m C=[(A[i]-B[i])%Mod for i in range(m)] C.extend(A[m:]) C.extend([-b%Mod for b in B[m:]]) return C def Times(self,A, k): """ [k*A[i]] を求める. """ return [k*a%Mod for a in A] #参考元 https://judge.yosupo.jp/submission/72676 def NTT(self, A): """ A に Mod を法とする数論変換を施す ※ Mod はグローバル変数から指定 References: https://github.com/atcoder/ac-library/blob/master/atcoder/convolution.hpp https://judge.yosupo.jp/submission/72676 """ N=len(A) H=(N-1).bit_length() l=0 I=self.root[2] rate2=self.rate2; rate3=self.rate3 while l=2: i=Q.popleft(); j=Q.popleft() A[i]=self.Convolution(A[i], A[j]) A[j]=None Q.append(i) i=Q.popleft() return A[i] def Inverse(self, F, length=None): if length==None: M=len(F) else: M=length if M<=0: return [] if self.is_sparse(F): """ 愚直に漸化式を用いて求める. 計算量: F にある係数が非零の項の個数を K, 求める最大次数を N として, O(NK) 時間 """ d,f=self.coefficients_list(F) G=[0]*M alpha=pow(F[0], Mod-2, Mod) G[0]=alpha for i in range(1, M): for j in range(1, len(d)): if d[j]<=i: G[i]+=f[j]*G[i-d[j]]%Mod else: break G[i]%=Mod G[i]=(-alpha*G[i])%Mod del G[M:] else: """ FFTの理論を応用して求める. 計算量: 求めたい項の個数をNとして, O(N log N) Reference: https://judge.yosupo.jp/submission/42413 """ N=len(F) r=pow(F[0],Mod-2,Mod) m=1 G=[r] while m>=1 if Q[0]==1: return P[0] else: return P[0]*pow(Q[0],Mod-2,Mod)%Mod def Nth_Term_of_Linearly_Recurrent_Sequence(A, C, N, offset=0): """ A[i]=C[0]*A[i-1]+C[1]*A[i-2]+...+C[d-1]*A[i-d] で表される数列 (A[i]) の第 N 項を求める. A=(A[0], ..., A[d-1]): 最初の d 項 C=(C[0], ..., C[d-1]): 線形漸化式 N: 求める項数 offset: ずらす項数 (初項が第 offset 項になる) """ assert len(A)==len(C) d=len(A) N-=offset if N<0: return 0 elif N