{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE MagicHash #-}
module Math.NumberTheory.Roots.Cubes
( integerCubeRoot
, integerCubeRoot'
, exactCubeRoot
, isCube
, isCube'
, isPossibleCube
) where
import Data.Bits (finiteBitSize, (.&.))
import GHC.Exts (Int#, Ptr(..), int2Double#, double2Int#, isTrue#, (/##), (**##), (<#))
import Numeric.Natural (Natural)
#ifdef MIN_VERSION_integer_gmp
import GHC.Exts (quotInt#, (*#), (-#))
import GHC.Integer.GMP.Internals (Integer(..), shiftLInteger, shiftRInteger, sizeofBigNat#)
import GHC.Integer.Logarithms (integerLog2#)
#define IS S#
#define IP Jp#
#define bigNatSize sizeofBigNat
#else
import GHC.Exts (minusWord#, timesWord#, quotWord#)
import GHC.Num.BigNat (bigNatSize#)
import GHC.Num.Integer (Integer(..), integerLog2#, integerShiftR#, integerShiftL#)
#endif
import Math.NumberTheory.Utils.BitMask (indexBitSet)
{-# SPECIALISE integerCubeRoot :: Int -> Int,
Word -> Word,
Integer -> Integer,
Natural -> Natural
#-}
integerCubeRoot :: Integral a => a -> a
integerCubeRoot :: a -> a
integerCubeRoot a
0 = a
0
integerCubeRoot a
n
| a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0 = a -> a
forall a. Integral a => a -> a
integerCubeRoot' a
n
| Bool
otherwise =
let m :: a
m = a -> a
forall a. Num a => a -> a
negate a
n
r :: a
r = if a
m a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0
then a -> a
forall a. Num a => a -> a
negate (a -> a) -> (Integer -> a) -> Integer -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> a
forall a. Num a => Integer -> a
fromInteger (Integer -> a) -> Integer -> a
forall a b. (a -> b) -> a -> b
$ Integer -> Integer
forall a. Integral a => a -> a
integerCubeRoot' (Integer -> Integer
forall a. Num a => a -> a
negate (Integer -> Integer) -> Integer -> Integer
forall a b. (a -> b) -> a -> b
$ a -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n)
else a -> a
forall a. Num a => a -> a
negate (a -> a
forall a. Integral a => a -> a
integerCubeRoot' a
m)
in if a
ra -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
*a
r a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
n then a
r else (a
ra -> a -> a
forall a. Num a => a -> a -> a
-a
1)
{-# RULES
"integerCubeRoot'/Int" integerCubeRoot' = cubeRootInt'
"integerCubeRoot'/Word" integerCubeRoot' = cubeRootWord
"integerCubeRoot'/Integer" integerCubeRoot' = cubeRootIgr
#-}
{-# INLINE [1] integerCubeRoot' #-}
integerCubeRoot' :: Integral a => a -> a
integerCubeRoot' :: a -> a
integerCubeRoot' a
0 = a
0
integerCubeRoot' a
n = a -> a -> a
forall a. Integral a => a -> a -> a
newton3 a
n (a -> a
forall a. Integral a => a -> a
approxCuRt a
n)
{-# SPECIALISE exactCubeRoot :: Int -> Maybe Int,
Word -> Maybe Word,
Integer -> Maybe Integer,
Natural -> Maybe Natural
#-}
exactCubeRoot :: Integral a => a -> Maybe a
exactCubeRoot :: a -> Maybe a
exactCubeRoot a
0 = a -> Maybe a
forall a. a -> Maybe a
Just a
0
exactCubeRoot a
n
| a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 =
if a
m a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0
then (Integer -> a) -> Maybe Integer -> Maybe a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (a -> a
forall a. Num a => a -> a
negate (a -> a) -> (Integer -> a) -> Integer -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> a
forall a. Num a => Integer -> a
fromInteger) (Maybe Integer -> Maybe a) -> Maybe Integer -> Maybe a
forall a b. (a -> b) -> a -> b
$ Integer -> Maybe Integer
forall a. Integral a => a -> Maybe a
exactCubeRoot (Integer -> Integer
forall a. Num a => a -> a
negate (Integer -> Integer) -> Integer -> Integer
forall a b. (a -> b) -> a -> b
$ a -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n)
else (a -> a) -> Maybe a -> Maybe a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> a
forall a. Num a => a -> a
negate (a -> Maybe a
forall a. Integral a => a -> Maybe a
exactCubeRoot a
m)
| a -> Bool
forall a. Integral a => a -> Bool
isPossibleCube a
n Bool -> Bool -> Bool
&& a
ra -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
*a
r a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
n = a -> Maybe a
forall a. a -> Maybe a
Just a
r
| Bool
otherwise = Maybe a
forall a. Maybe a
Nothing
where
m :: a
m = a -> a
forall a. Num a => a -> a
negate a
n
r :: a
r = a -> a
forall a. Integral a => a -> a
integerCubeRoot' a
n
{-# SPECIALISE isCube :: Int -> Bool,
Word -> Bool,
Integer -> Bool,
Natural -> Bool
#-}
isCube :: Integral a => a -> Bool
isCube :: a -> Bool
isCube a
0 = Bool
True
isCube a
n
| a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0 = a -> Bool
forall a. Integral a => a -> Bool
isCube' a
n
| a
m a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0 = a -> Bool
forall a. Integral a => a -> Bool
isCube' a
m
| Bool
otherwise = Integer -> Bool
forall a. Integral a => a -> Bool
isCube' (Integer -> Integer
forall a. Num a => a -> a
negate (a -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n) :: Integer)
where
m :: a
m = a -> a
forall a. Num a => a -> a
negate a
n
{-# SPECIALISE isCube' :: Int -> Bool,
Word -> Bool,
Integer -> Bool,
Natural -> Bool
#-}
isCube' :: Integral a => a -> Bool
isCube' :: a -> Bool
isCube' !a
n = a -> Bool
forall a. Integral a => a -> Bool
isPossibleCube a
n
Bool -> Bool -> Bool
&& (a
ra -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
*a
r a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
n)
where
r :: a
r = a -> a
forall a. Integral a => a -> a
integerCubeRoot' a
n
{-# SPECIALISE isPossibleCube :: Int -> Bool,
Word -> Bool,
Integer -> Bool,
Natural -> Bool
#-}
isPossibleCube :: Integral a => a -> Bool
isPossibleCube :: a -> Bool
isPossibleCube a
n'
= Ptr Word -> Int -> Bool
indexBitSet Ptr Word
mask512 (Integer -> Int
forall a. Num a => Integer -> a
fromInteger (Integer
n Integer -> Integer -> Integer
forall a. Bits a => a -> a -> a
.&. Integer
511))
Bool -> Bool -> Bool
&& Ptr Word -> Int -> Bool
indexBitSet Ptr Word
mask837 (Integer -> Int
forall a. Num a => Integer -> a
fromInteger (Integer
n Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
837))
Bool -> Bool -> Bool
&& Ptr Word -> Int -> Bool
indexBitSet Ptr Word
mask637 (Integer -> Int
forall a. Num a => Integer -> a
fromInteger (Integer
n Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
637))
Bool -> Bool -> Bool
&& Ptr Word -> Int -> Bool
indexBitSet Ptr Word
mask703 (Integer -> Int
forall a. Num a => Integer -> a
fromInteger (Integer
n Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
703))
where
n :: Integer
n = a -> Integer
forall a. Integral a => a -> Integer
toInteger a
n'
cubeRootInt' :: Int -> Int
cubeRootInt' :: Int -> Int
cubeRootInt' Int
0 = Int
0
cubeRootInt' Int
n
| Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
c Bool -> Bool -> Bool
|| Int
c Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1
| Int
0 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
d Bool -> Bool -> Bool
&& Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n = Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
| Bool
otherwise = Int
r
where
x :: Double
x = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n :: Double
r :: Int
r = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
truncate (Double
x Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
3))
c :: Int
c = Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
r
d :: Int
d = Int
cInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
cubeRootWordLimit :: Word
cubeRootWordLimit :: Word
cubeRootWordLimit = if Word -> Int
forall b. FiniteBits b => b -> Int
finiteBitSize (Word
0 :: Word) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
64 then Word
2642245 else Word
1625
cubeRootWord :: Word -> Word
cubeRootWord :: Word -> Word
cubeRootWord Word
0 = Word
0
cubeRootWord Word
w
| Word
r Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
> Word
cubeRootWordLimit = Word
cubeRootWordLimit
| Word
w Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
< Word
c = Word
rWord -> Word -> Word
forall a. Num a => a -> a -> a
-Word
1
| Word
c Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
< Word
w Bool -> Bool -> Bool
&& Word
e Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
< Word
w Bool -> Bool -> Bool
&& Word
c Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
< Word
e = Word
rWord -> Word -> Word
forall a. Num a => a -> a -> a
+Word
1
| Bool
otherwise = Word
r
where
r :: Word
r = Double -> Word
forall a b. (RealFrac a, Integral b) => a -> b
truncate ((Word -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Word
w) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
3) :: Double)
c :: Word
c = Word
rWord -> Word -> Word
forall a. Num a => a -> a -> a
*Word
rWord -> Word -> Word
forall a. Num a => a -> a -> a
*Word
r
d :: Word
d = Word
3Word -> Word -> Word
forall a. Num a => a -> a -> a
*Word
rWord -> Word -> Word
forall a. Num a => a -> a -> a
*(Word
rWord -> Word -> Word
forall a. Num a => a -> a -> a
+Word
1)
e :: Word
e = Word
cWord -> Word -> Word
forall a. Num a => a -> a -> a
+Word
d
cubeRootIgr :: Integer -> Integer
cubeRootIgr :: Integer -> Integer
cubeRootIgr Integer
0 = Integer
0
cubeRootIgr Integer
n = Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
newton3 Integer
n (Integer -> Integer
forall a. Integral a => a -> a
approxCuRt Integer
n)
{-# SPECIALISE newton3 :: Integer -> Integer -> Integer #-}
newton3 :: Integral a => a -> a -> a
newton3 :: a -> a -> a
newton3 a
n a
a = a -> a
go (a -> a
step a
a)
where
step :: a -> a
step a
k = (a
2a -> a -> a
forall a. Num a => a -> a -> a
*a
k a -> a -> a
forall a. Num a => a -> a -> a
+ a
n a -> a -> a
forall a. Integral a => a -> a -> a
`quot` (a
ka -> a -> a
forall a. Num a => a -> a -> a
*a
k)) a -> a -> a
forall a. Integral a => a -> a -> a
`quot` a
3
go :: a -> a
go a
k
| a
m a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
k = a -> a
go a
m
| Bool
otherwise = a
k
where
m :: a
m = a -> a
step a
k
{-# SPECIALISE approxCuRt :: Integer -> Integer #-}
approxCuRt :: Integral a => a -> a
approxCuRt :: a -> a
approxCuRt a
0 = a
0
approxCuRt a
n = Integer -> a
forall a. Num a => Integer -> a
fromInteger (Integer -> a) -> Integer -> a
forall a b. (a -> b) -> a -> b
$ Integer -> Integer
appCuRt (a -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n)
appCuRt :: Integer -> Integer
appCuRt :: Integer -> Integer
appCuRt (IS Int#
i#) = case Double# -> Int#
double2Int# (Int# -> Double#
int2Double# Int#
i# Double# -> Double# -> Double#
**## (Double#
1.0## Double# -> Double# -> Double#
/## Double#
3.0##)) of
Int#
r# -> Int# -> Integer
IS Int#
r#
appCuRt n :: Integer
n@(IP bn#)
| Int# -> Bool
isTrue# ((bigNatSize# bn#) <# thresh#) =
Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
floor (Integer -> Double
forall a. Num a => Integer -> a
fromInteger Integer
n Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1.0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
3.0) :: Double)
| Bool
otherwise = case Integer -> Int#
integerLog2# Integer
n of
#ifdef MIN_VERSION_integer_gmp
Int#
l# -> case (Int#
l# Int# -> Int# -> Int#
`quotInt#` Int#
3#) Int# -> Int# -> Int#
-# Int#
51# of
Int#
h# -> case Integer -> Int# -> Integer
shiftRInteger Integer
n (Int#
3# Int# -> Int# -> Int#
*# Int#
h#) of
Integer
m -> case Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
floor (Integer -> Double
forall a. Num a => Integer -> a
fromInteger Integer
m Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1.0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
3.0) :: Double) of
Integer
r -> Integer -> Int# -> Integer
shiftLInteger Integer
r Int#
h#
#else
l# -> case (l# `quotWord#` 3##) `minusWord#` 51## of
h# -> case integerShiftR# n (3## `timesWord#` h#) of
m -> case floor (fromInteger m ** (1.0/3.0) :: Double) of
r -> integerShiftL# r h#
#endif
where
thresh# :: Int#
thresh# :: Int#
thresh# = if Word -> Int
forall b. FiniteBits b => b -> Int
finiteBitSize (Word
0 :: Word) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
64 then Int#
5# else Int#
9#
appCuRt Integer
_ = [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error [Char]
"integerCubeRoot': negative argument"
mask512 :: Ptr Word
mask512 :: Ptr Word
mask512 = Addr# -> Ptr Word
forall a. Addr# -> Ptr a
Ptr Addr#
"\171\171\170\171\170\171\170\171\171\171\170\171\170\171\170\171\170\171\170\171\170\171\170\171\171\171\170\171\170\171\170\171\170\171\170\171\170\171\170\171\171\171\170\171\170\171\170\171\170\171\170\171\170\171\170\171\171\171\170\171\170\171\170\171"#
mask837 :: Ptr Word
mask837 :: Ptr Word
mask837 = Addr# -> Ptr Word
forall a. Addr# -> Ptr a
Ptr Addr#
"\ETX\SOH\NUL\b\b@@@\SOH\NUL\NUL\n\NUL0\DLE \NUL\NUL\NUL\EOT\b\EOT\NULP\NUL\NUL\128\ETX\NUL\STX\DLE\NUL\NUL\128@\129\NUL\NUL\NUL\EOT \NUL\160\NUL\NUL\NUL\ENQ\NUL\DLE\b0\NUL\NUL\NUL\ETX\EOT\STX\NUL(\NUL\NUL@\SOH\NUL\SOH\b\NUL\NUL@\160@\NUL\NUL\NUL\STX\DLE\NULp\NUL\NUL\128\STX\NUL\b\EOT\b\NUL\NUL\NUL\SOH\STX\ETX\NUL\DC4\NUL\NUL\160\128\128\NUL\EOT\EOT\NUL \DLE"#
mask637 :: Ptr Word
mask637 :: Ptr Word
mask637 = Addr# -> Ptr Word
forall a. Addr# -> Ptr a
Ptr Addr#
"\ETX!\NUL\b\EOT\NUL\NUL\STX\SOH@\b\DC4\b\SOH@ \NUL\NUL\DLE\b\NULB\160@\CAN\NUL\STX\SOH\NUL\128@\NUL\DLE\STX\ENQB@\DLE\b\NUL\NUL\EOT\130\128\DLE(\DLE\STX\128@\NUL\NUL \DLE\NUL\134@\129\DLE\NUL\EOT\STX\NUL\NUL\129\NUL \EOT\n\132\NUL \DLE\NUL\NUL\b\EOT\NUL!\DLE"#
mask703 :: Ptr Word
mask703 :: Ptr Word
mask703 = Addr# -> Ptr Word
forall a. Addr# -> Ptr a
Ptr Addr#
"\ETX\t\NUL\140` \NUL\NUL\DC1\b\DLE\SOH\128\NUL\NUL&@\DLE\NUL\128\NUL\b\b\128\NUL\NUL\SOH!\DLE\DLE\NUL\SOH\f\n\STX`\NUL\ETX\SOH\NUL\f@\160\NUL\NUL\ENQ\STX0\NUL\128\192\NUL\ACK@P0\128\NUL\b\b\132\128\NUL\NUL\SOH\DLE\DLE\NUL\SOH\NUL\b\STXd\NUL\NUL\SOH\128\b\DLE\136\NUL\NUL\EOT\ACK1\NUL\144@"#