#691684 (Python3) No.42 貯金箱の溜息

• 問題: No.42 貯金箱の溜息
• ユーザー: vwxyz
• 提出日時: 2021-08-21 07:40:23
• 言語: Python3 (3.12.2 + numpy 1.26.4 + scipy 1.12.0)
• 結果: AC
• 実行時間: 754 ms / 5,000 ms
• コード長: 20,417 bytes
• コンパイル時間: 162 ms
• コンパイル使用メモリ: 14,720 KB
• 実行使用メモリ: 13,952 KB
• 最終ジャッジ日時: 2024-10-14 16:19:18
• 合計ジャッジ時間: 3,864 ms
• ジャッジサーバーID （参考情報）: judge3 / judge4

入力	結果	実行時間; 実行使用メモリ
testcase_00	AC	754 ms; 13,824 KB
testcase_01	AC	736 ms; 13,952 KB
testcase_02	AC	741 ms; 13,952 KB

ソースコード

import bisect
import copy
import decimal
import fractions
import functools
import heapq
import itertools
import math
import random
import sys
from collections import Counter,deque,defaultdict
from functools import lru_cache,reduce
from heapq import heappush,heappop,heapify,heappushpop,_heappop_max,_heapify_max
def _heappush_max(heap,item):
    heap.append(item)
    heapq._siftdown_max(heap, 0, len(heap)-1)
def _heappushpop_max(heap, item):
    if heap and item < heap[0]:
        item, heap[0] = heap[0], item
        heapq._siftup_max(heap, 0)
    return item
from math import degrees, gcd as GCD
read=sys.stdin.read
readline=sys.stdin.readline
readlines=sys.stdin.readlines

def Extended_Euclid(n,m):
    stack=[]
    while m:
        stack.append((n,m))
        n,m=m,n%m
    if n>=0:
        x,y=1,0
    else:
        x,y=-1,0
    for i in range(len(stack)-1,-1,-1):
        n,m=stack[i]
        x,y=y,x-(n//m)*y
    return x,y

class MOD:
    def __init__(self,mod):
        self.mod=mod

    def Pow(self,a,n):
        a%=self.mod
        if n>=0:
            return pow(a,n,self.mod)
        else:
            assert math.gcd(a,self.mod)==1
            x=Extended_Euclid(a,self.mod)[0]
            return pow(x,-n,self.mod)

    def Build_Fact(self,N):
        assert N>=0
        self.factorial=[1]
        for i in range(1,N+1):
            self.factorial.append((self.factorial[-1]*i)%self.mod)
        self.factorial_inv=[None]*(N+1)
        self.factorial_inv[-1]=self.Pow(self.factorial[-1],-1)
        for i in range(N-1,-1,-1):
            self.factorial_inv[i]=(self.factorial_inv[i+1]*(i+1))%self.mod
        return self.factorial_inv

    def Fact(self,N):
        return self.factorial[N]

    def Fact_Inv(self,N):
        return self.factorial_inv[N]

    def Comb(self,N,K):
        if K<0 or K>N:
            return 0
        s=self.factorial[N]
        s=(s*self.factorial_inv[K])%self.mod
        s=(s*self.factorial_inv[N-K])%self.mod
        return s

class Lagrange_Interpolation:
    def __init__(self,X=False,Y=False,x0=None,xd=None,mod=0):
        self.degree=len(Y)-1
        self.mod=mod
        assert self.degree<self.mod
        if x0!=None and xd!=None:
            assert xd>0
            self.X=[(x0+i*xd)%self.mod for i in range(self.degree+1)]
            fact_inve=1
            for i in range(1,self.degree+1):
                fact_inve*=i*xd
                fact_inve%=self.mod
            fact_inve=MOD(self.mod).Pow(fact_inve,-1)
            self.coefficient=[y for y in Y]
            for i in range(self.degree-1,-1,-2):
                self.coefficient[i]*=-1
            for i in range(self.degree,-1,-1):
                self.coefficient[i]*=fact_inve
                self.coefficient[i]%=self.mod
                self.coefficient[self.degree-i]*=fact_inve
                self.coefficient[self.degree-i]%=self.mod
                fact_inve*=i*xd
                fact_inve%=self.mod        
        else:
            self.X=X
            assert len(self.X)==self.degree+1
            self.coefficient=[y for y in Y]
            for i in range(self.degree+1):
                for j in range(self.degree+1):
                    if i==j:
                        continue
                    self.coefficient[i]*=X[i]-X[j]
                    self.coefficient%=self.mod

    def __getitem__(self,N):
        N%=self.mod
        XX=[N-x for x in self.X]
        XX_left=[1]*(self.degree+2)
        for i in range(1,self.degree+2):
            XX_left[i]=XX_left[i-1]*XX[i-1]%self.mod
        XX_right=[1]*(self.degree+2)
        for i in range(self.degree,-1,-1):
            XX_right[i]=XX_right[i+1]*XX[i]%self.mod
        return sum(XX_left[i]*XX_right[i+1]*self.coefficient[i] for i in range(self.degree+1))%self.mod

def NTT(polynomial1,polynomial2):
    prim_root=3
    prim_root_inve=MOD(mod).Pow(prim_root,-1)
    def DFT(polynomial,inverse=False):
        dft=polynomial+[0]*((1<<n)-len(polynomial))
        if inverse:
            for bit in range(1,n+1):
                a=1<<bit-1
                x=pow(prim_root,mod-1>>bit,mod)
                U=[1]
                for _ in range(a):
                    U.append(U[-1]*x%mod)
                for i in range(1<<n-bit):
                    for j in range(a):
                        s=i*2*a+j
                        t=s+a
                        dft[s],dft[t]=(dft[s]+dft[t]*U[j])%mod,(dft[s]-dft[t]*U[j])%mod
        else:
            for bit in range(n,0,-1):
                a=1<<bit-1
                x=pow(prim_root_inve,mod-1>>bit,mod)
                U=[1]
                for _ in range(a):
                    U.append(U[-1]*x%mod)
                for i in range(1<<n-bit):
                    for j in range(a):
                        s=i*2*a+j
                        t=s+a
                        dft[s],dft[t]=(dft[s]+dft[t])%mod,U[j]*(dft[s]-dft[t])%mod
        return dft

    N=len(polynomial1)+len(polynomial2)-1
    n=(N-1).bit_length()
    ntt=[x*y%mod for x,y in zip(DFT(polynomial1),DFT(polynomial2))]
    ntt=DFT(ntt,inverse=True)
    x=pow((mod+1)//2,n)
    ntt=[ntt[i]*x%mod for i in range(N)]
    return ntt

def FFT(polynomial1,polynomial2,digit=10**5):
    def DFT(polynomial,n,inverse=False):
        N=len(polynomial)
        if inverse:
            primitive_root=[math.cos(-i*2*math.pi/(1<<n))+math.sin(-i*2*math.pi/(1<<n))*1j for i in range(1<<n)]
        else:
            primitive_root=[math.cos(i*2*math.pi/(1<<n))+math.sin(i*2*math.pi/(1<<n))*1j for i in range(1<<n)]
        dft=polynomial+[0]*((1<<n)-N)
        if inverse:
            for bit in range(1,n+1):
                a=1<<bit-1
                for i in range(1<<n-bit):
                    for j in range(a):
                        s=i*2*a+j
                        t=s+a
                        dft[s],dft[t]=dft[s]+dft[t]*primitive_root[j<<n-bit],dft[s]-dft[t]*primitive_root[j<<n-bit]
        else:
            for bit in range(n,0,-1):
                a=1<<bit-1
                for i in range(1<<n-bit):
                    for j in range(a):
                        s=i*2*a+j
                        t=s+a
                        dft[s],dft[t]=dft[s]+dft[t],primitive_root[j<<n-bit]*(dft[s]-dft[t])
        return dft
    
    def FFT_(polynomial1,polynomial2):
        N1=len(polynomial1)
        N2=len(polynomial2)
        N=N1+N2-1
        n=(N-1).bit_length()
        fft=[x*y for x,y in zip(DFT(polynomial1,n),DFT(polynomial2,n))]
        fft=DFT(fft,n,inverse=True)
        fft=[round((fft[i]/(1<<n)).real) for i in range(N)]
        return fft
    
    N1=len(polynomial1)
    N2=len(polynomial2)
    N=N1+N2-1
    polynomial11,polynomial12=[None]*N1,[None]*N1
    polynomial21,polynomial22=[None]*N2,[None]*N2
    for i in range(N1):
        polynomial11[i],polynomial12[i]=divmod(polynomial1[i],digit)
    for i in range(N2):
        polynomial21[i],polynomial22[i]=divmod(polynomial2[i],digit)
    polynomial=[0]*(N)
    a=digit**2-digit
    for i,x in enumerate(FFT_(polynomial11,polynomial21)):
        polynomial[i]+=x*a
    a=digit-1
    for i,x in enumerate(FFT_(polynomial12,polynomial22)):
        polynomial[i]-=x*a
    for i,x in enumerate(FFT_([x1+x2 for x1,x2 in zip(polynomial11,polynomial12)],[x1+x2 for x1,x2 in zip(polynomial21,polynomial22)])):
        polynomial[i]+=x*digit
    return polynomial

def Primitive_Root(p):
    if p==2:
        return 1
    if p==167772161:
        return 3
    if p==469762049:
        return 3
    if p==754974721:
        return 11
    if p==998244353:
        return 3
    if p==10**9+7:
        return 5
    divisors=[2]
    pp=(p-1)//2
    while pp%2==0:
        pp//=2
    for d in range(3,pp+1,2):
        if d**2>pp:
            break
        if pp%d==0:
            divisors.append(d)
            while pp%d==0:
                pp//=d
    if pp>1:
        divisors.append(pp)
    primitive_root=2
    while True:
        for d in divisors:
            if pow(primitive_root,(p-1)//d,p)==1:
                break
        else:
            return primitive_root
        primitive_root+=1

class Polynomial:
    def __init__(self,polynomial,max_degree=-1,eps=1e-12,mod=0):
        self.max_degree=max_degree
        if self.max_degree!=-1 and len(polynomial)>self.max_degree+1:
            self.polynomial=polynomial[:self.max_degree+1]
        else:
            self.polynomial=polynomial
        self.mod=mod
        self.eps=eps

    def __eq__(self,other):
        if type(other)!=Polynomial:
            return False
        if len(self.polynomial)!=len(other.polynomial):
            return False
        for i in range(len(self.polynomial)):
            if abs(self.polynomial[i]-other.polynomial[i])>self.eps:
                return False
        return True

    def __ne__(self,other):
        if type(other)!=Polynomial:
            return True
        if len(self.polynomial)!=len(other.polynomial):
            return True
        for i in range(len(self.polynomial)):
            if abs(self.polynomial[i]-other.polynomial[i])>self.eps:
                return True
        return False

    def __add__(self,other):
        assert type(other)==Polynomial
        summ=[0]*max(len(self.polynomial),len(other.polynomial))
        for i in range(len(self.polynomial)):
            summ[i]+=self.polynomial[i]
        for i in range(len(other.polynomial)):
            summ[i]+=other.polynomial[i]
        if self.mod:
            for i in range(len(summ)):
                summ[i]%=self.mod
        while summ and abs(summ[-1])<self.eps:
            summ.pop()
        summ=Polynomial(summ,max_degree=self.max_degree,eps=self.eps,mod=self.mod)
        return summ

    def __sub__(self,other):
        assert type(other)==Polynomial
        diff=[0]*max(len(self.polynomial),len(other.polynomial))
        for i in range(len(self.polynomial)):
            diff[i]+=self.polynomial[i]
        for i in range(len(other.polynomial)):
            diff[i]-=other.polynomial[i]
        if self.mod:
            for i in range(len(diff)):
                diff[i]%=self.mod
        while diff and abs(diff[-1])<self.eps:
            diff.pop()
        diff=Polynomial(diff,max_degree=self.max_degree,eps=self.eps,mod=self.mod)
        return diff

    def __mul__(self,other):
        if type(other)==Polynomial:
            if self.max_degree==-1:
                prod=[0]*(len(self.polynomial)+len(other.polynomial)-1)
                for i in range(len(self.polynomial)):
                    for j in range(len(other.polynomial)):
                        prod[i+j]+=self.polynomial[i]*other.polynomial[j]
            else:
                prod=[0]*min(len(self.polynomial)+len(other.polynomial)-1,self.max_degree+1)
                for i in range(len(self.polynomial)):
                    for j in range(min(len(other.polynomial),self.max_degree+1-i)):
                        prod[i+j]+=self.polynomial[i]*other.polynomial[j]
            if self.mod:
                for i in range(len(prod)):
                    prod[i]%=self.mod
            while prod and abs(prod[-1])<self.eps:
                prod.pop()
        else:
            if self.mod:
                prod=[x*other%self.mod for x in self.polynomial]
            else:
                prod=[x*other for x in self.polynomial]
            while prod and abs(prod[-1])<self.eps:
                prod.pop()
        prod=Polynomial(prod,max_degree=self.max_degree,eps=self.eps,mod=self.mod)
        return prod

    def __matmul__(self,other):
        assert type(other)==Polynomial
        if self.mod:
            prod=NTT(self.polynomial,other.polynomial)
        else:
            prod=FFT(self.polynomial,other.polynomial)
        if self.max_degree!=-1 and len(prod)>self.max_degree+1:
            prod=prod[:self.max_degree+1]
            while prod and abs(prod[-1])<self.eps:
                prod.pop()
        prod=Polynomial(prod,max_degree=self.max_degree,eps=self.eps,mod=self.mod)
        return prod

    def __truediv__(self,other):
        if type(other)==Polynomial:
            assert other.polynomial
            for n in range(len(other.polynomial)):
                if self.eps<abs(other.polynomial[n]):
                    break
            assert len(self.polynomial)>n
            for i in range(n):
                assert abs(self.polynomial[i])<self.eps
            self_polynomial=self.polynomial[n:]
            other_polynomial=other.polynomial[n:]
            if self.mod:
                inve=MOD(self.mod).Pow(other_polynomial[0],-1)
            else:
                inve=1/other_polynomial[0]
            quot=[]
            for i in range(len(self_polynomial)-len(other_polynomial)+1):
                if self.mod:
                    quot.append(self_polynomial[i]*inve%self.mod)
                else:
                    quot.append(self_polynomial[i]*inve)
                for j in range(len(other_polynomial)):
                    self_polynomial[i+j]-=other_polynomial[j]*quot[-1]
                    if self.mod:
                        self_polynomial[i+j]%=self.mod
            for i in range(len(self_polynomial)-len(other_polynomial)+1,len(self_polynomial)):
                if self.eps<abs(self_polynomial[i]):
                    assert self.max_degree!=-1
                    self_polynomial=self_polynomial[-len(other_polynomial)+1:]
                    while len(quot)<=self.max_degree:
                        self_polynomial.append(0)
                        if self.mod:
                            quot.append(self_polynomial[0]*inve%self.mod)
                            self_polynomial=[(self_polynomial[i]-other_polynomial[i]*quot[-1])%self.mod for i in range(1,len(self_polynomial))]
                        else:
                            quot.append(self_polynomial[0]*inve)
                            self_polynomial=[(self_polynomial[i]-other_polynomial[i]*quot[-1]) for i in range(1,len(self_polynomial))]
                    break
            quot=Polynomial(quot,max_degree=self.max_degree,eps=self.eps,mod=self.mod)
        else:
            assert self.eps<abs(other)
            if self.mod:
                inve=MOD(self.mod).Pow(other,-1)
                quot=Polynomial([x*inve%self.mod for x in self.polynomial],max_degree=self.max_degree,eps=self.eps,mod=self.mod)
            else:
                quot=Polynomial([x/other for x in self.polynomial],max_degree=self.max_degree,eps=self.eps,mod=self.mod)
        return quot

    def __floordiv__(self,other):
        assert type(other)==Polynomial
        quot=[0]*(len(self.polynomial)-len(other.polynomial)+1)
        rema=[x for x in self.polynomial]
        if self.mod:
            inve=MOD(self.mod).Pow(other.polynomial[-1],-1)
            for i in range(len(self.polynomial)-len(other.polynomial),-1,-1):
                quot[i]=rema[i+len(other.polynomial)-1]*inve%self.mod
                for j in range(len(other.polynomial)):
                    rema[i+j]-=quot[i]*other.polynomial[j]
                    rema[i+j]%=self.mod
        else:
            inve=1/other.polynomial[-1]
            for i in range(len(self.polynomial)-len(other.polynomial),-1,-1):
                quot[i]=rema[i+len(other.polynomial)-1]*inve
                for j in range(len(other.polynomial)):
                    rema[i+j]-=quot[i]*other.polynomial[j]
        quot=Polynomial(quot,max_degree=self.max_degree,eps=self.eps,mod=self.mod)
        return quot

    def __mod__(self,other):
        assert type(other)==Polynomial
        quot=[0]*(len(self.polynomial)-len(other.polynomial)+1)
        rema=[x for x in self.polynomial]
        if self.mod:
            inve=MOD(self.mod).Pow(other.polynomial[-1],-1)
            for i in range(len(self.polynomial)-len(other.polynomial),-1,-1):
                quot[i]=rema[i+len(other.polynomial)-1]*inve%self.mod
                for j in range(len(other.polynomial)):
                    rema[i+j]-=quot[i]*other.polynomial[j]
                    rema[i+j]%=self.mod
                rema.pop()
        else:
            inve=1/other.polynomial[-1]
            for i in range(len(self.polynomial)-len(other.polynomial),-1,-1):
                quot[i]=rema[i+len(other.polynomial)-1]*inve
                for j in range(len(other.polynomial)):
                    rema[i+j]-=quot[i]*other.polynomial[j]
                rema.pop()
        rema=Polynomial(rema,max_degree=self.max_degree,eps=self.eps,mod=self.mod)
        return rema

    def __divmod__(self,other):
        assert type(other)==Polynomial
        quot=[0]*(len(self.polynomial)-len(other.polynomial)+1)
        rema=[x for x in self.polynomial]
        if self.mod:
            inve=MOD(self.mod).Pow(other.polynomial[-1],-1)
            for i in range(len(self.polynomial)-len(other.polynomial),-1,-1):
                quot[i]=rema[i+len(other.polynomial)-1]*inve%self.mod
                for j in range(len(other.polynomial)):
                    rema[i+j]-=quot[i]*other.polynomial[j]
                    rema[i+j]%=self.mod
                rema.pop()
        else:
            inve=1/other.polynomial[-1]
            for i in range(len(self.polynomial)-len(other.polynomial),-1,-1):
                quot[i]=rema[i+len(other.polynomial)-1]*inve
                for j in range(len(other.polynomial)):
                    rema[i+j]-=quot[i]*other.polynomial[j]
                rema.pop()
        quot=Polynomial(quot,max_degree=self.max_degree,eps=self.eps,mod=self.mod)
        rema=Polynomial(rema,max_degree=self.max_degree,eps=self.eps,mod=self.mod)
        return quot,rema

    def __neg__(self):
        if self.mod:
            nega=Polynomial([(-x)%self.mod for x in self.polynomial],max_degree=self.max_degree,eps=self.eps,mod=self.mod)
        else:
            nega=Polynomial([-x for x in self.polynomial],max_degree=self.max_degree,eps=self.eps,mod=self.mod)
        return nega

    def __pos__(self):
        posi=Polynomial([x for x in self.polynomial],max_degree=self.max_degree,eps=self.eps,mod=self.mod)
        return posi

    def __bool__(self):
        return self.polynomial

    def __getitem__(self,n):
        if n<=len(self.polynomial)-1:
            return self.polynomial[n]
        else:
            return 0

    def __setitem__(self,n,x):
        if self.mod:
            x%=self.mod
        if self.max_degree==-1 or n<=self.max_degree:
            if n<=len(self.polynomial)-1:
                self.polynomial[n]=x
            elif self.eps<abs(x):
                self.polynomial+=[0]*(n-len(self.polynomial))+[x]

    def __call__(self,x):
        retu=0
        pow_x=1
        for i in range(len(self.polynomial)):
            retu+=pow_x*self.polynomial[i]
            pow_x*=x
            if self.mod:
                retu%=self.mod
                pow_x%=self.mod
        return retu

    def __str__(self):
        return "["+", ".join(map(str,self.polynomial))+"]"

def Bostan_Mori(poly_deno,poly_nume,N,mod=0,fft=False,ntt=False):
    if type(poly_deno)==Polynomial:
        poly_deno=poly_deno.polynomial
    if type(poly_nume)==Polynomial:
        poly_nume=poly_nume.polynomial
    if ntt:
        convolve=NTT
    elif fft:
        convolve=FFT
    else:
        def convolve(poly_deno,poly_nume):
            conv=[0]*(len(poly_deno)+len(poly_nume)-1)
            for i in range(len(poly_deno)):
                for j in range(len(poly_nume)):
                    conv[i+j]+=poly_deno[i]*poly_nume[j]
            if mod:
                for i in range(len(conv)):
                    conv[i]%=mod
            return conv
    while N:
        poly_nume_=[-x if i%2 else x for i,x in enumerate(poly_nume)]
        if N%2:
            poly_deno=convolve(poly_deno,poly_nume_)[1::2]
        else:
            poly_deno=convolve(poly_deno,poly_nume_)[::2]
        poly_nume=convolve(poly_nume,poly_nume_)[::2]
        if fft and mod:
            for i in range(len(poly_deno)):
                poly_deno[i]%=mod
            for i in range(len(poly_nume)):
                poly_nume[i]%=mod
        N//=2
    return poly_deno[0]

mod=10**9+9
P=Polynomial([1],max_degree=3000-1,mod=mod)
for money in (1,5,10,50,100,500):
    lst=[0]*(money+1)
    lst[0]=1
    lst[money]=-1
    P/=Polynomial(lst)
T=int(readline())
for _ in range(T):
    M=int(readline())
    r=M%500
    Y=[P[i] for i in range(r,3000,500)]
    LI=Lagrange_Interpolation(Y=Y,x0=r,xd=500,mod=mod)
    ans=LI[M]
    print(ans)