-- |
-- Module:      Data.Chimera.WheelMapping
-- Copyright:   (c) 2017 Andrew Lelechenko
-- Licence:     MIT
-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
--
-- Helpers for mapping to <http://mathworld.wolfram.com/RoughNumber.html rough numbers>
-- and back. This has various applications in number theory.
--
-- __Example__
--
-- Let @isPrime@ be an expensive predicate,
-- which checks whether its argument is a prime number.
-- We can memoize it as usual:
--
-- > isPrimeCache1 :: UChimera Bool
-- > isPrimeCache1 = tabulate isPrime
-- >
-- > isPrime1 :: Word -> Bool
-- > isPrime1 = index isPrimeCache1
--
-- But one may argue that since the only even prime number is 2,
-- it is quite wasteful to cache @isPrime@ for even arguments.
-- So we can save half the space by memoizing it for odd
-- numbers only:
--
-- > isPrimeCache2 :: UChimera Bool
-- > isPrimeCache2 = tabulate (isPrime . (\n -> 2 * n + 1))
-- >
-- > isPrime2 :: Word -> Bool
-- > isPrime2 n
-- >   | n == 2    = True
-- >   | even n    = False
-- >   | otherwise = index isPrimeCache2 ((n - 1) `quot` 2)
--
-- Here @\\n -> 2 * n + 1@ maps n to the (n+1)-th odd number,
-- and @\\n -> (n - 1) \`quot\` 2@ takes it back. These functions
-- are available below as 'fromWheel2' and 'toWheel2'.
--
-- Odd numbers are the simplest example of numbers, lacking
-- small prime factors (so called
-- <http://mathworld.wolfram.com/RoughNumber.html rough numbers>).
-- Removing numbers, having small prime factors, is sometimes
-- called <https://en.wikipedia.org/wiki/Wheel_factorization wheel sieving>.
--
-- One can go further and exclude not only even numbers,
-- but also integers, divisible by 3.
-- To do this we need a function which maps n to the (n+1)-th number coprime with 2 and 3
-- (thus, with 6) and its inverse: namely, 'fromWheel6' and 'toWheel6'. Then write
--
-- > isPrimeCache6 :: UChimera Bool
-- > isPrimeCache6 = tabulate (isPrime . fromWheel6)
-- >
-- > isPrime6 :: Word -> Bool
-- > isPrime6 n
-- >   | n `elem` [2, 3] = True
-- >   | n `gcd` 6 /= 1  = False
-- >   | otherwise       = index isPrimeCache6 (toWheel6 n)
--
-- Thus, the wheel of 6 saves more space, improving memory locality.
--
-- (If you need to reduce memory consumption even further,
-- consider using 'Data.Bit.Bit' wrapper,
-- which provides an instance of unboxed vector,
-- packing one boolean per bit instead of one boolean per byte for 'Bool')
--

{-# LANGUAGE BangPatterns  #-}
{-# LANGUAGE CPP           #-}
{-# LANGUAGE MagicHash     #-}
{-# LANGUAGE UnboxedTuples #-}

module Data.Chimera.WheelMapping
  ( fromWheel2
  , toWheel2
  , fromWheel6
  , toWheel6
  , fromWheel30
  , toWheel30
  , fromWheel210
  , toWheel210
  ) where

import Data.Bits
import Data.Chimera.Compat
import GHC.Exts

bits :: Int
bits :: Int
bits = Word -> Int
fbs (Word
0 :: Word)

-- | Left inverse for 'fromWheel2'. Monotonically non-decreasing function.
--
-- prop> toWheel2 . fromWheel2 == id
toWheel2 :: Word -> Word
toWheel2 :: Word -> Word
toWheel2 Word
i = Word
i Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftR` Int
1
{-# INLINE toWheel2 #-}

-- | 'fromWheel2' n is the (n+1)-th positive odd number.
-- Sequence <https://oeis.org/A005408 A005408>.
--
-- prop> map fromWheel2 [0..] == [ n | n <- [0..], n `gcd` 2 == 1 ]
--
-- >>> map fromWheel2 [0..9]
-- [1,3,5,7,9,11,13,15,17,19]
fromWheel2 :: Word -> Word
fromWheel2 :: Word -> Word
fromWheel2 Word
i = Word
i Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftL` Int
1 Word -> Word -> Word
forall a. Num a => a -> a -> a
+ Word
1
{-# INLINE fromWheel2 #-}

-- | Left inverse for 'fromWheel6'. Monotonically non-decreasing function.
--
-- prop> toWheel6 . fromWheel6 == id
toWheel6 :: Word -> Word
toWheel6 :: Word -> Word
toWheel6 i :: Word
i@(W# Word#
i#) = case Int
bits of
  Int
64 -> Word# -> Word
W# Word#
z1# Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftR` Int
1
  Int
_  -> Word
i Word -> Word -> Word
forall a. Integral a => a -> a -> a
`quot` Word
3
  where
    m# :: Word#
m# = Word#
12297829382473034411## -- (2^65+1) / 3
    !(# Word#
z1#, Word#
_ #) = Word# -> Word# -> (# Word#, Word# #)
timesWord2# Word#
m# Word#
i#

{-# INLINE toWheel6 #-}

-- | 'fromWheel6' n is the (n+1)-th positive number, not divisible by 2 or 3.
-- Sequence <https://oeis.org/A007310 A007310>.
--
-- prop> map fromWheel6 [0..] == [ n | n <- [0..], n `gcd` 6 == 1 ]
--
-- >>> map fromWheel6 [0..9]
-- [1,5,7,11,13,17,19,23,25,29]
fromWheel6 :: Word -> Word
fromWheel6 :: Word -> Word
fromWheel6 Word
i = Word
i Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftL` Int
1 Word -> Word -> Word
forall a. Num a => a -> a -> a
+ Word
i Word -> Word -> Word
forall a. Num a => a -> a -> a
+ (Word
i Word -> Word -> Word
forall a. Bits a => a -> a -> a
.&. Word
1) Word -> Word -> Word
forall a. Num a => a -> a -> a
+ Word
1
{-# INLINE fromWheel6 #-}

-- | Left inverse for 'fromWheel30'. Monotonically non-decreasing function.
--
-- prop> toWheel30 . fromWheel30 == id
toWheel30 :: Word -> Word
toWheel30 :: Word -> Word
toWheel30 i :: Word
i@(W# Word#
i#) = Word
q Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftL` Int
3 Word -> Word -> Word
forall a. Num a => a -> a -> a
+ (Word
r Word -> Word -> Word
forall a. Num a => a -> a -> a
+ Word
r Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftR` Int
4) Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftR` Int
2
  where
    (Word
q, Word
r) = case Int
bits of
      Int
64 -> (Word
q64, Word
r64)
      Int
_  -> Word
i Word -> Word -> (Word, Word)
forall a. Integral a => a -> a -> (a, a)
`quotRem` Word
30

    m# :: Word#
m# = Word#
9838263505978427529## -- (2^67+7) / 15
    !(# Word#
z1#, Word#
_ #) = Word# -> Word# -> (# Word#, Word# #)
timesWord2# Word#
m# Word#
i#
    q64 :: Word
q64 = Word# -> Word
W# Word#
z1# Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftR` Int
4
    r64 :: Word
r64 = Word
i Word -> Word -> Word
forall a. Num a => a -> a -> a
- Word
q64 Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftL` Int
5 Word -> Word -> Word
forall a. Num a => a -> a -> a
+ Word
q64 Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftL` Int
1

{-# INLINE toWheel30 #-}

-- | 'fromWheel30' n is the (n+1)-th positive number, not divisible by 2, 3 or 5.
-- Sequence <https://oeis.org/A007775 A007775>.
--
-- prop> map fromWheel30 [0..] == [ n | n <- [0..], n `gcd` 30 == 1 ]
--
-- >>> map fromWheel30 [0..9]
-- [1,7,11,13,17,19,23,29,31,37]
fromWheel30 :: Word -> Word
fromWheel30 :: Word -> Word
fromWheel30 Word
i = ((Word
i Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftL` Int
2 Word -> Word -> Word
forall a. Num a => a -> a -> a
- Word
i Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftR` Int
2) Word -> Word -> Word
forall a. Bits a => a -> a -> a
.|. Word
1)
              Word -> Word -> Word
forall a. Num a => a -> a -> a
+ ((Word
i Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftL` Int
1 Word -> Word -> Word
forall a. Num a => a -> a -> a
- Word
i Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftR` Int
1) Word -> Word -> Word
forall a. Bits a => a -> a -> a
.&. Word
2)
{-# INLINE fromWheel30 #-}

-- | Left inverse for 'fromWheel210'. Monotonically non-decreasing function.
--
-- prop> toWheel210 . fromWheel210 == id
toWheel210 :: Word -> Word
toWheel210 :: Word -> Word
toWheel210 i :: Word
i@(W# Word#
i#) = Word
q Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftL` Int
5 Word -> Word -> Word
forall a. Num a => a -> a -> a
+ Word
q Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftL` Int
4 Word -> Word -> Word
forall a. Num a => a -> a -> a
+ Word# -> Word
W# Word#
tableEl#
  where
    !(Word
q, W# Word#
r#) = case Int
bits of
      Int
64 -> (Word
q64, Word
r64)
      Int
_  -> Word
i Word -> Word -> (Word, Word)
forall a. Integral a => a -> a -> (a, a)
`quotRem` Word
210

    m# :: Word#
m# = Word#
5621864860559101445## -- (2^69+13) / 105
    !(# Word#
z1#, Word#
_ #) = Word# -> Word# -> (# Word#, Word# #)
timesWord2# Word#
m# Word#
i#
    q64 :: Word
q64 = Word# -> Word
W# Word#
z1# Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftR` Int
6
    r64 :: Word
r64 = Word
i Word -> Word -> Word
forall a. Num a => a -> a -> a
- Word
q64 Word -> Word -> Word
forall a. Num a => a -> a -> a
* Word
210

    tableEl# :: Word#
tableEl# =
#if MIN_VERSION_base(4,16,0)
      word8ToWord#
#endif
      (Addr# -> Int# -> Word#
indexWord8OffAddr# Addr#
table# (Word# -> Int#
word2Int# Word#
r#))

    table# :: Addr#
    table# :: Addr#
table# = Addr#
"\NUL\NUL\NUL\NUL\NUL\NUL\NUL\NUL\NUL\NUL\NUL\SOH\SOH\STX\STX\STX\STX\ETX\ETX\EOT\EOT\EOT\EOT\ENQ\ENQ\ENQ\ENQ\ENQ\ENQ\ACK\ACK\a\a\a\a\a\a\b\b\b\b\t\t\n\n\n\n\v\v\v\v\v\v\f\f\f\f\f\f\r\r\SO\SO\SO\SO\SO\SO\SI\SI\SI\SI\DLE\DLE\DC1\DC1\DC1\DC1\DC1\DC1\DC2\DC2\DC2\DC2\DC3\DC3\DC3\DC3\DC3\DC3\DC4\DC4\DC4\DC4\DC4\DC4\DC4\DC4\NAK\NAK\NAK\NAK\SYN\SYN\ETB\ETB\ETB\ETB\CAN\CAN\EM\EM\EM\EM\SUB\SUB\SUB\SUB\SUB\SUB\SUB\SUB\ESC\ESC\ESC\ESC\ESC\ESC\FS\FS\FS\FS\GS\GS\GS\GS\GS\GS\RS\RS\US\US\US\US      !!\"\"\"\"\"\"######$$$$%%&&&&''''''(())))))****++,,,,--........../"#
    -- map Data.Char.chr [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 13, 13, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 16, 16, 17, 17, 17, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 20, 20, 20, 21, 21, 21, 21, 22, 22, 23, 23, 23, 23, 24, 24, 25, 25, 25, 25, 26, 26, 26, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 27, 28, 28, 28, 28, 29, 29, 29, 29, 29, 29, 30, 30, 31, 31, 31, 31, 32, 32, 32, 32, 32, 32, 33, 33, 34, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 35, 36, 36, 36, 36, 37, 37, 38, 38, 38, 38, 39, 39, 39, 39, 39, 39, 40, 40, 41, 41, 41, 41, 41, 41, 42, 42, 42, 42, 43, 43, 44, 44, 44, 44, 45, 45, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 47]

{-# INLINE toWheel210 #-}

-- | 'fromWheel210' n is the (n+1)-th positive number, not divisible by 2, 3, 5 or 7.
-- Sequence <https://oeis.org/A008364 A008364>.
--
-- prop> map fromWheel210 [0..] == [ n | n <- [0..], n `gcd` 210 == 1 ]
--
-- >>> map fromWheel210 [0..9]
-- [1,11,13,17,19,23,29,31,37,41]
fromWheel210 :: Word -> Word
fromWheel210 :: Word -> Word
fromWheel210 i :: Word
i@(W# Word#
i#) = Word
q Word -> Word -> Word
forall a. Num a => a -> a -> a
* Word
210 Word -> Word -> Word
forall a. Num a => a -> a -> a
+ Word# -> Word
W# Word#
tableEl#
  where
    !(Word
q, W# Word#
r#) = case Int
bits of
      Int
64 -> (Word
q64, Word
r64)
      Int
_  -> Word
i Word -> Word -> (Word, Word)
forall a. Integral a => a -> a -> (a, a)
`quotRem` Word
48

    m# :: Word#
m# = Word#
12297829382473034411## -- (2^65+1) / 3
    !(# Word#
z1#, Word#
_ #) = Word# -> Word# -> (# Word#, Word# #)
timesWord2# Word#
m# Word#
i#
    q64 :: Word
q64 = Word# -> Word
W# Word#
z1# Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftR` Int
5
    r64 :: Word
r64 = Word
i Word -> Word -> Word
forall a. Num a => a -> a -> a
- Word
q64 Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftL` Int
5 Word -> Word -> Word
forall a. Num a => a -> a -> a
- Word
q64 Word -> Int -> Word
forall a. Bits a => a -> Int -> a
`shiftL` Int
4

    tableEl# :: Word#
tableEl# =
#if MIN_VERSION_base(4,16,0)
      word8ToWord#
#endif
      (Addr# -> Int# -> Word#
indexWord8OffAddr# Addr#
table# (Word# -> Int#
word2Int# Word#
r#))

    table# :: Addr#
    table# :: Addr#
table# = Addr#
"\SOH\v\r\DC1\DC3\ETB\GS\US%)+/5;=CGIOSYaegkmqy\DEL\131\137\139\143\149\151\157\163\167\169\173\179\181\187\191\193\197\199\209"#
    -- map Data.Char.chr [1, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 121, 127, 131, 137, 139, 143, 149, 151, 157, 163, 167, 169, 173, 179, 181, 187, 191, 193, 197, 199, 209]

{-# INLINE fromWheel210 #-}