    -- |
    -- Module:      Math.NumberTheory.SmoothNumbers
    -- Copyright:   (c) 2018 Frederick Schneider, 2018-2019 Andrew Lelechenko
    -- Licence:     MIT
    -- Maintainer:  Frederick Schneider <frederick.schneider2011@gmail.com>
    --
    -- A <https://en.wikipedia.org/wiki/Smooth_number smooth number>
    -- is an number, which can be represented as a product of powers of elements
    -- from a given set (smooth basis). E. g., 48 = 3 * 4 * 4 is smooth
    -- over a set {3, 4}, and 24 is not.
    --

    {-# LANGUAGE FlexibleContexts    #-}
    {-# LANGUAGE ScopedTypeVariables #-}

    module Math.NumberTheory.SmoothNumbers
      ( SmoothBasis
      , unSmoothBasis
      , fromList
      , isSmooth
      , smoothOver
      , smoothOver'
      ) where

    import Prelude hiding (div, mod, gcd, (+))
    import Data.Euclidean
    import Data.List (nub)
    import Data.Maybe
    import Data.Semiring

    -- | An abstract representation of a smooth basis.
    newtype SmoothBasis a = SmoothBasis
      { forall a. SmoothBasis a -> [a]
    unSmoothBasis :: [a]
      -- ^ Unwrap elements of a smooth basis.
      }
      deriving (Int -> SmoothBasis a -> ShowS
    [SmoothBasis a] -> ShowS
    SmoothBasis a -> String
    (Int -> SmoothBasis a -> ShowS)
    -> (SmoothBasis a -> String)
    -> ([SmoothBasis a] -> ShowS)
    -> Show (SmoothBasis a)
    forall a. Show a => Int -> SmoothBasis a -> ShowS
    forall a. Show a => [SmoothBasis a] -> ShowS
    forall a. Show a => SmoothBasis a -> String
    forall a.
    (Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
    $cshowsPrec :: forall a. Show a => Int -> SmoothBasis a -> ShowS
    showsPrec :: Int -> SmoothBasis a -> ShowS
    $cshow :: forall a. Show a => SmoothBasis a -> String
    show :: SmoothBasis a -> String
    $cshowList :: forall a. Show a => [SmoothBasis a] -> ShowS
    showList :: [SmoothBasis a] -> ShowS
    Show)

    -- | Build a 'SmoothBasis' from a list of numbers,
    -- sanitizing it from duplicates, zeros and units.
    --
    -- >>> fromList [2, 3]
    -- SmoothBasis {unSmoothBasis = [2,3]}
    -- >>> fromList [2, 2]
    -- SmoothBasis {unSmoothBasis = [2]}
    -- >>> fromList [1, 3]
    -- SmoothBasis {unSmoothBasis = [3]}
    fromList :: (Eq a, GcdDomain a) => [a] -> SmoothBasis a
    fromList :: forall a. (Eq a, GcdDomain a) => [a] -> SmoothBasis a
    fromList
      = [a] -> SmoothBasis a
    forall a. [a] -> SmoothBasis a
    SmoothBasis
      ([a] -> SmoothBasis a) -> ([a] -> [a]) -> [a] -> SmoothBasis a
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (a -> Bool) -> [a] -> [a]
    forall a. (a -> Bool) -> [a] -> [a]
    filter (\a
    x -> Bool -> Bool
    not (a -> Bool
    forall a. (Eq a, Semiring a) => a -> Bool
    isZero a
    x) Bool -> Bool -> Bool
    && Maybe a -> Bool
    forall a. Maybe a -> Bool
    isNothing (a
    forall a. Semiring a => a
    one a -> a -> Maybe a
    forall a. GcdDomain a => a -> a -> Maybe a
    `divide` a
    x))
      ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [a]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . [a] -> [a]
    forall a. Eq a => [a] -> [a]
    nub

    -- | A generalization of 'smoothOver',
    -- suitable for non-'Integral' domains.
    -- The first argument must be an appropriate
    -- <https://en.wikipedia.org/wiki/Ideal_norm Ideal norm> function,
    -- like 'Math.NumberTheory.Quadratic.GaussianIntegers.norm'
    -- or 'Math.NumberTheory.Quadratic.EisensteinIntegers.norm'.
    --
    -- This routine is more efficient than filtering with 'isSmooth'.
    --
    -- >>> import Math.NumberTheory.Quadratic.GaussianIntegers
    -- >>> take 10 (smoothOver' norm (fromList [1+ι,3]))
    -- [1,1+ι,2,2+2*ι,3,4,3+3*ι,4+4*ι,6,8]
    smoothOver'
      :: (Eq a, Num a, Ord b)
      => (a -> b) -- ^ <https://en.wikipedia.org/wiki/Ideal_norm Ideal norm>
      -> SmoothBasis a
      -> [a]
    smoothOver' :: forall a b.
    (Eq a, Num a, Ord b) =>
    (a -> b) -> SmoothBasis a -> [a]
    smoothOver' a -> b
    norm (SmoothBasis [a]
    pl) =
      (a -> [a] -> [a]) -> [a] -> [a] -> [a]
    forall a b. (a -> b -> b) -> b -> [a] -> b
    forall (t :: * -> *) a b.
    Foldable t =>
    (a -> b -> b) -> b -> t a -> b
    foldr
      (\a
    p [a]
    l -> ([a] -> [a] -> [a]) -> [a] -> [[a]] -> [a]
    forall a b. (a -> b -> b) -> b -> [a] -> b
    forall (t :: * -> *) a b.
    Foldable t =>
    (a -> b -> b) -> b -> t a -> b
    foldr [a] -> [a] -> [a]
    skipHead [] ([[a]] -> [a]) -> [[a]] -> [a]
    forall a b. (a -> b) -> a -> b
    $ ([a] -> [a]) -> [a] -> [[a]]
    forall a. (a -> a) -> a -> [a]
    iterate ((a -> a) -> [a] -> [a]
    forall a b. (a -> b) -> [a] -> [b]
    map (a -> a
    forall a. Num a => a -> a
    abs (a -> a) -> (a -> a) -> a -> a
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (a -> a -> a
    forall a. Num a => a -> a -> a
    Prelude.* a
    p))) [a]
    l)
      [a
    1]
      [a]
    pl
      where
        skipHead :: [a] -> [a] -> [a]
    skipHead []      [a]
    b = [a]
    b
        skipHead (a
    h : [a]
    t) [a]
    b = a
    h a -> [a] -> [a]
    forall a. a -> [a] -> [a]
    : [a] -> [a] -> [a]
    merge [a]
    t [a]
    b

        merge :: [a] -> [a] -> [a]
    merge [a]
    a [] = [a]
    a
        merge [] [a]
    b = [a]
    b
        merge a :: [a]
    a@(a
    ah : [a]
    at) b :: [a]
    b@(a
    bh : [a]
    bt)
          | a -> b
    norm a
    bh b -> b -> Bool
    forall a. Ord a => a -> a -> Bool
    < a -> b
    norm a
    ah = a
    bh a -> [a] -> [a]
    forall a. a -> [a] -> [a]
    : [a] -> [a] -> [a]
    merge [a]
    a [a]
    bt
          | a
    ah a -> a -> Bool
    forall a. Eq a => a -> a -> Bool
    == a
    bh          = a
    ah a -> [a] -> [a]
    forall a. a -> [a] -> [a]
    : [a] -> [a] -> [a]
    merge [a]
    at [a]
    bt
          | Bool
    otherwise         = a
    ah a -> [a] -> [a]
    forall a. a -> [a] -> [a]
    : [a] -> [a] -> [a]
    merge [a]
    at [a]
    b

    -- | Generate an infinite ascending list of
    -- <https://en.wikipedia.org/wiki/Smooth_number smooth numbers>
    -- over a given smooth basis.
    --
    -- This routine is more efficient than filtering with 'isSmooth'.
    --
    -- >>> take 10 (smoothOver (fromList [2, 5]))
    -- [1,2,4,5,8,10,16,20,25,32]
    smoothOver :: Integral a => SmoothBasis a -> [a]
    smoothOver :: forall a. Integral a => SmoothBasis a -> [a]
    smoothOver = (a -> a) -> SmoothBasis a -> [a]
    forall a b.
    (Eq a, Num a, Ord b) =>
    (a -> b) -> SmoothBasis a -> [a]
    smoothOver' a -> a
    forall a. Num a => a -> a
    abs

    -- | Check that a given number is smooth under a given 'SmoothBasis'.
    --
    -- >>> isSmooth (fromList [2,3]) 12
    -- True
    -- >>> isSmooth (fromList [2,3]) 15
    -- False
    isSmooth :: (Eq a, GcdDomain a) => SmoothBasis a -> a -> Bool
    isSmooth :: forall a. (Eq a, GcdDomain a) => SmoothBasis a -> a -> Bool
    isSmooth SmoothBasis a
    prs a
    x = Bool -> Bool
    not (a -> Bool
    forall a. (Eq a, Semiring a) => a -> Bool
    isZero a
    x) Bool -> Bool -> Bool
    && [a] -> a -> Bool
    forall a. (Eq a, GcdDomain a) => [a] -> a -> Bool
    go (SmoothBasis a -> [a]
    forall a. SmoothBasis a -> [a]
    unSmoothBasis SmoothBasis a
    prs) a
    x
      where
        go :: (Eq a, GcdDomain a) => [a] -> a -> Bool
        go :: forall a. (Eq a, GcdDomain a) => [a] -> a -> Bool
    go [] a
    n = Maybe a -> Bool
    forall a. Maybe a -> Bool
    isJust (a
    forall a. Semiring a => a
    one a -> a -> Maybe a
    forall a. GcdDomain a => a -> a -> Maybe a
    `divide` a
    n)
        go pps :: [a]
    pps@(a
    p:[a]
    ps) a
    n = case a
    n a -> a -> Maybe a
    forall a. GcdDomain a => a -> a -> Maybe a
    `divide` a
    p of
          Maybe a
    Nothing -> [a] -> a -> Bool
    forall a. (Eq a, GcdDomain a) => [a] -> a -> Bool
    go [a]
    ps a
    n
          Just a
    q  -> [a] -> a -> Bool
    forall a. (Eq a, GcdDomain a) => [a] -> a -> Bool
    go [a]
    pps a
    q Bool -> Bool -> Bool
    || [a] -> a -> Bool
    forall a. (Eq a, GcdDomain a) => [a] -> a -> Bool
    go [a]
    ps a
    n