{-# LANGUAGE BangPatterns #-} {-# OPTIONS_GHC -O2 #-} module Main where import Control.Monad (replicateM, replicateM_) import Data.ByteString.Char8 qualified as BS import Data.Maybe (fromJust) tuplify2 (x : y : _) = (x, y) tuplify2 _ = undefined -- Input functions with ByteString readInt = fst . fromJust . BS.readInteger readIntTuple = tuplify2 . map readInt . BS.words readIntList = map readInt . BS.words getInt = readInt <$> BS.getLine getIntList = readIntList <$> BS.getLine getIntNList n = map readIntList <$> replicateM (fromIntegral n) BS.getLine getIntMatrix = map readIntList . BS.lines <$> BS.getContents getIntTuple = readIntTuple <$> BS.getLine getIntNTuples n = map readIntTuple <$> replicateM (fromIntegral n) BS.getLine getIntTuples = map readIntTuple . BS.lines <$> BS.getContents -- | -- Max Weighted Floor (mwf) -- mwf(n,m,a,b,c,d) = max_{0 <= x < n} a*x + b*floor((c*x + d)/m) -- を返す非再帰(反復)実装。 -- -- 前提: -- - n > 0, m > 0 -- -- 計算量/メモリ: -- - 時間: O(log m)(ユークリッド互除法的再帰による構造縮約) -- - 追加メモリ: O(1) mwf :: Integer -> Integer -> Integer -> Integer -> Integer -> Integer -> Integer mwf n0 m0 a0 b0 c0 d0 = let !sum0 = 0 !max0 = b0 * (d0 `div` m0) -- x = 0 のとき in loop max0 sum0 n0 m0 a0 b0 c0 d0 where loop :: Integer -> Integer -> Integer -> Integer -> Integer -> Integer -> Integer -> Integer -> Integer loop !maxAcc !sumAcc !n !m !a !b !c !d = let (q1, c') = c `divMod` m !a' = a + b * q1 (q2, d') = d `divMod` m !sumAcc' = sumAcc + b * q2 !n1 = n - 1 !ymax = (c' * n1 + d') `div` m !maxAcc' = max maxAcc (max sumAcc' (sumAcc' + a' * n1 + b * ymax)) in if ymax == 0 || (a >= 0 && b >= 0) || (a <= 0 && b <= 0) then maxAcc' else if a' >= 0 then loop maxAcc' sumAcc' ymax c' b a' m (m - d' - 1) else loop maxAcc' (sumAcc' + a' + b) ymax c' b a' m (m - d' - 1) -- | -- max_{l <= x < r} a*x + b*floor((c*x + d)/m) を返す。 -- -- 既存の mwf(n,m,...)(0 <= x < n)を用いる。 -- 前提: l < r, m > 0 mwfLr :: Integer -> Integer -> Integer -> Integer -> Integer -> Integer -> Integer -> Integer mwfLr l r m a b c d | l >= r || m <= 0 = error "mwfLR: require l < r and m > 0" | otherwise = let n = r - l (q, d') = (c * l + d) `divMod` m in a * l + b * q + mwf n m a b c d' -- | -- Δ(D,A,B,x) = floor(x/D) - floor( (floor(x/A) * floor(A*B/D)) / B ) において、 -- u_min(D,A,B,K) = min { u >= 0 | Δ(D,A,B,u*D) > K } を半開区間二分探索 [0, A'BK+2) で求め、 -- x_min(D,A,B,K) = min { x >= 0 | Δ(D,A,B,x) > K } = u_min(D,A,B,K)*D を返す(解なしは -1)。 -- -- 前提: -- * D > 0, A > 0, B > 0, K >= 0(整数) computeXmin :: Integer -> Integer -> Integer -> Integer -> Integer computeXmin d a b k | d <= 0 || a <= 0 || b <= 0 || k < 0 = error "computeXmin: invalid input" | otherwise = let g = gcd d a d' = d `div` g a' = a `div` g (m', r') = (a' * b) `divMod` d' tDelta = b * k in if r' == 0 && d' * k + 1 >= a' then (-1) -- 解なし else -- 半開区間 [0, hi) で二分探索 let lo0 = 0 hi0 = a' * b * k + 2 loop !lo !hi | lo + 1 >= hi = lo | otherwise = let mid = (lo + hi) `div` 2 fMid = mwf mid a' b (-m') d' 0 in if fMid <= tDelta then loop mid hi else loop lo mid uMin = loop lo0 hi0 in d * uMin -- | -- 検算用 Δ(D,A,B,x)。 deltaVal :: Integer -> Integer -> Integer -> Integer -> Integer deltaVal d a b x = let p = x `div` d m = (a * b) `div` d q = ((x `div` a) * m) `div` b in p - q -- | -- 入出力(複数ケース) solve :: IO () solve = do t <- getInt replicateM_ (fromIntegral t :: Int) $ do [d, a, b, k] <- getIntList print (computeXmin d a b k) main :: IO () main = solve