{-
	Copyright (C) 2011 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 <http://www.gnu.org/licenses/>.
-}
{- |
 [@AUTHOR@]	Dr. Alistair Ward

 [@DESCRIPTION@]	Exports functions involving integral powers.
-}

module Factory.Math.Power(
-- * Functions
	square,
	squaresFrom,
	maybeSquareNumber,
	cube,
	cubeRoot,
	raiseModulo,
-- ** Predicates
	isPerfectPower
) where

import qualified	Data.Set

-- | Mainly for convenience.
{-# INLINE square #-}
square :: Num n => n -> n
square	= (^ (2 :: Int))

-- | Just for convenience.
cube :: Num n => n -> n
cube	= (^ (3 :: Int))

{- |
	* Iteratively generate sequential /squares/, from the specified initial value,
	based on the fact that @(x + 1)^2 = x^2 + 2 * x + 1@.

	* The initial value doesn't need to be either positive or integral.
-}
squaresFrom :: Num n => n -> [(n, n)]
squaresFrom from	= iterate (\(x, y) -> (x + 1, y + 2 * x + 1)) (from, square from)

-- | Just for convenience.
cubeRoot :: Double -> Double
cubeRoot	= (** recip 3)

{- |
	* Raise an arbitrary number to the specified positive integral power, using /modular/ arithmetic.

	* Implements exponentiation as a sequence of either /squares/ or multiplications by the base;
	<http://en.wikipedia.org/wiki/Exponentiation_by_squaring>.

	* <http://en.wikipedia.org/wiki/Modular_exponentiation>.
-}
raiseModulo :: (Integral i, Integral power)
	=> i	-- ^ Base.
	-> power
	-> i	-- ^ Modulus.
	-> i	-- ^ Result.
raiseModulo _ _ 0	= error "Factory.Math.Power.raiseModulo:\tzero modulus."
raiseModulo _ _ 1	= 0
raiseModulo _ 0 modulus	= 1 `mod` modulus
raiseModulo base power modulus
	| base < 0		= (`mod` modulus) . (if odd power then negate else id) $ raiseModulo (negate base) power modulus	--Recurse.
	| power < 0		= error $ "Factory.Math.Power.raiseModulo:\tnegative power; " ++ show power
	| first `elem` [0, 1]	= first
	| otherwise		= slave power
	where
		first	= base `mod` modulus

		slave 1	= first
		slave e	= (`mod` modulus) . (if r == 0 {-even-} then id else (* base)) . square $ slave q {-recurse-}	where
			(q, r)	= e `quotRem` 2

{- |
	* Returns @(Just . sqrt)@ if the specified integer is a /square number/ (AKA /perfect square/).

	* <http://en.wikipedia.org/wiki/Square_number>.

	* <http://mathworld.wolfram.com/SquareNumber.html>.

	* @(square . sqrt)@ is expensive, so the modulus of the operand is tested first, in an attempt to prove it isn't a /perfect square/.
	The set of tests, and the valid moduli within each test, are ordered to maximize the rate of failure-detection.
-}
maybeSquareNumber :: Integral i => i -> Maybe i
maybeSquareNumber i
--	| i < 0					= Nothing	--This function is performance-sensitive, but this test is neither strictly nor frequently required.
	| all (\(modulus, valid) -> mod i modulus `elem` valid) [
--							--Distribution of moduli amongst perfect squares	Cumulative failure-detection.
		(16,	[0,1,4,9]),			--All moduli are equally likely.			75%
		(9,	[0,1,4,7]),			--Zero occurs 33%, the others only 22%.			88%
		(17,	[1,2,4,8,9,13,15,16,0]),	--Zero only occurs 5.8%, the others 11.8%.		94%
-- These additional tests, aren't always cost-effective.
		(13,	[1,3,4,9,10,12,0]),		--Zero only occurs 7.7%, the others 15.4%.		97%
		(7,	[1,2,4,0]),			--Zero only occurs 14.3%, the others 28.6%.		98%
		(5,	[1,4,0])			--Zero only occurs 20%, the others 40%.			99%

--	] && fromIntegral iSqrt == sqrt'	= Just iSqrt	--CAVEAT: erroneously True for 187598574531033120, whereas 187598574531033121 is square.
	] && square iSqrt == i			= Just iSqrt
	| otherwise				= Nothing
	where
		sqrt' :: Double
		sqrt'	= sqrt $ fromIntegral i

		iSqrt	= round sqrt'

{- |
	* An integer @(> 1)@ which can be expressed as an integral power @(> 1)@ of a smaller /natural/ number.

	* CAVEAT: /zero/ and /one/ are normally excluded from this set.

	* <http://en.wikipedia.org/wiki/Perfect_power>.

	* <http://mathworld.wolfram.com/PerfectPower.html>.
-}
isPerfectPower :: Integral i => i -> Bool
isPerfectPower i
	| i < square 2	= False
	| otherwise	= i `Data.Set.member` foldr (
		\n set	-> if n `Data.Set.member` set
			then set
--			else Data.Set.union set . Data.Set.fromList . takeWhile (<= i) . iterate (* n) $ square n	--TODO: test relative speed.
			else foldr Data.Set.insert set . takeWhile (<= i) . iterate (* n) $ square n
	) Data.Set.empty [2 .. round $ sqrt (fromIntegral i :: Double)]