{-# LANGUAGE LambdaCase #-}
{-
	Copyright (C) 2011-2017 Dr. Alistair Ward

	This program is free software: you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation, either version 3 of the License, or
	(at your option) any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program.  If not, see <https://www.gnu.org/licenses/>.
-}
{- |
 [@AUTHOR@]	Dr. Alistair Ward

 [@DESCRIPTION@]

	* Implements several different prime-factorisation algorithms.

	* <https://www.tug.org/texinfohtml/coreutils.html#factor-invocation>.
-}

module Factory.Math.Implementations.PrimeFactorisation(
-- * Types
-- ** Data-types
	Algorithm(
--		DixonsMethod,
		FermatsMethod,
		TrialDivision
	)
-- * Functions
--	factoriseByDixonsMethod
--	factoriseByFermatsMethod
--	factoriseByTrialDivision
) where

import			Control.Arrow((&&&))
import qualified	Control.Arrow
import qualified	Control.DeepSeq
import qualified	Control.Parallel.Strategies
import qualified	Data.Default
import qualified	Data.Maybe
import qualified	Data.Numbers.Primes
import qualified	Factory.Data.Exponential	as Data.Exponential
import			Factory.Data.Exponential((<^))
import qualified	Factory.Data.PrimeFactors	as Data.PrimeFactors
import qualified	Factory.Math.PerfectPower	as Math.PerfectPower
import qualified	Factory.Math.Power		as Math.Power
import qualified	Factory.Math.PrimeFactorisation	as Math.PrimeFactorisation
import qualified	ToolShed.Data.Pair

-- | The algorithms by which prime-factorisation has been implemented.
data Algorithm
	= DixonsMethod	-- ^ <https://en.wikipedia.org/wiki/Dixon%27s_factorization_method>.
	| FermatsMethod	-- ^ <https://en.wikipedia.org/wiki/Fermat%27s_factorization_method>.
	| TrialDivision	-- ^ <https://en.wikipedia.org/wiki/Trial_division>.
	deriving (Algorithm -> Algorithm -> Bool
(Algorithm -> Algorithm -> Bool)
-> (Algorithm -> Algorithm -> Bool) -> Eq Algorithm
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Algorithm -> Algorithm -> Bool
$c/= :: Algorithm -> Algorithm -> Bool
== :: Algorithm -> Algorithm -> Bool
$c== :: Algorithm -> Algorithm -> Bool
Eq, ReadPrec [Algorithm]
ReadPrec Algorithm
Int -> ReadS Algorithm
ReadS [Algorithm]
(Int -> ReadS Algorithm)
-> ReadS [Algorithm]
-> ReadPrec Algorithm
-> ReadPrec [Algorithm]
-> Read Algorithm
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [Algorithm]
$creadListPrec :: ReadPrec [Algorithm]
readPrec :: ReadPrec Algorithm
$creadPrec :: ReadPrec Algorithm
readList :: ReadS [Algorithm]
$creadList :: ReadS [Algorithm]
readsPrec :: Int -> ReadS Algorithm
$creadsPrec :: Int -> ReadS Algorithm
Read, Int -> Algorithm -> ShowS
[Algorithm] -> ShowS
Algorithm -> String
(Int -> Algorithm -> ShowS)
-> (Algorithm -> String)
-> ([Algorithm] -> ShowS)
-> Show Algorithm
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Algorithm] -> ShowS
$cshowList :: [Algorithm] -> ShowS
show :: Algorithm -> String
$cshow :: Algorithm -> String
showsPrec :: Int -> Algorithm -> ShowS
$cshowsPrec :: Int -> Algorithm -> ShowS
Show)

instance Data.Default.Default Algorithm	where
	def :: Algorithm
def	= Algorithm
TrialDivision

instance Math.PrimeFactorisation.Algorithmic Algorithm	where
	primeFactors :: Algorithm -> base -> Factors base Int
primeFactors = \case
		Algorithm
DixonsMethod	-> base -> Factors base Int
forall base exponent. base -> Factors base exponent
factoriseByDixonsMethod
		Algorithm
FermatsMethod	-> Factors base Int -> Factors base Int
forall base exponent.
(Ord base, Num exponent, Ord exponent) =>
Factors base exponent -> Factors base exponent
Data.PrimeFactors.reduce (Factors base Int -> Factors base Int)
-> (base -> Factors base Int) -> base -> Factors base Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. base -> Factors base Int
forall base exponent.
(NFData base, NFData exponent, Integral base, Num exponent) =>
base -> Factors base exponent
factoriseByFermatsMethod
		Algorithm
TrialDivision	-> base -> Factors base Int
forall base exponent.
(Integral base, Num exponent) =>
base -> Factors base exponent
factoriseByTrialDivision

-- | <https://en.wikipedia.org/wiki/Dixon%27s_factorization_method>.
factoriseByDixonsMethod :: base -> Data.PrimeFactors.Factors base exponent
factoriseByDixonsMethod :: base -> Factors base exponent
factoriseByDixonsMethod	= base -> Factors base exponent
forall a. HasCallStack => a
undefined

{- |
	* <https://en.wikipedia.org/wiki/Fermat%27s_factorization_method>.

	* <https://mathworld.wolfram.com/FermatsFactorizationMethod.html>.

	* <https://en.wikipedia.org/wiki/Congruence_of_squares>.

	*	@i = f1 * f2@							Assume a non-trivial factorisation, ie. one in which both factors exceed one.
	=>	@i = (larger + smaller) * (larger - smaller)@			Represent the co-factors as a sum and difference.
	=>	@i = larger^2 - smaller^2@					Which has an integral solution if @i@ is neither /even/ nor a /perfect square/.
	=>	@sqrt (larger^2 - i) = smaller@					Search for /larger/, which results in an integral value for /smaller/.

	* Given that the smaller factor /f2/, can't be less than 3 (/i/ isn't /even/), then the larger /f1/, can't be greater than @(i `div` 3)@.
	So:	@(f2 >= 3) && (f1 <= i `div` 3)@				Two equations which can be used to solve for /larger/.
	=>	@(larger - smaller >= 3) && (larger + smaller <= i `div` 3)@	Add these to eliminate /smaller/.
	=>	@larger <= (i + 9) `div` 6@					The upper bound of the search-space.

	* This algorithm works best when there's a factor close to the /square-root/.
-}
factoriseByFermatsMethod :: (
	Control.DeepSeq.NFData	base,
	Control.DeepSeq.NFData	exponent,
	Integral		base,
	Num			exponent
 ) => base -> Data.PrimeFactors.Factors base exponent
factoriseByFermatsMethod :: base -> Factors base exponent
factoriseByFermatsMethod base
i
	| base
i base -> base -> Bool
forall a. Ord a => a -> a -> Bool
<= base
3				= [base -> Exponential base exponent
forall exponent base.
Num exponent =>
base -> Exponential base exponent
Data.Exponential.rightIdentity base
i]
	| base -> Bool
forall a. Integral a => a -> Bool
even base
i				= base -> Exponential base exponent
forall exponent base.
Num exponent =>
base -> Exponential base exponent
Data.Exponential.rightIdentity base
2 Exponential base exponent
-> Factors base exponent -> Factors base exponent
forall a. a -> [a] -> [a]
: base -> Factors base exponent
forall base exponent.
(NFData base, NFData exponent, Integral base, Num exponent) =>
base -> Factors base exponent
factoriseByFermatsMethod (base
i base -> base -> base
forall a. Integral a => a -> a -> a
`div` base
2) {-recurse-}
	| Maybe base -> Bool
forall a. Maybe a -> Bool
Data.Maybe.isJust Maybe base
maybeSquareNumber	= (Exponential base exponent -> exponent -> Exponential base exponent
forall exponent base.
Num exponent =>
Exponential base exponent -> exponent -> Exponential base exponent
<^ exponent
2) (Exponential base exponent -> Exponential base exponent)
-> Factors base exponent -> Factors base exponent
forall a b. (a -> b) -> [a] -> [b]
`map` base -> Factors base exponent
forall base exponent.
(NFData base, NFData exponent, Integral base, Num exponent) =>
base -> Factors base exponent
factoriseByFermatsMethod (Maybe base -> base
forall a. HasCallStack => Maybe a -> a
Data.Maybe.fromJust Maybe base
maybeSquareNumber) {-recurse-}
	| [(base, base)] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(base, base)]
factors				= [base -> Exponential base exponent
forall exponent base.
Num exponent =>
base -> Exponential base exponent
Data.Exponential.rightIdentity base
i]	-- Prime.
	| Bool
otherwise				= (Factors base exponent
 -> Factors base exponent -> Factors base exponent)
-> (Factors base exponent, Factors base exponent)
-> Factors base exponent
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Factors base exponent
-> Factors base exponent -> Factors base exponent
forall a. [a] -> [a] -> [a]
(++) ((Factors base exponent, Factors base exponent)
 -> Factors base exponent)
-> ((base, base) -> (Factors base exponent, Factors base exponent))
-> (base, base)
-> Factors base exponent
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Strategy (Factors base exponent, Factors base exponent)
-> (Factors base exponent, Factors base exponent)
-> (Factors base exponent, Factors base exponent)
forall a. Strategy a -> a -> a
Control.Parallel.Strategies.withStrategy (
		Strategy (Factors base exponent)
-> Strategy (Factors base exponent)
-> Strategy (Factors base exponent, Factors base exponent)
forall a b. Strategy a -> Strategy b -> Strategy (a, b)
Control.Parallel.Strategies.parTuple2 Strategy (Factors base exponent)
forall a. NFData a => Strategy a
Control.Parallel.Strategies.rdeepseq Strategy (Factors base exponent)
forall a. NFData a => Strategy a
Control.Parallel.Strategies.rdeepseq	-- CAVEAT: unproductive on the size of integers tested so far.
	) ((Factors base exponent, Factors base exponent)
 -> (Factors base exponent, Factors base exponent))
-> ((base, base) -> (Factors base exponent, Factors base exponent))
-> (base, base)
-> (Factors base exponent, Factors base exponent)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (base -> Factors base exponent)
-> (base, base) -> (Factors base exponent, Factors base exponent)
forall a b. (a -> b) -> (a, a) -> (b, b)
ToolShed.Data.Pair.mirror base -> Factors base exponent
forall base exponent.
(NFData base, NFData exponent, Integral base, Num exponent) =>
base -> Factors base exponent
factoriseByFermatsMethod ((base, base) -> Factors base exponent)
-> (base, base) -> Factors base exponent
forall a b. (a -> b) -> a -> b
$ [(base, base)] -> (base, base)
forall a. [a] -> a
head [(base, base)]
factors
	where
--		maybeSquareNumber :: Integral i => Maybe i
		maybeSquareNumber :: Maybe base
maybeSquareNumber	= base -> Maybe base
forall i. Integral i => i -> Maybe i
Math.PerfectPower.maybeSquareNumber base
i

--		factors :: Integral i => [i]
		factors :: [(base, base)]
factors	= ((base, Maybe base) -> (base, base))
-> [(base, Maybe base)] -> [(base, base)]
forall a b. (a -> b) -> [a] -> [b]
map (
			(
				(base -> base -> base) -> (base, base) -> base
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry base -> base -> base
forall a. Num a => a -> a -> a
(+) ((base, base) -> base)
-> ((base, base) -> base) -> (base, base) -> (base, base)
forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& (base -> base -> base) -> (base, base) -> base
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry (-)	-- Construct the co-factors as the sum and difference of /larger/ and /smaller/.
			) ((base, base) -> (base, base))
-> ((base, Maybe base) -> (base, base))
-> (base, Maybe base)
-> (base, base)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Maybe base -> base) -> (base, Maybe base) -> (base, base)
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (d, b) (d, c)
Control.Arrow.second Maybe base -> base
forall a. HasCallStack => Maybe a -> a
Data.Maybe.fromJust
		 ) ([(base, Maybe base)] -> [(base, base)])
-> (Double -> [(base, Maybe base)]) -> Double -> [(base, base)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((base, Maybe base) -> Bool)
-> [(base, Maybe base)] -> [(base, Maybe base)]
forall a. (a -> Bool) -> [a] -> [a]
filter (
			Maybe base -> Bool
forall a. Maybe a -> Bool
Data.Maybe.isJust (Maybe base -> Bool)
-> ((base, Maybe base) -> Maybe base) -> (base, Maybe base) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (base, Maybe base) -> Maybe base
forall a b. (a, b) -> b
snd	-- Search for a perfect square.
		 ) ([(base, Maybe base)] -> [(base, Maybe base)])
-> (Double -> [(base, Maybe base)])
-> Double
-> [(base, Maybe base)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((base, base) -> (base, Maybe base))
-> [(base, base)] -> [(base, Maybe base)]
forall a b. (a -> b) -> [a] -> [b]
map (
			(base -> Maybe base) -> (base, base) -> (base, Maybe base)
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (d, b) (d, c)
Control.Arrow.second ((base -> Maybe base) -> (base, base) -> (base, Maybe base))
-> (base -> Maybe base) -> (base, base) -> (base, Maybe base)
forall a b. (a -> b) -> a -> b
$ base -> Maybe base
forall i. Integral i => i -> Maybe i
Math.PerfectPower.maybeSquareNumber {-hotspot-} (base -> Maybe base) -> (base -> base) -> base -> Maybe base
forall b c a. (b -> c) -> (a -> b) -> a -> c
. base -> base -> base
forall a. Num a => a -> a -> a
subtract base
i	-- Associate the corresponding value of /smaller/.
		 ) ([(base, base)] -> [(base, Maybe base)])
-> (Double -> [(base, base)]) -> Double -> [(base, Maybe base)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((base, base) -> Bool) -> [(base, base)] -> [(base, base)]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (
			(base -> base -> Bool
forall a. Ord a => a -> a -> Bool
<= (base
i base -> base -> base
forall a. Num a => a -> a -> a
+ base
9) base -> base -> base
forall a. Integral a => a -> a -> a
`div` base
6) (base -> Bool) -> ((base, base) -> base) -> (base, base) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (base, base) -> base
forall a b. (a, b) -> a
fst	-- Terminate the search at the maximum value of /larger/.
		 ) ([(base, base)] -> [(base, base)])
-> (Double -> [(base, base)]) -> Double -> [(base, base)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. base -> [(base, base)]
forall n. (Enum n, Num n) => n -> [(n, n)]
Math.Power.squaresFrom {-hotspot-} (base -> [(base, base)])
-> (Double -> base) -> Double -> [(base, base)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> base
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Double -> [(base, base)]) -> Double -> [(base, base)]
forall a b. (a -> b) -> a -> b
$ Double -> Double
forall a. Floating a => a -> a
sqrt (base -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral base
i :: Double)	-- Start the search at the minimum value of /larger/.

{- |
	* Decomposes the specified integer, into a product of /prime/-factors,
	using <https://mathworld.wolfram.com/DirectSearchFactorization.html>, AKA <https://en.wikipedia.org/wiki/Trial_division>.

	* This works best when the factors are small.
-}
factoriseByTrialDivision :: (Integral base, Num exponent) => base -> Data.PrimeFactors.Factors base exponent
factoriseByTrialDivision :: base -> Factors base exponent
factoriseByTrialDivision	= [base] -> base -> Factors base exponent
forall t exponent.
(Num exponent, Integral t) =>
[t] -> t -> Factors t exponent
slave [base]
forall int. Integral int => [int]
Data.Numbers.Primes.primes where
	slave :: [t] -> t -> Factors t exponent
slave [t]
primes t
i
		| [t] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [t]
primeCandidates	= [t -> Exponential t exponent
forall exponent base.
Num exponent =>
base -> Exponential base exponent
Data.Exponential.rightIdentity t
i]
		| Bool
otherwise		= t -> Exponential t exponent
forall exponent base.
Num exponent =>
base -> Exponential base exponent
Data.Exponential.rightIdentity t
lowestPrimeFactor Exponential t exponent -> Factors t exponent -> Factors t exponent
forall base exponent.
(Ord base, Num exponent) =>
Exponential base exponent
-> Factors base exponent -> Factors base exponent
`Data.PrimeFactors.insert'` [t] -> t -> Factors t exponent
slave [t]
primeCandidates (t
i t -> t -> t
forall a. Integral a => a -> a -> a
`quot` t
lowestPrimeFactor)
		where
			primeCandidates :: [t]
primeCandidates	= (t -> Bool) -> [t] -> [t]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (
				(t -> t -> Bool
forall a. Eq a => a -> a -> Bool
/= t
0) (t -> Bool) -> (t -> t) -> t -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (t
i t -> t -> t
forall a. Integral a => a -> a -> a
`rem`)
			 ) ([t] -> [t]) -> [t] -> [t]
forall a b. (a -> b) -> a -> b
$ (t -> Bool) -> [t] -> [t]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (
				t -> t -> Bool
forall a. Ord a => a -> a -> Bool
<= t -> t
forall i. Integral i => i -> i
Math.PrimeFactorisation.maxBoundPrimeFactor t
i
			 ) [t]
primes

			lowestPrimeFactor :: t
lowestPrimeFactor	= [t] -> t
forall a. [a] -> a
head [t]
primeCandidates