module Math.NumberTheory.Powers.Fourth
( integerFourthRoot
, integerFourthRoot'
, exactFourthRoot
, isFourthPower
, isFourthPower'
, isPossibleFourthPower
) where
#include "MachDeps.h"
import GHC.Base
import GHC.Integer
import GHC.Integer.GMP.Internals
import Data.Array.Unboxed
import Data.Array.ST
import Data.Array.Base (unsafeAt, unsafeWrite)
import Data.Bits
import Data.Word
import Math.NumberTheory.Logarithms.Internal (integerLog2#)
integerFourthRoot :: Integral a => a -> a
integerFourthRoot n
| n < 0 = error "integerFourthRoot: negative argument"
| otherwise = integerFourthRoot' n
integerFourthRoot' :: Integral a => a -> a
integerFourthRoot' 0 = 0
integerFourthRoot' n = newton4 n (approxBiSqrt n)
exactFourthRoot :: Integral a => a -> Maybe a
exactFourthRoot 0 = Just 0
exactFourthRoot n
| n < 0 = Nothing
| isPossibleFourthPower n && r2*r2 == n = Just r
| otherwise = Nothing
where
r = integerFourthRoot' n
r2 = r*r
isFourthPower :: Integral a => a -> Bool
isFourthPower 0 = True
isFourthPower n = n > 0 && isFourthPower' n
isFourthPower' :: Integral a => a -> Bool
isFourthPower' n = isPossibleFourthPower n && r2*r2 == n
where
r = integerFourthRoot' n
r2 = r*r
isPossibleFourthPower :: Integral a => a -> Bool
isPossibleFourthPower n =
biSqRes256 `unsafeAt` (fromIntegral n .&. 255)
&& biSqRes425 `unsafeAt` (fromIntegral (n `rem` 425))
&& biSqRes377 `unsafeAt` (fromIntegral (n `rem` 377))
newton4 :: Integral a => a -> a -> a
newton4 n a = go (step a)
where
step k = (3*k + n `quot` (k*k*k)) `quot` 4
go k
| m < k = go m
| otherwise = k
where
m = step k
approxBiSqrt :: Integral a => a -> a
approxBiSqrt = fromInteger . appBiSqrt . fromIntegral
#if WORD_SIZE_IN_BITS == 64
#define THRESH 7
#else
#define THRESH 13
#endif
appBiSqrt :: Integer -> Integer
appBiSqrt (S# i#) = S# (double2Int# (sqrtDouble# (sqrtDouble# (int2Double# i#))))
appBiSqrt n@(J# s# _)
| s# <# THRESH# = floor (sqrt . sqrt $ fromInteger n :: Double)
| otherwise = case integerLog2# n of
l# -> case uncheckedIShiftRA# l# 2# -# 47# of
h# -> case shiftRInteger n (4# *# h#) of
m -> case floor (sqrt $ sqrt $ fromInteger m :: Double) of
r -> shiftLInteger r h#
biSqRes256 :: UArray Int Bool
biSqRes256 = runSTUArray $ do
ar <- newArray (0,255) False
let note 257 = return ar
note i = unsafeWrite ar i True >> note (i+16)
unsafeWrite ar 0 True
unsafeWrite ar 16 True
note 1
biSqRes425 :: UArray Int Bool
biSqRes425 = runSTUArray $ do
ar <- newArray (0,424) False
let note 154 = return ar
note i = unsafeWrite ar ((i*i*i*i) `rem` 425) True >> note (i+1)
note 0
biSqRes377 :: UArray Int Bool
biSqRes377 = runSTUArray $ do
ar <- newArray (0,376) False
let note 144 = return ar
note i = unsafeWrite ar ((i*i*i*i) `rem` 377) True >> note (i+1)
note 0
biSqrtInt :: Int -> Int
biSqrtInt 0 = 0
biSqrtInt n
#if WORD_SIZE_IN_BITS == 64
| r > 55108 = 55108
#else
| r > 215 = 215
#endif
| n < r4 = r1
| otherwise = r
where
x :: Double
x = fromIntegral n
r = truncate (sqrt (sqrt x))
r2 = r*r
r4 = r2*r2
biSqrtWord :: Word -> Word
biSqrtWord 0 = 0
biSqrtWord n
#if WORD_SIZE_IN_BITS == 64
| r > 65535 = 65535
#else
| r > 255 = 255
#endif
| n < r4 = r1
| otherwise = r
where
x :: Double
x = fromIntegral n
r = truncate (sqrt (sqrt x))
r2 = r*r
r4 = r2*r2