-- |
-- 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