{- 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)]