|
| 1 | +-- https://github.com/minoki/my-atcoder-solutions |
| 2 | +{-# LANGUAGE BangPatterns #-} |
| 3 | +{-# LANGUAGE DataKinds #-} |
| 4 | +{-# LANGUAGE NoStarIsType #-} |
| 5 | +{-# LANGUAGE TypeFamilies #-} |
| 6 | +{-# LANGUAGE TypeOperators #-} |
| 7 | +import Control.Monad |
| 8 | +import Control.Monad.ST |
| 9 | +import qualified Data.ByteString.Char8 as BS |
| 10 | +import Data.Char (isSpace) |
| 11 | +import Data.Coerce |
| 12 | +import Data.Int (Int64) |
| 13 | +import Data.List (foldl', tails, unfoldr) |
| 14 | +import qualified Data.Vector.Unboxing as U |
| 15 | +import qualified Data.Vector.Unboxing.Mutable as UM |
| 16 | +import GHC.TypeNats (type (+), KnownNat, Nat, |
| 17 | + type (^), natVal) |
| 18 | + |
| 19 | +type Poly = U.Vector (IntMod (10^9 + 7)) |
| 20 | +type PolyM s = UM.MVector s (IntMod (10^9 + 7)) |
| 21 | + |
| 22 | +sum' :: KnownNat m => [IntMod m] -> IntMod m |
| 23 | +sum' = fromIntegral . foldl' (\x y -> x + unwrapN y) 0 |
| 24 | +{-# INLINE sum' #-} |
| 25 | + |
| 26 | +-- 多項式は |
| 27 | +-- U.fromList [a,b,c,...,z] = a + b * X + c * X^2 + ... + z * X^(k-1) |
| 28 | +-- により表す。 |
| 29 | + |
| 30 | +-- 多項式を X^k - X^(k-1) - ... - X - 1 で割った余りを返す。 |
| 31 | +reduceM :: Int -> PolyM s -> ST s (PolyM s) |
| 32 | +reduceM !k !v = loop (UM.length v) |
| 33 | + where loop !l | l <= k = return (UM.take l v) |
| 34 | + | otherwise = do b <- UM.read v (l - 1) |
| 35 | + forM_ [l - k - 1 .. l - 2] $ \i -> do |
| 36 | + UM.modify v (+ b) i |
| 37 | + loop (l - 1) |
| 38 | + |
| 39 | +-- 多項式の積を X^k - X^(k-1) - ... - X - 1 で割った余りを返す。 |
| 40 | +mulP :: Int -> Poly -> Poly -> Poly |
| 41 | +mulP !k !v !w = {- U.force $ -} U.create $ do |
| 42 | + let !vl = U.length v |
| 43 | + !wl = U.length w |
| 44 | + s <- UM.new (vl + wl - 1) |
| 45 | + forM_ [0 .. vl + wl - 2] $ \i -> do |
| 46 | + let !x = sum' [(v U.! (i-j)) * (w U.! j) | j <- [max 0 (i - vl + 1) .. min (wl - 1) i]] |
| 47 | + UM.write s i x |
| 48 | + reduceM k s |
| 49 | + |
| 50 | +-- 多項式に X をかけたものを X^k - X^(k-1) - ... - X - 1 で割った余りを返す。 |
| 51 | +mulByX :: Int -> Poly -> Poly |
| 52 | +mulByX !k !v |
| 53 | + | U.length v == k = let !v_k = v U.! (k-1) |
| 54 | + in U.generate k $ \i -> if i == 0 then |
| 55 | + v_k |
| 56 | + else |
| 57 | + v_k + (v U.! (i - 1)) |
| 58 | + | otherwise = U.cons 0 v |
| 59 | + |
| 60 | +-- X の(mod X^k - X^(k-1) - ... - X - 1 での)n 乗 |
| 61 | +powX :: Int -> Int -> Poly |
| 62 | +powX !k !n = doPowX n |
| 63 | + where |
| 64 | + doPowX 0 = U.fromList [1] -- 1 |
| 65 | + doPowX 1 = U.fromList [0,1] -- X |
| 66 | + doPowX i = case i `quotRem` 2 of |
| 67 | + (j,0) -> let !f = doPowX j -- X^j mod P |
| 68 | + in mulP k f f |
| 69 | + (j,_) -> let !f = doPowX j -- X^j mod P |
| 70 | + in mulByX k (mulP k f f) |
| 71 | + |
| 72 | +main :: IO () |
| 73 | +main = do |
| 74 | + [k,n] <- unfoldr (BS.readInt . BS.dropWhile isSpace) <$> BS.getLine |
| 75 | + -- 2 <= k <= 1000 |
| 76 | + -- 1 <= n <= 10^9 |
| 77 | + if n <= k then |
| 78 | + print 1 |
| 79 | + else do |
| 80 | + let f = powX k (n - k) -- X^(n-k) mod X^k - X^(k-1) - ... - X - 1 |
| 81 | + let seq = replicate k 1 ++ map (sum . take k) (tails seq) -- 数列 |
| 82 | + print $ sum $ zipWith (*) (U.toList f) (drop (k-1) seq) |
| 83 | + |
| 84 | +-- |
| 85 | +-- Modular Arithmetic |
| 86 | +-- |
| 87 | + |
| 88 | +newtype IntMod (m :: Nat) = IntMod { unwrapN :: Int64 } deriving (Eq) |
| 89 | + |
| 90 | +instance Show (IntMod m) where |
| 91 | + show (IntMod x) = show x |
| 92 | + |
| 93 | +instance KnownNat m => Num (IntMod m) where |
| 94 | + t@(IntMod x) + IntMod y |
| 95 | + | x + y >= modulus = IntMod (x + y - modulus) |
| 96 | + | otherwise = IntMod (x + y) |
| 97 | + where modulus = fromIntegral (natVal t) |
| 98 | + t@(IntMod x) - IntMod y |
| 99 | + | x >= y = IntMod (x - y) |
| 100 | + | otherwise = IntMod (x - y + modulus) |
| 101 | + where modulus = fromIntegral (natVal t) |
| 102 | + t@(IntMod x) * IntMod y = IntMod ((x * y) `rem` modulus) |
| 103 | + where modulus = fromIntegral (natVal t) |
| 104 | + fromInteger n = let result = IntMod (fromInteger (n `mod` fromIntegral modulus)) |
| 105 | + modulus = natVal result |
| 106 | + in result |
| 107 | + abs = undefined; signum = undefined |
| 108 | + {-# SPECIALIZE instance Num (IntMod 1000000007) #-} |
| 109 | + |
| 110 | +fromIntegral_Int64_IntMod :: KnownNat m => Int64 -> IntMod m |
| 111 | +fromIntegral_Int64_IntMod n = result |
| 112 | + where |
| 113 | + result | 0 <= n && n < modulus = IntMod n |
| 114 | + | otherwise = IntMod (n `mod` modulus) |
| 115 | + modulus = fromIntegral (natVal result) |
| 116 | + |
| 117 | +{-# RULES |
| 118 | +"fromIntegral/Int->IntMod" fromIntegral = fromIntegral_Int64_IntMod . (fromIntegral :: Int -> Int64) :: Int -> IntMod (10^9 + 7) |
| 119 | +"fromIntegral/Int64->IntMod" fromIntegral = fromIntegral_Int64_IntMod :: Int64 -> IntMod (10^9 + 7) |
| 120 | + #-} |
| 121 | + |
| 122 | +instance U.Unboxable (IntMod m) where |
| 123 | + type Rep (IntMod m) = Int64 |
0 commit comments