diff --git a/app/Main.hs b/app/Main.hs index d75498a..3bf7683 100644 --- a/app/Main.hs +++ b/app/Main.hs @@ -1,33 +1,46 @@ module Main (main) where import Control.Exception -import Data.Vector qualified as V import Poly import System.Random import System.TimeIt +import Control.Monad main :: IO () main = do do - let f :: Poly Int = Poly (V.fromList [1, 2, 3]) - let g :: Poly Int = Poly (V.fromList [4, 5]) + let f :: Poly Int = makePoly [1, 2, 3] + let g :: Poly Int = makePoly [4, 5] putStrLn $ "f: " <> show f <> ", g: " <> show g putStrLn $ "f + g: " <> show (f + g) putStrLn $ "Naive f * g: " <> show (f * g) putStrLn $ "Karatsuba f * g: " <> show (normalize $ karatsubaMult f g) + putStrLn "" + experimentFor 250 experimentFor 500 experimentFor 1000 + karatsubaFor 3000 where experimentFor n = do setStdGen $ mkStdGen 10 - let randomPoly size = Poly <$> V.replicateM size (randomRIO (-100, 100)) + let randomPoly size = makePoly <$> replicateM size (randomRIO (-100, 100)) putStrLn $ "Size " <> show n f :: Poly Int <- randomPoly n g :: Poly Int <- randomPoly n putStrLn "naive:" _ <- timeIt $ evaluate (f * g) putStrLn "Karatsuba:" - _ <- timeIt $ evaluate (f * g) + _ <- timeIt $ evaluate (karatsubaMult f g) + putStrLn "Finished" + + karatsubaFor n = do + setStdGen $ mkStdGen 10 + let randomPoly size = makePoly <$> replicateM size (randomRIO (-100, 100)) + putStrLn $ "Size " <> show n + f :: Poly Int <- randomPoly n + g :: Poly Int <- randomPoly n + putStrLn "Karatsuba:" + _ <- timeIt $ evaluate (karatsubaMult f g) putStrLn "Finished" diff --git a/src/Poly.hs b/src/Poly.hs index cca0eee..4214554 100644 --- a/src/Poly.hs +++ b/src/Poly.hs @@ -2,64 +2,67 @@ module Poly where import Data.List import Data.Maybe -import Data.Vector qualified as V +import Data.Vector.Unboxed qualified as V -- Zip two vectors while padding 0s on the shorter vector. -vecZipPad0With :: (Num a) => (a -> a -> a) -> V.Vector a -> V.Vector a -> V.Vector a +vecZipPad0With :: (V.Unbox a, Num a) => (a -> a -> a) -> V.Vector a -> V.Vector a -> V.Vector a vecZipPad0With f xs ys = V.generate (max (V.length xs) (V.length ys)) $ \i -> fromMaybe 0 (xs V.!? i) `f` fromMaybe 0 (ys V.!? i) -- | Polynomial type. -- --- >>> Poly (V.fromList [1 .. 5]) +-- >>> Poly (V.fromList [1 :: Int .. 5]) -- 1 X^0 + 2 X^1 + 3 X^2 + 4 X^3 + 5 X^4 --- >>> Poly (V.fromList [1, 2]) * Poly (V.fromList [3, 4, 5]) +-- >>> Poly (V.fromList [1 :: Int, 2]) * Poly (V.fromList [3 :: Int, 4, 5]) -- 3 X^0 + 10 X^1 + 13 X^2 + 10 X^3 --- >>> Poly (V.fromList [1, 2]) * Poly (V.fromList []) +-- >>> Poly (V.fromList [1 :: Int, 2]) * Poly (V.fromList []) -- 0 X^0 newtype Poly a = Poly (V.Vector a) deriving (Eq) +makePoly :: (V.Unbox a) => [a] -> Poly a +makePoly = Poly . V.fromList + -- | Degree, assuming top term is nonzero -degree :: Poly a -> Int -degree (Poly f) = length f - 1 +degree :: (V.Unbox a) => Poly a -> Int +degree (Poly f) = V.length f - 1 -- | Shift up polynomial by X^n -shiftUp :: (Num a) => Int -> Poly a -> Poly a +shiftUp :: (V.Unbox a, Num a) => Int -> Poly a -> Poly a shiftUp n (Poly f) = Poly $ V.replicate n 0 <> f -- | Shift down polynomial by X^n -shiftDown :: Int -> Poly a -> Poly a +shiftDown :: (V.Unbox a) => Int -> Poly a -> Poly a shiftDown n (Poly f) = Poly $ V.drop n f -- | Remainder under X^n -remXn :: Int -> Poly a -> Poly a +remXn :: (V.Unbox a) => Int -> Poly a -> Poly a remXn n (Poly f) = Poly $ V.take n f -- | Normalize polynomial, removing leading 0s --- --- >>> normalize $ Poly (V.fromList [1, 0, 0]) +-- +-- >>> normalize $ Poly (V.fromList [1 :: Int, 0, 0]) -- 1 X^0 --- --- >>> normalize $ Poly (V.fromList [1, 2, 3, 0]) +-- +-- >>> normalize $ Poly (V.fromList [1 :: Int, 2, 3, 0]) -- 1 X^0 + 2 X^1 + 3 X^2 -normalize :: (Eq a, Num a) => Poly a -> Poly a +normalize :: (Eq a, Num a, V.Unbox a) => Poly a -> Poly a normalize (Poly f) = Poly remain where (_, remain) = V.spanR (== 0) f -- | This Num instance implements the classical multiplication. -instance (Num a) => Num (Poly a) where +instance (Num a, V.Unbox a) => Num (Poly a) where (+) :: Poly a -> Poly a -> Poly a Poly f + Poly g = Poly $ vecZipPad0With (+) f g (-) :: Poly a -> Poly a -> Poly a Poly f - Poly g = Poly $ vecZipPad0With (-) f g (*) :: Poly a -> Poly a -> Poly a - Poly f * Poly g = sum (Poly <$> mults) + Poly f * Poly g = sum (map Poly mults) where - mults = V.imap (\i fi -> V.map (fi *) (V.replicate i 0 <> g)) f + mults = zipWith (\i fi -> V.map (fi *) (V.replicate i 0 <> g)) [0 ..] (V.toList f) negate :: Poly a -> Poly a negate (Poly f) = Poly $ V.map negate f abs :: Poly a -> Poly a @@ -69,11 +72,10 @@ instance (Num a) => Num (Poly a) where fromInteger :: Integer -> Poly a fromInteger = Poly . V.singleton . fromInteger -instance (Show a) => Show (Poly a) where - show :: (Show a) => Poly a -> String - show (Poly p) = intercalate " + " . V.toList $ V.imap (\i coeff -> show coeff <> " X^" <> show i) p +instance (V.Unbox a, Show a) => Show (Poly a) where + show (Poly p) = intercalate " + " $ zipWith (\i coeff -> show coeff <> " X^" <> show i) [0 :: Int ..] (V.toList p) -karatsubaMult :: (Num a) => Poly a -> Poly a -> Poly a +karatsubaMult :: (Num a, V.Unbox a) => Poly a -> Poly a -> Poly a karatsubaMult a b = atLog degBound a b where degBound = fromJust $ find (> max (degree a) (degree b)) [2 ^ i | i <- [0 :: Int ..]]