{-# LANGUAGE NoImplicitPrelude #-}
{-
Some of these functions might be moved to NumericPrelude.

Wikipedia: (primitive) roots of unity modulo n
   (primitive) roots must be units and all units are (primitive) roots
   maximum possible order for primitive roots - Carmichael
   all possible orders: divisor of Carmichael (proof? statement already in Carmichael-function-article)
   sum of primitive roots that vanishes
   order of primitive root is a divisor of each possible exponent
      proof with GCD and diophantine in exponent
   check for primitive root: fast exponentiation,
      primitivity: check exponents that are prime divisors
   how to find a primitive root: just try
   sum of powers of a primitive root is zero
   number of primitive roots of given order
      g(n,k) > 0 if k|lambda(n)
      g(n,k) = 0 else
      g(n,1) = 1
      g(4,2) = 1
      g(2^n,2) = 3 for n>=3  ((-1) is always a square root of 1)
      g(2^n,2^k) = 2^k for k>=2 && k<n-1
      g(n,2) = 1 for n>=3 and n in http://oeis.org/A033948
      sum(g(n,k), k\in\N) = phi(n)
      There are only a few patterns that occur as rows of g,
      but a row of g (i.e. g(n)) does functionally depend on
      either lambda(n) nor phi(n)
      lambda(14) = 6   nozeros(g(14)) = [1,1,2,2]   (f ~ [1,2,3,6])
      lambda(21) = 6   nozeros(g(21)) = [1,3,2,6]   (f ~ [1,4,3,12])
      phi(13) = 12   nozeros(g(13)) = [1,1,2,2,2,4]   (f ~ [1,2,3,4,6,12])
      phi(21) = 12   nozeros(g(21)) = [1,3,2,6]       (f ~ [1,4,3,12])
      However length(nozeros(f(n))) = numberofdivisors(lambda(n))
      numberofdivisors=A000005
   number of roots of given order
      easier to compute
      k|m => f(n,k) | f(n,m)
      g(n,k) = f(n,k) - sum(f(n,d), d|k and k/d prime) + ...
         inclusion-exclusion-principle
      better to write the other round:
      f(n,k) = sum(g(n,d), d|k) - this is Dirichlet convolution
      RUNM says f(n,k) is multiplicative
         list it in multiplicative function
      f(n,1) = 1 for n>=2
      f(n,lambda(n)) = phi(n)
      f(n,a·b) = f(n,a)·f(n,b) if a and b are coprime (conjecture)
      f(n,lcm(a,b)) = lcm(f(n,a),f(n,b)) (conjecture)
      If this conjecture is true, then we only need to know f(n,p^i).
      The following conjecture is wrong:
         for prime p it is   f(n,p^i) = gcd(lambda(n),p^i)
      counterexamples
         f(8,2) = 4, lambda(8)=2
         f(63,3) = 9, lambda(63)=6
         f(275,5) = 25, lambda(275)=20
         f(1247,7) = 49, lambda(1247)=84
      It seems to be:
         for prime p it is   f(n,p^i) = p^j for some j
   How to find a modulus where there is a primitive root of order o?
      just try numbers from the sequence o+1, 2*o+1, 3*o+1
      Because of [[Dirichlet's theorem on arithmetic progressions]]
      you will somewhen find a prime p,
      and its Carmichael value is p-1, which is a multiple of o.
      In this ring even 'o' is a unit.
   How to find a modulus where there are primitive roots of orders o1,..,ok?
      Just search for a modulus with roots of order lcm(o1,...,ok).
      The smallest such modulus should also be the smallest one
      that has primitive roots of orders o1,..,ok?
      Proof: If a ring has primitive roots of orders o1,..,ok
      then all orders divide the carmichael value of that ring,
      thus lcm(o1,...,ok) divides the carmichael value of that ring,
      thus there is a primitive root of order lcm(o1,...,ok).
-}
module Synthesizer.Basic.NumberTheory (
   fermatFactors,
   uniquePrimeFactors,
   primeFactors,
   multiplicativeGenerator,
   Order (Order, getOrder),
   PrimitiveRoot(primitiveRootCandidates, maximumOrderOfPrimitiveRootsOfUnity),
   primitiveRootsOfUnity,
   lcmMulti,
   primitiveRootsOfUnityFullOrbit,
   primitiveRootsOfOrbit,
   hasPrimitiveRootOfUnityNaive,
   ordersOfPrimitiveRootsOfUnityTest,
   orderOfOrbit,
   hasPrimitiveRootOfUnityInteger,
   ordersOfPrimitiveRootsOfUnityInteger,
   ordersOfRootsOfUnityInteger,
   ordersOfRootsOfUnityIntegerCondensed,
   rootsOfUnityPower,
   ringsWithPrimitiveRootOfUnityAndUnit,
   ringsWithPrimitiveRootsOfUnityAndUnitsNaive,
   ringWithPrimitiveRootsOfUnityAndUnits,
   ringWithPrimitiveRootsOfUnity,
   is3Smooth,
   is5Smooth,
   numbers3Smooth,
   numbers5Smooth,
   ceilingPowerOfTwo,
   ceilingPower,
   ceilingLog,
   powerOfTwoFactors,
   divideByMaximumPower,
   ceiling3Smooth,
   ceiling5Smooth,
   isPrime,
   raderWorstCases,
   fastFourierRing,

   -- for testing
   multiplicativeGeneratorSet,
   multiplicativeGeneratorDivisors,
   primitiveRootsOfUnityPower,
   primitiveRootsOfUnityNaive,
   primitiveRootsOfUnityFullOrbitTest,
   maximumOrderOfPrimitiveRootsOfUnityNaive,
   maximumOrderOfPrimitiveRootsOfUnityInteger,
   divideByMaximumPowerRecursive,
   numbers3SmoothCorec,
   numbers3SmoothFoldr,
   numbers3SmoothSet,
   numbers5SmoothCorec,
   numbers5SmoothFoldr,
   numbers5SmoothSet,
   ceiling3SmoothScan,
   ceiling5SmoothScan,
   ceiling3SmoothNaive,
   ceiling5SmoothNaive,
   ceiling3SmoothTrace,
   ceiling5SmoothTrace,
   ) where

import qualified Synthesizer.State.Signal as SigS

import qualified Data.Set as Set
import qualified Data.Map as Map

import qualified Algebra.Ring as Ring
import qualified Algebra.Units as Units
import qualified Algebra.PrincipalIdealDomain as PID
import qualified Algebra.IntegralDomain as Integral
import qualified Algebra.ZeroTestable as ZeroTestable

import qualified Number.ResidueClass.Check as RC
import Number.ResidueClass.Check ((/:), )

import qualified Number.FixedPoint as FP
import Data.Bits (Bits, (.&.), (.|.), shiftR, )

import qualified Data.List.HT as ListHT
import Data.List (unfoldr, mapAccumL, genericDrop, genericSplitAt, )
import Data.Tuple.HT (mapFst, mapSnd, mapPair, swap, )
import Data.Maybe.HT (toMaybe, )

import Test.QuickCheck (Arbitrary(arbitrary), )

import NumericPrelude.Numeric
import NumericPrelude.Base


{- |
The first pair member in @powerOfTwoFactors n@
is the maximum factor of @n@, that is a power of two.
-}
powerOfTwoFactors ::
   (Bits a, Integral.C a) => a -> (a, a)
powerOfTwoFactors :: forall a. (Bits a, C a) => a -> (a, a)
powerOfTwoFactors a
n =
   let powerOfTwo :: a
powerOfTwo = a
n a -> a -> a
forall a. Bits a => a -> a -> a
.&. (-a
n)
   in  (a
powerOfTwo, a -> a -> a
forall a. C a => a -> a -> a
div a
n a
powerOfTwo)


{- |
List all factorizations of an odd number
where the first factor is at most the second factor
and the first factors are in descending order.
-}
fermatFactors :: Integer -> [(Integer,Integer)]
fermatFactors :: Integer -> [(Integer, Integer)]
fermatFactors Integer
n =
   let root :: Integer
root = Integer -> Integer -> Integer
FP.sqrt Integer
1 Integer
n
   in  ((Integer, Integer) -> (Integer, Integer))
-> [(Integer, Integer)] -> [(Integer, Integer)]
forall a b. (a -> b) -> [a] -> [b]
map (\(Integer
a,Integer
b) -> (Integer
bInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
a,Integer
bInteger -> Integer -> Integer
forall a. C a => a -> a -> a
+Integer
a)) ([(Integer, Integer)] -> [(Integer, Integer)])
-> [(Integer, Integer)] -> [(Integer, Integer)]
forall a b. (a -> b) -> a -> b
$
       [(Integer, Integer)]
-> [(Integer, Integer)] -> [(Integer, Integer)]
forall a b c. Ord a => [(a, b)] -> [(a, c)] -> [(b, c)]
mergeAndFilter
          ([Integer] -> [Integer] -> [(Integer, Integer)]
forall a b. [a] -> [b] -> [(a, b)]
zip ((Integer -> Integer -> Integer)
-> Integer -> [Integer] -> [Integer]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Integer -> Integer -> Integer
forall a. C a => a -> a -> a
(+) Integer
n [Integer
1,Integer
3..]) [Integer
0 .. Integer -> Integer -> Integer
forall a. C a => a -> a -> a
div (Integer
nInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
1) Integer
2])
          ([Integer] -> [Integer] -> [(Integer, Integer)]
forall a b. [a] -> [b] -> [(a, b)]
zip ((Integer -> Integer -> Integer)
-> Integer -> [Integer] -> [Integer]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Integer -> Integer -> Integer
forall a. C a => a -> a -> a
(+) (Integer
rootInteger -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
root) ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer
2Integer -> Integer -> Integer
forall a. C a => a -> a -> a
+) (Integer
2Integer -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
rootInteger -> Integer -> Integer
forall a. C a => a -> a -> a
+Integer
1)) [Integer
root..])

mergeAndFilter :: (Ord a) => [(a,b)] -> [(a,c)] -> [(b,c)]
mergeAndFilter :: forall a b c. Ord a => [(a, b)] -> [(a, c)] -> [(b, c)]
mergeAndFilter ((a
a0,b
b):[(a, b)]
a0s) ((a
a1,c
c):[(a, c)]
a1s) =
   case a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare a
a0 a
a1 of
      Ordering
LT -> [(a, b)] -> [(a, c)] -> [(b, c)]
forall a b c. Ord a => [(a, b)] -> [(a, c)] -> [(b, c)]
mergeAndFilter [(a, b)]
a0s ((a
a1,c
c)(a, c) -> [(a, c)] -> [(a, c)]
forall a. a -> [a] -> [a]
:[(a, c)]
a1s)
      Ordering
GT -> [(a, b)] -> [(a, c)] -> [(b, c)]
forall a b c. Ord a => [(a, b)] -> [(a, c)] -> [(b, c)]
mergeAndFilter ((a
a0,b
b)(a, b) -> [(a, b)] -> [(a, b)]
forall a. a -> [a] -> [a]
:[(a, b)]
a0s) [(a, c)]
a1s
      Ordering
EQ -> (b
b,c
c) (b, c) -> [(b, c)] -> [(b, c)]
forall a. a -> [a] -> [a]
: [(a, b)] -> [(a, c)] -> [(b, c)]
forall a b c. Ord a => [(a, b)] -> [(a, c)] -> [(b, c)]
mergeAndFilter [(a, b)]
a0s [(a, c)]
a1s
mergeAndFilter [(a, b)]
_ [(a, c)]
_ = []



multiplicativeGenerator :: Integer -> Integer
multiplicativeGenerator :: Integer -> Integer
multiplicativeGenerator = Integer -> Integer
multiplicativeGeneratorDivisors

{- |
Argument must be a prime.
Usage of Set for efficient filtering of candidates seems to be overkill,
since the multiplicative generator seems to be small in most cases,
i.e. 2 or 3.

Smallest multiplicative generators for primes:
<http://oeis.org/A001918>

Especially large generators:
$ filter ((>31) . snd) $ map (\n -> (n, multiplicativeGenerator n)) $ tail NumberTheory.primes
[(36721,37),(48889,34),(51361,37),(55441,38),(63361,37),(64609,35),(71761,44),(88321,34),(92401,34),(93481,35),(95471,43),(97441,37),(104711,43),(110881,69)

$ filter ((>63) . snd) $ map (\n -> (n, multiplicativeGenerator n)) $ tail NumberTheory.primes
[(110881,69),(760321,73)

A solution with medium complexity
could at least observe the least 64 numbers using a Word64.
-}
multiplicativeGeneratorSet :: Integer -> Integer
multiplicativeGeneratorSet :: Integer -> Integer
multiplicativeGeneratorSet Integer
p =
   let search :: Set Integer -> Integer
search Set Integer
candidates =
          case Set Integer -> Maybe (Integer, Set Integer)
forall a. Set a -> Maybe (a, Set a)
Set.minView Set Integer
candidates of
             Maybe (Integer, Set Integer)
Nothing -> [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error ([Char] -> Integer) -> [Char] -> Integer
forall a b. (a -> b) -> a -> b
$ Integer -> [Char]
forall a. Show a => a -> [Char]
show Integer
p [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
" is not a prime"
             Just (Integer
x,Set Integer
rest) ->
                case T Integer -> Set Integer
forall a. Ord a => T a -> Set a
orbitSet (T Integer -> Set Integer) -> T Integer -> Set Integer
forall a b. (a -> b) -> a -> b
$ Integer -> Integer -> T Integer
forall a. C a => a -> a -> T a
orbit Integer
p Integer
x of
                   Set Integer
new ->
                      -- fromIntegral (Set.size new) == p-2
                      if Set Integer
new Set Integer -> Set Integer -> Bool
forall a. Eq a => a -> a -> Bool
== [Integer] -> Set Integer
forall a. Ord a => [a] -> Set a
Set.fromList [Integer
1..Integer
pInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
1]
                        then Integer
x
                        else Set Integer -> Integer
search (Set Integer -> Set Integer -> Set Integer
forall a. Ord a => Set a -> Set a -> Set a
Set.difference Set Integer
rest Set Integer
new)
   in  Set Integer -> Integer
search (Set Integer -> Integer) -> Set Integer -> Integer
forall a b. (a -> b) -> a -> b
$ [Integer] -> Set Integer
forall a. Ord a => [a] -> Set a
Set.fromList [Integer
1..Integer
pInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
1]

multiplicativeGeneratorDivisors :: Integer -> Integer
multiplicativeGeneratorDivisors :: Integer -> Integer
multiplicativeGeneratorDivisors Integer
p =
   [Integer] -> Integer
forall a. HasCallStack => [a] -> a
head ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$ Integer -> Order -> [Integer]
forall a. (PrimitiveRoot a, Eq a) => a -> Order -> [a]
primitiveRootsOfUnity Integer
p (Integer -> Order
Order (Integer -> Order) -> Integer -> Order
forall a b. (a -> b) -> a -> b
$ Integer
pInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
1)


newtype Order = Order {Order -> Integer
getOrder :: Integer}
   deriving (Int -> Order -> [Char] -> [Char]
[Order] -> [Char] -> [Char]
Order -> [Char]
(Int -> Order -> [Char] -> [Char])
-> (Order -> [Char]) -> ([Order] -> [Char] -> [Char]) -> Show Order
forall a.
(Int -> a -> [Char] -> [Char])
-> (a -> [Char]) -> ([a] -> [Char] -> [Char]) -> Show a
$cshowsPrec :: Int -> Order -> [Char] -> [Char]
showsPrec :: Int -> Order -> [Char] -> [Char]
$cshow :: Order -> [Char]
show :: Order -> [Char]
$cshowList :: [Order] -> [Char] -> [Char]
showList :: [Order] -> [Char] -> [Char]
Show, Order -> Order -> Bool
(Order -> Order -> Bool) -> (Order -> Order -> Bool) -> Eq Order
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Order -> Order -> Bool
== :: Order -> Order -> Bool
$c/= :: Order -> Order -> Bool
/= :: Order -> Order -> Bool
Eq, Eq Order
Eq Order =>
(Order -> Order -> Ordering)
-> (Order -> Order -> Bool)
-> (Order -> Order -> Bool)
-> (Order -> Order -> Bool)
-> (Order -> Order -> Bool)
-> (Order -> Order -> Order)
-> (Order -> Order -> Order)
-> Ord Order
Order -> Order -> Bool
Order -> Order -> Ordering
Order -> Order -> Order
forall a.
Eq a =>
(a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
$ccompare :: Order -> Order -> Ordering
compare :: Order -> Order -> Ordering
$c< :: Order -> Order -> Bool
< :: Order -> Order -> Bool
$c<= :: Order -> Order -> Bool
<= :: Order -> Order -> Bool
$c> :: Order -> Order -> Bool
> :: Order -> Order -> Bool
$c>= :: Order -> Order -> Bool
>= :: Order -> Order -> Bool
$cmax :: Order -> Order -> Order
max :: Order -> Order -> Order
$cmin :: Order -> Order -> Order
min :: Order -> Order -> Order
Ord)

instance Arbitrary Order where
   arbitrary :: Gen Order
arbitrary = (Integer -> Order) -> Gen Integer -> Gen Order
forall a b. (a -> b) -> Gen a -> Gen b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Integer -> Order
Order (Integer -> Order) -> (Integer -> Integer) -> Integer -> Order
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer
1Integer -> Integer -> Integer
forall a. C a => a -> a -> a
+) (Integer -> Integer) -> (Integer -> Integer) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Integer
forall a. C a => a -> a
abs) Gen Integer
forall a. Arbitrary a => Gen a
arbitrary

instance Enum Order where
   succ :: Order -> Order
succ (Order Integer
n) = Integer -> Order
Order (Integer
nInteger -> Integer -> Integer
forall a. C a => a -> a -> a
+Integer
1)
   pred :: Order -> Order
pred (Order Integer
n) = Integer -> Order
Order (Integer
nInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
1)
   fromEnum :: Order -> Int
fromEnum (Order Integer
n) = Integer -> Int
forall a. Enum a => a -> Int
fromEnum Integer
n
   toEnum :: Int -> Order
toEnum Int
n = Integer -> Order
Order (Int -> Integer
forall a. Enum a => Int -> a
toEnum Int
n)
   enumFrom :: Order -> [Order]
enumFrom (Order Integer
from) =
      (Integer -> Order) -> [Integer] -> [Order]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> Order
Order ([Integer] -> [Order]) -> [Integer] -> [Order]
forall a b. (a -> b) -> a -> b
$ Integer -> [Integer]
forall a. Enum a => a -> [a]
enumFrom Integer
from
   enumFromThen :: Order -> Order -> [Order]
enumFromThen (Order Integer
from) (Order Integer
thn) =
      (Integer -> Order) -> [Integer] -> [Order]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> Order
Order ([Integer] -> [Order]) -> [Integer] -> [Order]
forall a b. (a -> b) -> a -> b
$ Integer -> Integer -> [Integer]
forall a. Enum a => a -> a -> [a]
enumFromThen Integer
from Integer
thn
   enumFromTo :: Order -> Order -> [Order]
enumFromTo (Order Integer
from) (Order Integer
to) =
      (Integer -> Order) -> [Integer] -> [Order]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> Order
Order ([Integer] -> [Order]) -> [Integer] -> [Order]
forall a b. (a -> b) -> a -> b
$ Integer -> Integer -> [Integer]
forall a. Enum a => a -> a -> [a]
enumFromTo Integer
from Integer
to
   enumFromThenTo :: Order -> Order -> Order -> [Order]
enumFromThenTo (Order Integer
from) (Order Integer
thn) (Order Integer
to) =
      (Integer -> Order) -> [Integer] -> [Order]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> Order
Order ([Integer] -> [Order]) -> [Integer] -> [Order]
forall a b. (a -> b) -> a -> b
$ Integer -> Integer -> Integer -> [Integer]
forall a. Enum a => a -> a -> a -> [a]
enumFromThenTo Integer
from Integer
thn Integer
to

countOrder :: [a] -> Order
countOrder :: forall a. [a] -> Order
countOrder = (Order -> a -> Order) -> Order -> [a] -> Order
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (\Order
o a
_ -> Order -> Order
forall a. Enum a => a -> a
succ Order
o) (Integer -> Order
Order Integer
0)

dividesOrder :: Order -> Order -> Bool
dividesOrder :: Order -> Order -> Bool
dividesOrder (Order Integer
k) (Order Integer
n) =
   Integer -> Integer -> Bool
forall a. (C a, C a) => a -> a -> Bool
divides Integer
k Integer
n


-- class Integral.C a => PrimitiveRoot a where
class PID.C a => PrimitiveRoot a where
   primitiveRootCandidates :: a -> [a]
   maximumOrderOfPrimitiveRootsOfUnity :: a -> Order

instance PrimitiveRoot Integer where
   primitiveRootCandidates :: Integer -> [Integer]
primitiveRootCandidates Integer
modu = [Integer
1 .. Integer
moduInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
1]
   maximumOrderOfPrimitiveRootsOfUnity :: Integer -> Order
maximumOrderOfPrimitiveRootsOfUnity =
      Integer -> Order
maximumOrderOfPrimitiveRootsOfUnityInteger

{-
For 'ordersOfPrimitiveRootsOfUnityInteger'
and the connection to Euler's totient function
we have chosen to have

> primitiveRootsOfUnity m 1 == [1].
-}
primitiveRootsOfUnity ::
   (PrimitiveRoot a, Eq a) => a -> Order -> [a]
primitiveRootsOfUnity :: forall a. (PrimitiveRoot a, Eq a) => a -> Order -> [a]
primitiveRootsOfUnity =
   a -> Order -> [a]
forall a. (PrimitiveRoot a, Eq a) => a -> Order -> [a]
primitiveRootsOfUnityPower

{-
First check, that element x is a root of unity.
If x is not primitive,
this means there is a non-maximal exponent y with x^y=1.
This y must be a divisor of order.
Thus it is enough to check all possibilities of order/q as exponents,
where q is a prime divisor of order.
Computing a single power is much faster
than computing all powers up to the maximum power.

Verifying that a ring has no primitive root of the wanted order
takes a long time if we do it by exhaustive search.
In the case of a=Integer we could first check,
whether the considered residue ring has a primitive root of wanted order
using the Carmichael function.
We could certainly count the number of primitive roots
and stop searching if we reach that count.
-}
primitiveRootsOfUnityPower ::
   (PrimitiveRoot a, Eq a) => a -> Order -> [a]
primitiveRootsOfUnityPower :: forall a. (PrimitiveRoot a, Eq a) => a -> Order -> [a]
primitiveRootsOfUnityPower a
modu (Order Integer
order) =
   let greatDivisors :: [Integer]
greatDivisors = (Integer -> Integer) -> [Integer] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Integer -> Integer -> Integer
forall a. C a => a -> a -> a
div Integer
order) ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$ Integer -> [Integer]
forall a. (C a, Bits a, C a, Ord a) => a -> [a]
uniquePrimeFactors Integer
order
   in  (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
filter
          (\a
n ->
             let pow :: Integer -> a
pow Integer
y = T a -> a
forall a. T a -> a
RC.representative (T a -> a) -> T a -> a
forall a b. (a -> b) -> a -> b
$ (a
n a -> a -> T a
forall a. C a => a -> a -> T a
/: a
modu) T a -> Integer -> T a
forall a. C a => a -> Integer -> a
^ Integer
y
             in  a -> a -> Bool
forall a. C a => a -> a -> Bool
PID.coprime a
n a
modu
                 Bool -> Bool -> Bool
&&
                 Integer -> a
pow Integer
order a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
forall a. C a => a
one
                 Bool -> Bool -> Bool
&&
                 (Integer -> Bool) -> [Integer] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Integer
y -> Integer -> a
pow Integer
y a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
forall a. C a => a
one) [Integer]
greatDivisors) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$
       a -> [a]
forall a. PrimitiveRoot a => a -> [a]
primitiveRootCandidates a
modu

primitiveRootsOfUnityNaive ::
   (PrimitiveRoot a, Eq a) => a -> Order -> [a]
primitiveRootsOfUnityNaive :: forall a. (PrimitiveRoot a, Eq a) => a -> Order -> [a]
primitiveRootsOfUnityNaive a
_ (Order Integer
0) = []
primitiveRootsOfUnityNaive a
modu (Order Integer
expo) =
   (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
filter
      (\a
n ->
         let ([a]
prefix,a
end:[a]
_) =
                Integer -> [a] -> ([a], [a])
forall i a. Integral i => i -> [a] -> ([a], [a])
genericSplitAt (Integer
expoInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
1) ([a] -> ([a], [a])) -> [a] -> ([a], [a])
forall a b. (a -> b) -> a -> b
$ T a -> [a]
forall y. T y -> [y]
SigS.toList (T a -> [a]) -> T a -> [a]
forall a b. (a -> b) -> a -> b
$ a -> a -> T a
forall a. C a => a -> a -> T a
orbit a
modu a
n
         in  (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (a
1a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/=) [a]
prefix Bool -> Bool -> Bool
&& a
enda -> a -> Bool
forall a. Eq a => a -> a -> Bool
==a
1) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$
   a -> [a]
forall a. PrimitiveRoot a => a -> [a]
primitiveRootCandidates a
modu

orbitSet :: Ord a => SigS.T a -> Set.Set a
orbitSet :: forall a. Ord a => T a -> Set a
orbitSet T a
list =
   (a -> (Set a -> Set a) -> Set a -> Set a)
-> (Set a -> Set a) -> T a -> Set a -> Set a
forall x acc. (x -> acc -> acc) -> acc -> T x -> acc
SigS.foldR
      (\a
new Set a -> Set a
cont Set a
seen ->
         if a -> Set a -> Bool
forall a. Ord a => a -> Set a -> Bool
Set.member a
new Set a
seen
           then Set a
seen
           else Set a -> Set a
cont (a -> Set a -> Set a
forall a. Ord a => a -> Set a -> Set a
Set.insert a
new Set a
seen))
      Set a -> Set a
forall a. a -> a
id T a
list Set a
forall a. Set a
Set.empty

orbit :: (Integral.C a) => a -> a -> SigS.T a
orbit :: forall a. C a => a -> a -> T a
orbit a
p a
x = (a -> a) -> a -> T a
forall a. (a -> a) -> a -> T a
SigS.iterate (\a
y -> a -> a -> a
forall a. C a => a -> a -> a
mod (a
xa -> a -> a
forall a. C a => a -> a -> a
*a
y) a
p) a
x


{- |
Does not emit values in ascending order
and may return duplicates (e.g. primitiveRootsOfUnityFullOrbit 70000 10).
I hoped it would be faster than the other implementations
since it eliminates non-roots in large blocks.
However it seems that managing the root candidates in a Set
reduces performance significantly.

The idea:
Start with a seed that is a unit.
Compute its orbit until a one is reached.
Since it is a unit, we always encounter a one.
We do not need to check for non-unit seeds,
since (gcd modu seed) will be a divisor of all seed powers.
In the orbit all numbers are powers of each other.
Now finding the roots is a matter of solving
a Diophantine equation of the exponents.
In one such step all powers in an orbit are classified as roots or non-roots
and thus we can remove them all from the set of root candidates
and continue with the remaining candidates.
Duplicates can occur if a seed
in a later iteration is found again as power of another seed.
-}
primitiveRootsOfUnityFullOrbit ::
   (PrimitiveRoot a, Ord a) => a -> Order -> [a]
primitiveRootsOfUnityFullOrbit :: forall a. (PrimitiveRoot a, Ord a) => a -> Order -> [a]
primitiveRootsOfUnityFullOrbit a
modu Order
expo =
   let search :: Set a -> Maybe ([a], Set a)
search Set a
candidates =
          (((a, Set a) -> ([a], Set a))
 -> Maybe (a, Set a) -> Maybe ([a], Set a))
-> Maybe (a, Set a)
-> ((a, Set a) -> ([a], Set a))
-> Maybe ([a], Set a)
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((a, Set a) -> ([a], Set a))
-> Maybe (a, Set a) -> Maybe ([a], Set a)
forall a b. (a -> b) -> Maybe a -> Maybe b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Set a -> Maybe (a, Set a)
forall a. Set a -> Maybe (a, Set a)
Set.minView Set a
candidates) (((a, Set a) -> ([a], Set a)) -> Maybe ([a], Set a))
-> ((a, Set a) -> ([a], Set a)) -> Maybe ([a], Set a)
forall a b. (a -> b) -> a -> b
$ \(a
x,Set a
rest) ->
          ([a] -> Set a) -> ([a], [a]) -> ([a], Set a)
forall b c a. (b -> c) -> (a, b) -> (a, c)
mapSnd (Set a -> Set a -> Set a
forall a. Ord a => Set a -> Set a -> Set a
Set.difference Set a
rest (Set a -> Set a) -> ([a] -> Set a) -> [a] -> Set a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> Set a
forall a. Ord a => [a] -> Set a
Set.fromList) (([a], [a]) -> ([a], Set a)) -> ([a], [a]) -> ([a], Set a)
forall a b. (a -> b) -> a -> b
$
          a -> Order -> a -> ([a], [a])
forall a. (PrimitiveRoot a, Ord a) => a -> Order -> a -> ([a], [a])
primitiveRootsOfOrbit a
modu Order
expo a
x
   in  [[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([[a]] -> [a]) -> [[a]] -> [a]
forall a b. (a -> b) -> a -> b
$ (Set a -> Maybe ([a], Set a)) -> Set a -> [[a]]
forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr Set a -> Maybe ([a], Set a)
search (Set a -> [[a]]) -> Set a -> [[a]]
forall a b. (a -> b) -> a -> b
$ [a] -> Set a
forall a. Ord a => [a] -> Set a
Set.fromList ([a] -> Set a) -> [a] -> Set a
forall a b. (a -> b) -> a -> b
$
       -- needed for modules with many divisors
       (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
filter (a -> a -> Bool
forall a. C a => a -> a -> Bool
PID.coprime a
modu) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$
       a -> [a]
forall a. PrimitiveRoot a => a -> [a]
primitiveRootCandidates a
modu

primitiveRootsOfUnityFullOrbitTest ::
   (PrimitiveRoot a, Ord a) => a -> Order -> [(a,[a])]
primitiveRootsOfUnityFullOrbitTest :: forall a. (PrimitiveRoot a, Ord a) => a -> Order -> [(a, [a])]
primitiveRootsOfUnityFullOrbitTest a
modu Order
expo =
   let search :: Set a -> Maybe ((a, [a]), Set a)
search Set a
candidates =
          (((a, Set a) -> ((a, [a]), Set a))
 -> Maybe (a, Set a) -> Maybe ((a, [a]), Set a))
-> Maybe (a, Set a)
-> ((a, Set a) -> ((a, [a]), Set a))
-> Maybe ((a, [a]), Set a)
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((a, Set a) -> ((a, [a]), Set a))
-> Maybe (a, Set a) -> Maybe ((a, [a]), Set a)
forall a b. (a -> b) -> Maybe a -> Maybe b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Set a -> Maybe (a, Set a)
forall a. Set a -> Maybe (a, Set a)
Set.minView Set a
candidates) (((a, Set a) -> ((a, [a]), Set a)) -> Maybe ((a, [a]), Set a))
-> ((a, Set a) -> ((a, [a]), Set a)) -> Maybe ((a, [a]), Set a)
forall a b. (a -> b) -> a -> b
$ \(a
x,Set a
rest) ->
          ([a] -> (a, [a]), [a] -> Set a) -> ([a], [a]) -> ((a, [a]), Set a)
forall a c b d. (a -> c, b -> d) -> (a, b) -> (c, d)
mapPair ((,) a
x,
                   Set a -> Set a -> Set a
forall a. Ord a => Set a -> Set a -> Set a
Set.difference Set a
rest (Set a -> Set a) -> ([a] -> Set a) -> [a] -> Set a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> Set a
forall a. Ord a => [a] -> Set a
Set.fromList) (([a], [a]) -> ((a, [a]), Set a))
-> ([a], [a]) -> ((a, [a]), Set a)
forall a b. (a -> b) -> a -> b
$
          a -> Order -> a -> ([a], [a])
forall a. (PrimitiveRoot a, Ord a) => a -> Order -> a -> ([a], [a])
primitiveRootsOfOrbit a
modu Order
expo a
x
   in  (Set a -> Maybe ((a, [a]), Set a)) -> Set a -> [(a, [a])]
forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr Set a -> Maybe ((a, [a]), Set a)
search (Set a -> [(a, [a])]) -> Set a -> [(a, [a])]
forall a b. (a -> b) -> a -> b
$ [a] -> Set a
forall a. Ord a => [a] -> Set a
Set.fromList ([a] -> Set a) -> [a] -> Set a
forall a b. (a -> b) -> a -> b
$
       (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
filter (a -> a -> Bool
forall a. C a => a -> a -> Bool
PID.coprime a
modu) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$
       a -> [a]
forall a. PrimitiveRoot a => a -> [a]
primitiveRootCandidates a
modu

primitiveRootsOfOrbit ::
   (PrimitiveRoot a, Ord a) => a -> Order -> a -> ([a], [a])
primitiveRootsOfOrbit :: forall a. (PrimitiveRoot a, Ord a) => a -> Order -> a -> ([a], [a])
primitiveRootsOfOrbit a
modu (Order Integer
expo) a
x =
   let orb :: [a]
orb = (a
1a -> [a] -> [a]
forall a. a -> [a] -> [a]
:) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (a
1a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/=) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (\a
y -> a -> a -> a
forall a. C a => a -> a -> a
mod (a
xa -> a -> a
forall a. C a => a -> a -> a
*a
y) a
modu) a
x
       (Order Integer
orbitSize) = [a] -> Order
forall a. [a] -> Order
countOrder [a]
orb
   in  (if Integer
expoInteger -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==Integer
0
          then []
          else
            {-
            size = length orb
            Search for m and k with 0<k and 0<m and m<size
            and expo*m = size*k
            such that for all l with 0<l and l<k
            m does not divide size*l.
            To this end we ask for every m
            for the smallest r such that size divides r*m.
            If r=expo then x^m is a primitive root of order expo.
            If r<expo then x^m has order smaller than expo.
            The searched r is div size (gcd size m).
            However expo = div size (gcd size m) implies,
            that expo is a divisor of size.
                expo = div size (gcd size m)
             => div size expo = gcd size m
                s = gcd size m
            We have to consider for m
            only the multiples of s.
            Then divide both sides of the equation by s, yielding
                1 = gcd expo m'
            -}
            case Integer -> Integer -> (Integer, Integer)
forall a. C a => a -> a -> (a, a)
divMod Integer
orbitSize Integer
expo of
               (Integer
s,Integer
0) ->
                  ((Integer, a) -> a) -> [(Integer, a)] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (Integer, a) -> a
forall a b. (a, b) -> b
snd ([(Integer, a)] -> [a]) -> [(Integer, a)] -> [a]
forall a b. (a -> b) -> a -> b
$ ((Integer, a) -> Bool) -> [(Integer, a)] -> [(Integer, a)]
forall a. (a -> Bool) -> [a] -> [a]
filter (Integer -> Integer -> Bool
forall a. C a => a -> a -> Bool
PID.coprime Integer
expo (Integer -> Bool)
-> ((Integer, a) -> Integer) -> (Integer, a) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer, a) -> Integer
forall a b. (a, b) -> a
fst) ([(Integer, a)] -> [(Integer, a)])
-> [(Integer, a)] -> [(Integer, a)]
forall a b. (a -> b) -> a -> b
$
                  [Integer] -> [a] -> [(Integer, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip
                     [Integer
0 .. Integer
expoInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
1]
                     -- (ListHT.sieve s $ orb)
                     (([a] -> a) -> [[a]] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map [a] -> a
forall a. HasCallStack => [a] -> a
head ([[a]] -> [a]) -> [[a]] -> [a]
forall a b. (a -> b) -> a -> b
$ ([a] -> [a]) -> [a] -> [[a]]
forall a. (a -> a) -> a -> [a]
iterate (Integer -> [a] -> [a]
forall i a. Integral i => i -> [a] -> [a]
genericDrop Integer
s) [a]
orb)
               (Integer, Integer)
_ -> [],
        [a]
orb)


hasPrimitiveRootOfUnityNaive ::
   (PrimitiveRoot a, Ord a) => a -> Order -> Bool
hasPrimitiveRootOfUnityNaive :: forall a. (PrimitiveRoot a, Ord a) => a -> Order -> Bool
hasPrimitiveRootOfUnityNaive a
modu Order
expo =
   ((a, Order) -> Bool) -> [(a, Order)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Order -> Order -> Bool
dividesOrder Order
expo (Order -> Bool) -> ((a, Order) -> Order) -> (a, Order) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a, Order) -> Order
forall a b. (a, b) -> b
snd) ([(a, Order)] -> Bool) -> [(a, Order)] -> Bool
forall a b. (a -> b) -> a -> b
$
   a -> [(a, Order)]
forall a. (PrimitiveRoot a, Ord a) => a -> [(a, Order)]
ordersOfPrimitiveRootsOfUnityTest a
modu

{-
This should be a maximum both with respect to magnitude and to divisibility.
-}
maximumOrderOfPrimitiveRootsOfUnityNaive ::
   (PrimitiveRoot a, Ord a) => a -> Order
maximumOrderOfPrimitiveRootsOfUnityNaive :: forall a. (PrimitiveRoot a, Ord a) => a -> Order
maximumOrderOfPrimitiveRootsOfUnityNaive =
   (Order -> Order -> Order) -> Order -> [Order] -> Order
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl Order -> Order -> Order
forall a. Ord a => a -> a -> a
max (Integer -> Order
Order Integer
1) ([Order] -> Order) -> (a -> [Order]) -> a -> Order
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((a, Order) -> Order) -> [(a, Order)] -> [Order]
forall a b. (a -> b) -> [a] -> [b]
map (a, Order) -> Order
forall a b. (a, b) -> b
snd ([(a, Order)] -> [Order]) -> (a -> [(a, Order)]) -> a -> [Order]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> [(a, Order)]
forall a. (PrimitiveRoot a, Ord a) => a -> [(a, Order)]
ordersOfPrimitiveRootsOfUnityTest

{- |
Computes a list of seeds and according maximum orders of roots of unity.
All divisors of those maximum orders are possible orders of roots of unity, too.
-}
ordersOfPrimitiveRootsOfUnityTest ::
   (PrimitiveRoot a, Ord a) => a -> [(a, Order)]
ordersOfPrimitiveRootsOfUnityTest :: forall a. (PrimitiveRoot a, Ord a) => a -> [(a, Order)]
ordersOfPrimitiveRootsOfUnityTest a
modu =
   let search :: Set a -> Maybe ((a, Order), Set a)
search Set a
candidates =
          (((a, Set a) -> ((a, Order), Set a))
 -> Maybe (a, Set a) -> Maybe ((a, Order), Set a))
-> Maybe (a, Set a)
-> ((a, Set a) -> ((a, Order), Set a))
-> Maybe ((a, Order), Set a)
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((a, Set a) -> ((a, Order), Set a))
-> Maybe (a, Set a) -> Maybe ((a, Order), Set a)
forall a b. (a -> b) -> Maybe a -> Maybe b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Set a -> Maybe (a, Set a)
forall a. Set a -> Maybe (a, Set a)
Set.minView Set a
candidates) (((a, Set a) -> ((a, Order), Set a)) -> Maybe ((a, Order), Set a))
-> ((a, Set a) -> ((a, Order), Set a)) -> Maybe ((a, Order), Set a)
forall a b. (a -> b) -> a -> b
$ \(a
x,Set a
rest) ->
          (Order -> (a, Order), [a] -> Set a)
-> (Order, [a]) -> ((a, Order), Set a)
forall a c b d. (a -> c, b -> d) -> (a, b) -> (c, d)
mapPair ((,) a
x,
                   Set a -> Set a -> Set a
forall a. Ord a => Set a -> Set a -> Set a
Set.difference Set a
rest (Set a -> Set a) -> ([a] -> Set a) -> [a] -> Set a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> Set a
forall a. Ord a => [a] -> Set a
Set.fromList) ((Order, [a]) -> ((a, Order), Set a))
-> (Order, [a]) -> ((a, Order), Set a)
forall a b. (a -> b) -> a -> b
$
          a -> a -> (Order, [a])
forall a. (PrimitiveRoot a, Ord a) => a -> a -> (Order, [a])
orderOfOrbit a
modu a
x
   in  (Set a -> Maybe ((a, Order), Set a)) -> Set a -> [(a, Order)]
forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr Set a -> Maybe ((a, Order), Set a)
search (Set a -> [(a, Order)]) -> Set a -> [(a, Order)]
forall a b. (a -> b) -> a -> b
$ [a] -> Set a
forall a. Ord a => [a] -> Set a
Set.fromList ([a] -> Set a) -> [a] -> Set a
forall a b. (a -> b) -> a -> b
$
       (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
filter (a -> a -> Bool
forall a. C a => a -> a -> Bool
PID.coprime a
modu) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$
       a -> [a]
forall a. PrimitiveRoot a => a -> [a]
primitiveRootCandidates a
modu

{- |
modu and x must be coprime.
If they are not,
then none of the numbers in the orbit is a root of unity
and the function enters an infinite loop.
-}
orderOfOrbit ::
   (PrimitiveRoot a, Ord a) => a -> a -> (Order, [a])
orderOfOrbit :: forall a. (PrimitiveRoot a, Ord a) => a -> a -> (Order, [a])
orderOfOrbit a
modu a
x =
   let cyc :: [a]
cyc = (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (a
forall a. C a => a
onea -> a -> Bool
forall a. Eq a => a -> a -> Bool
/=) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ T a -> [a]
forall y. T y -> [y]
SigS.toList (T a -> [a]) -> T a -> [a]
forall a b. (a -> b) -> a -> b
$ a -> a -> T a
forall a. C a => a -> a -> T a
orbit a
modu a
x
   in  (Order -> Order
forall a. Enum a => a -> a
succ (Order -> Order) -> Order -> Order
forall a b. (a -> b) -> a -> b
$ [a] -> Order
forall a. [a] -> Order
countOrder [a]
cyc, [a]
cyc)


{-
This test speeds up 'hasPrimitiveRootOfUnityNaive' considerably
by considering the prime factors of modu.
If modu is a prime, then the ring has a multiplicative generator,
i.e. a primitive root of unity of order modu-1.
That is, all primitive roots of unity are of an order that divides modu-1.
It seems that if modu is a power of a prime,
then the according ring has also multiplicative generator.
I think this is the reason for generalising the Rader transform
to signals of prime power length.
-}
hasPrimitiveRootOfUnityInteger ::
   Integer -> Order -> Bool
hasPrimitiveRootOfUnityInteger :: Integer -> Order -> Bool
hasPrimitiveRootOfUnityInteger Integer
modu Order
expo =
   Order -> Order -> Bool
dividesOrder Order
expo (Order -> Bool) -> Order -> Bool
forall a b. (a -> b) -> a -> b
$
   Integer -> Order
maximumOrderOfPrimitiveRootsOfUnityInteger Integer
modu

{-
Carmichael theorem:
If the integer residue rings with coprime moduli m0 and m1
have primitive roots of maximum order o0 and o1, respectively,
then the integer ring with modulus m0*m1
has maximum order (lcm o0 o1).
-}

{-
This is the Carmichael function.
<http://oeis.org/A002322>
-}
maximumOrderOfPrimitiveRootsOfUnityInteger ::
   Integer -> Order
maximumOrderOfPrimitiveRootsOfUnityInteger :: Integer -> Order
maximumOrderOfPrimitiveRootsOfUnityInteger =
   Integer -> Order
Order (Integer -> Order) -> (Integer -> Integer) -> Integer -> Order
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   [Integer] -> Integer
forall a. C a => [a] -> a
lcmMulti ([Integer] -> Integer)
-> (Integer -> [Integer]) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   ((Integer, Integer) -> Integer)
-> [(Integer, Integer)] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map
      (\(Integer
e,Integer
p) ->
         if Integer
p Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
2 Bool -> Bool -> Bool
&& Integer
e Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
2
           then Integer
pInteger -> Integer -> Integer
forall a. C a => a -> Integer -> a
^(Integer
eInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
2)
           else Integer
pInteger -> Integer -> Integer
forall a. C a => a -> Integer -> a
^(Integer
eInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
1) Integer -> Integer -> Integer
forall a. C a => a -> a -> a
* (Integer
pInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
1)) ([(Integer, Integer)] -> [Integer])
-> (Integer -> [(Integer, Integer)]) -> Integer -> [Integer]
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   ((Int, Integer) -> (Integer, Integer))
-> [(Int, Integer)] -> [(Integer, Integer)]
forall a b. (a -> b) -> [a] -> [b]
map ((Int -> Integer) -> (Int, Integer) -> (Integer, Integer)
forall a c b. (a -> c) -> (a, b) -> (c, b)
mapFst Int -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral) ([(Int, Integer)] -> [(Integer, Integer)])
-> (Integer -> [(Int, Integer)]) -> Integer -> [(Integer, Integer)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   Integer -> [(Int, Integer)]
forall a. (PrimitiveRoot a, Ord a) => a -> [(Int, a)]
primeFactors


{-
The sum of the sub-lists should equal the Euler totient function values
<http://oeis.org/A000010>.
-}
ordersOfPrimitiveRootsOfUnityInteger :: [[Int]]
ordersOfPrimitiveRootsOfUnityInteger :: [[Int]]
ordersOfPrimitiveRootsOfUnityInteger =
   ((Integer -> [Int]) -> [Integer] -> [[Int]])
-> [Integer] -> (Integer -> [Int]) -> [[Int]]
forall a b c. (a -> b -> c) -> b -> a -> c
flip (Integer -> [Int]) -> [Integer] -> [[Int]]
forall a b. (a -> b) -> [a] -> [b]
map [Integer
1..] ((Integer -> [Int]) -> [[Int]]) -> (Integer -> [Int]) -> [[Int]]
forall a b. (a -> b) -> a -> b
$ \Integer
modu ->
   let maxOrder :: Order
maxOrder = Integer -> Order
forall a. PrimitiveRoot a => a -> Order
maximumOrderOfPrimitiveRootsOfUnity (Integer
modu::Integer)
   in  (Order -> Int) -> [Order] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map ([Integer] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Integer] -> Int) -> (Order -> [Integer]) -> Order -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Order -> [Integer]
forall a. (PrimitiveRoot a, Eq a) => a -> Order -> [a]
primitiveRootsOfUnityPower Integer
modu) ([Order] -> [Int]) -> [Order] -> [Int]
forall a b. (a -> b) -> a -> b
$
--       filter (flip divides maxOrder) $
       [Integer -> Order
Order Integer
1 .. Order
maxOrder]

ordersOfRootsOfUnityInteger :: [[Int]]
ordersOfRootsOfUnityInteger :: [[Int]]
ordersOfRootsOfUnityInteger =
   ((Integer -> [Int]) -> [Integer] -> [[Int]])
-> [Integer] -> (Integer -> [Int]) -> [[Int]]
forall a b c. (a -> b -> c) -> b -> a -> c
flip (Integer -> [Int]) -> [Integer] -> [[Int]]
forall a b. (a -> b) -> [a] -> [b]
map [Integer
1..] ((Integer -> [Int]) -> [[Int]]) -> (Integer -> [Int]) -> [[Int]]
forall a b. (a -> b) -> a -> b
$ \Integer
modu ->
   (Order -> Int) -> [Order] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map ([Integer] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Integer] -> Int) -> (Order -> [Integer]) -> Order -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Order -> [Integer]
forall a. (PrimitiveRoot a, Eq a) => a -> Order -> [a]
rootsOfUnityPower (Integer
modu::Integer)) ([Order] -> [Int]) -> [Order] -> [Int]
forall a b. (a -> b) -> a -> b
$
   [Integer -> Order
Order Integer
1 ..]
{-
mapM_ print $ map (\n -> (n, maximumOrderOfPrimitiveRootsOfUnityInteger (fromIntegral n), take 30 $ ordersOfRootsOfUnityInteger !! (n-1))) [2..30]

mapM_ print $ map (\n -> (n, maximumOrderOfPrimitiveRootsOfUnityInteger (fromIntegral n), let row = ordersOfRootsOfUnityInteger !! (n-1) in map (row!!) $ map pred $ take 10 $ iterate (2*) 1)) [2..30]
-}

ordersOfRootsOfUnityIntegerCondensed :: [[Int]]
ordersOfRootsOfUnityIntegerCondensed :: [[Int]]
ordersOfRootsOfUnityIntegerCondensed =
   ((Integer -> [Int]) -> [Integer] -> [[Int]])
-> [Integer] -> (Integer -> [Int]) -> [[Int]]
forall a b c. (a -> b -> c) -> b -> a -> c
flip (Integer -> [Int]) -> [Integer] -> [[Int]]
forall a b. (a -> b) -> [a] -> [b]
map [Integer
1..] ((Integer -> [Int]) -> [[Int]]) -> (Integer -> [Int]) -> [[Int]]
forall a b. (a -> b) -> a -> b
$ \Integer
modu ->
   let maxOrder :: Order
maxOrder = Integer -> Order
forall a. PrimitiveRoot a => a -> Order
maximumOrderOfPrimitiveRootsOfUnity (Integer
modu::Integer)
   in  (Order -> Int) -> [Order] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map ([Integer] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Integer] -> Int) -> (Order -> [Integer]) -> Order -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Order -> [Integer]
forall a. (PrimitiveRoot a, Eq a) => a -> Order -> [a]
rootsOfUnityPower Integer
modu) ([Order] -> [Int]) -> [Order] -> [Int]
forall a b. (a -> b) -> a -> b
$
--       filter (flip divides maxOrder) $
       [Integer -> Order
Order Integer
1 .. Order
maxOrder]

rootsOfUnityPower ::
   (PrimitiveRoot a, Eq a) => a -> Order -> [a]
rootsOfUnityPower :: forall a. (PrimitiveRoot a, Eq a) => a -> Order -> [a]
rootsOfUnityPower a
modu (Order Integer
expo) =
   (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
filter
      (\a
n ->
         a -> a -> Bool
forall a. C a => a -> a -> Bool
PID.coprime a
n a
modu
         Bool -> Bool -> Bool
&&
         T a -> a
forall a. T a -> a
RC.representative ((a
n a -> a -> T a
forall a. C a => a -> a -> T a
/: a
modu) T a -> Integer -> T a
forall a. C a => a -> Integer -> a
^ Integer
expo) a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
forall a. C a => a
one) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$
   a -> [a]
forall a. PrimitiveRoot a => a -> [a]
primitiveRootCandidates a
modu

{-
Corollary from the Carmichael function properties:
If in Z_n there is a primitive root r of unity of order o,
then for every Z_{m \cdot n} there is also a primitive root of unity
with the same order.

Collary:
If in Z_n1 you have a primitive root of order o1,
and so on for Z_{n_k} and order ok,
then Z_{lcm(n1,...,nk)} has primitive roots for every of the order o1,...,on.

Conjecture:
If Z_n has a total number of m primitive roots of unity of order o,
then Z_{k*n} has at least m primitive roots of unity of order o.

Conjecture:
Precondition: In Z_n there is a primitive root of prime order o.
Claims:
a) There are natural numbers m and k with n = m*(k*o+1) or n = m*o.
b) The smallest such n is of the form k*o+1 with k>1.

Counterexample for a) and non-prime o: o=15, n=77
Counterexample for b) and non-prime o:
   o=20, n=25; o=27, n=81; o=54, n=81; o=55, n=121

Corollary from definition of Carmichael function:
For n>1, Z_{2^{n+2}} has primitive root of unity of order 2^n.
-}

{- |
Given an order find integer residue rings
where roots of unity of this order exist.
The way they are constructed also warrants,
that 'order' is a unit (i.e. invertible) in those rings.

The list is not exhaustive
but computes suggestions quickly.
The first found modulus is often the smallest one that exist,
but not always (smallest counter-example: Order 80).
However, the first modulus is not the smallest one
among the ones that only have the wanted primitive root,
but where 'order' is not necessarily a unit.
E.g.

> ringsWithPrimitiveRootOfUnityAndUnit 840 == 2521 : 3361 : ...

but the smallest modulus is 1763.

Most of the numbers are primes.
For these the recursion property of the Carmichael function
immediately yields that there are primitive roots of unity of order 'order'.
-}
ringsWithPrimitiveRootOfUnityAndUnit :: Order -> [Integer]
ringsWithPrimitiveRootOfUnityAndUnit :: Order -> [Integer]
ringsWithPrimitiveRootOfUnityAndUnit order :: Order
order@(Order Integer
k) =
   (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
filter ((Integer -> Order -> Bool) -> Order -> Integer -> Bool
forall a b c. (a -> b -> c) -> b -> a -> c
flip Integer -> Order -> Bool
hasPrimitiveRootOfUnityInteger Order
order) ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$
   (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer
kInteger -> Integer -> Integer
forall a. C a => a -> a -> a
+) Integer
1


ringsWithPrimitiveRootsOfUnityAndUnitsNaive :: [Order] -> [Integer] -> [Integer]
ringsWithPrimitiveRootsOfUnityAndUnitsNaive :: [Order] -> [Integer] -> [Integer]
ringsWithPrimitiveRootsOfUnityAndUnitsNaive [Order]
rootOrders [Integer]
units =
   (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
filter
      (\Integer
n ->
         (Order -> Bool) -> [Order] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Integer -> Order -> Bool
hasPrimitiveRootOfUnityInteger Integer
n) [Order]
rootOrders Bool -> Bool -> Bool
&&
         (Integer -> Bool) -> [Integer] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Integer -> Integer -> Bool
forall a. C a => a -> a -> Bool
PID.coprime Integer
n) [Integer]
units)
      [Integer
1..]


{-
It would be nice to have the Omega monad here
in order to enumerate all rings.
-}
ringWithPrimitiveRootsOfUnityAndUnits :: [Order] -> [Integer] -> Integer
ringWithPrimitiveRootsOfUnityAndUnits :: [Order] -> [Integer] -> Integer
ringWithPrimitiveRootsOfUnityAndUnits [Order]
rootOrders [Integer]
units =
   let p :: Integer
p = [Integer] -> Integer
forall a. C a => [a] -> a
lcmMulti [Integer]
units
   in  [Integer] -> Integer
forall a. C a => [a] -> a
lcmMulti ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$
       (Order -> Integer) -> [Order] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map ([Integer] -> Integer
forall a. HasCallStack => [a] -> a
head ([Integer] -> Integer) -> (Order -> [Integer]) -> Order -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
filter (Integer -> Integer -> Bool
forall a. C a => a -> a -> Bool
PID.coprime Integer
p) ([Integer] -> [Integer])
-> (Order -> [Integer]) -> Order -> [Integer]
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
            Order -> [Integer]
ringsWithPrimitiveRootOfUnityAndUnit) ([Order] -> [Integer]) -> [Order] -> [Integer]
forall a b. (a -> b) -> a -> b
$
       [Order]
rootOrders

{-
Search for an appriopriate modulus
using the recursive definition of the Carmichael function.
We chose the prime factors of the Carmichael function argument
such that we get at least the prime factors in the function value that we need.

The modulus constructed this way is usually not the smallest possible
and it also does not provide that 'n' is a unit in the residue ring.
The simpler function 'ringsWithPrimitiveRootOfUnityAndUnit'
will usually produce a smaller modulus.
-}
ringWithPrimitiveRootsOfUnity :: Order -> Integer
ringWithPrimitiveRootsOfUnity :: Order -> Integer
ringWithPrimitiveRootsOfUnity (Order Integer
n) =
   case Integer
n of
      Integer
0 -> Integer
2
      Integer
_ ->
         [Integer] -> Integer
forall a. C a => [a] -> a
product ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$ ((Int, Integer) -> Integer) -> [(Int, Integer)] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map ((Int -> Integer -> Integer) -> (Int, Integer) -> Integer
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Int -> Integer -> Integer
forall a b. (C a, C b) => b -> a -> a
ringPower) ([(Int, Integer)] -> [Integer]) -> [(Int, Integer)] -> [Integer]
forall a b. (a -> b) -> a -> b
$ (Map Integer Int, [(Int, Integer)]) -> [(Int, Integer)]
forall a b. (a, b) -> b
snd ((Map Integer Int, [(Int, Integer)]) -> [(Int, Integer)])
-> (Map Integer Int, [(Int, Integer)]) -> [(Int, Integer)]
forall a b. (a -> b) -> a -> b
$
         (Map Integer Int
 -> (Int, Integer) -> (Map Integer Int, (Int, Integer)))
-> Map Integer Int
-> [(Int, Integer)]
-> (Map Integer Int, [(Int, Integer)])
forall (t :: * -> *) s a b.
Traversable t =>
(s -> a -> (s, b)) -> s -> t a -> (s, t b)
mapAccumL
            (\Map Integer Int
factors (Int
e,Integer
p) ->
               if Int -> Integer -> Map Integer Int -> Int
forall k a. Ord k => a -> k -> Map k a -> a
Map.findWithDefault Int
0 Integer
p Map Integer Int
factors Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
e
                 then (Map Integer Int
factors, (Int
0,Integer
p))
                 else
                   if Integer
pInteger -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==Integer
2
                     then
                       (Map Integer Int
factors,
                        case Int
e of
                           Int
0 -> (Int
0,Integer
2)
                           Int
1 -> (Int
1,Integer
3)
                           Int
2 -> (Int
1,Integer
5)
                           Int
_ -> (Int
eInt -> Int -> Int
forall a. C a => a -> a -> a
+Int
2, Integer
2))
                     else
                       ((Int -> Int -> Int)
-> Map Integer Int -> Map Integer Int -> Map Integer Int
forall k a. Ord k => (a -> a -> a) -> Map k a -> Map k a -> Map k a
Map.unionWith Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Map Integer Int
factors (Map Integer Int -> Map Integer Int)
-> Map Integer Int -> Map Integer Int
forall a b. (a -> b) -> a -> b
$
                           [(Integer, Int)] -> Map Integer Int
forall k a. Ord k => [(k, a)] -> Map k a
Map.fromList ([(Integer, Int)] -> Map Integer Int)
-> [(Integer, Int)] -> Map Integer Int
forall a b. (a -> b) -> a -> b
$ ((Int, Integer) -> (Integer, Int))
-> [(Int, Integer)] -> [(Integer, Int)]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Integer) -> (Integer, Int)
forall a b. (a, b) -> (b, a)
swap ([(Int, Integer)] -> [(Integer, Int)])
-> [(Int, Integer)] -> [(Integer, Int)]
forall a b. (a -> b) -> a -> b
$ Integer -> [(Int, Integer)]
forall a. (PrimitiveRoot a, Ord a) => a -> [(Int, a)]
primeFactors (Integer -> [(Int, Integer)]) -> Integer -> [(Int, Integer)]
forall a b. (a -> b) -> a -> b
$ Integer
pInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
1,
                        (Int
eInt -> Int -> Int
forall a. C a => a -> a -> a
+Int
1, Integer
p)))
            Map Integer Int
forall k a. Map k a
Map.empty ([(Int, Integer)] -> (Map Integer Int, [(Int, Integer)]))
-> [(Int, Integer)] -> (Map Integer Int, [(Int, Integer)])
forall a b. (a -> b) -> a -> b
$
         [(Int, Integer)] -> [(Int, Integer)]
forall a. [a] -> [a]
reverse ([(Int, Integer)] -> [(Int, Integer)])
-> [(Int, Integer)] -> [(Int, Integer)]
forall a b. (a -> b) -> a -> b
$ Integer -> [(Int, Integer)]
forall a. (PrimitiveRoot a, Ord a) => a -> [(Int, a)]
primeFactors (Integer -> [(Int, Integer)]) -> Integer -> [(Int, Integer)]
forall a b. (a -> b) -> a -> b
$ [Integer] -> Integer
forall a. C a => [a] -> a
lcmMulti ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$
         Integer
n Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: (Integer -> Integer) -> [Integer] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Integer -> Integer -> Integer
forall a. C a => a -> a -> a
subtract Integer
1) (Integer -> [Integer]
partialPrimes Integer
n)

lcmMulti :: (PID.C a) => [a] -> a
lcmMulti :: forall a. C a => [a] -> a
lcmMulti = (a -> a -> a) -> a -> [a] -> a
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl a -> a -> a
forall a. C a => a -> a -> a
lcm a
forall a. C a => a
one


{- |
List all numbers that only contain prime factors 2 and 3 in ascending order.
<http://oeis.org/A003586>
-}
numbers3Smooth :: [Integer]
numbers3Smooth :: [Integer]
numbers3Smooth = [Integer]
numbers3SmoothCorec

numbers3SmoothCorec :: [Integer]
numbers3SmoothCorec :: [Integer]
numbers3SmoothCorec = Integer -> [Integer] -> [Integer]
forall a. (Ord a, C a) => a -> [a] -> [a]
mergePowers Integer
3 ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer
2Integer -> Integer -> Integer
forall a. C a => a -> a -> a
*) Integer
1

mergePowers :: (Ord a, Ring.C a) => a -> [a] -> [a]
mergePowers :: forall a. (Ord a, C a) => a -> [a] -> [a]
mergePowers a
_ [] = []
mergePowers a
p (a
x:[a]
xs) =
   let ys :: [a]
ys = a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> a -> Bool) -> [a] -> [a] -> [a]
forall a. (a -> a -> Bool) -> [a] -> [a] -> [a]
ListHT.mergeBy a -> a -> Bool
forall a. Ord a => a -> a -> Bool
(<=) [a]
xs ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
pa -> a -> a
forall a. C a => a -> a -> a
*) [a]
ys)
   in  [a]
ys

numbers3SmoothFoldr :: [Integer]
numbers3SmoothFoldr :: [Integer]
numbers3SmoothFoldr =
   ([Integer] -> [Integer] -> [Integer])
-> [Integer] -> [[Integer]] -> [Integer]
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr
      (\(Integer
x0:Integer
x1:[Integer]
xs) [Integer]
ys -> Integer
x0 Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: Integer
x1 Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: (Integer -> Integer -> Bool) -> [Integer] -> [Integer] -> [Integer]
forall a. (a -> a -> Bool) -> [a] -> [a] -> [a]
ListHT.mergeBy Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
(<=) [Integer]
xs [Integer]
ys)
      ([Char] -> [Integer]
forall a. HasCallStack => [Char] -> a
error [Char]
"numbers3SmoothFoldr: infinite list should not have an end") ([[Integer]] -> [Integer]) -> [[Integer]] -> [Integer]
forall a b. (a -> b) -> a -> b
$
   ([Integer] -> [Integer]) -> [Integer] -> [[Integer]]
forall a. (a -> a) -> a -> [a]
iterate ((Integer -> Integer) -> [Integer] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Integer
3Integer -> Integer -> Integer
forall a. C a => a -> a -> a
*)) ([Integer] -> [[Integer]]) -> [Integer] -> [[Integer]]
forall a b. (a -> b) -> a -> b
$
   (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer
2Integer -> Integer -> Integer
forall a. C a => a -> a -> a
*) Integer
1

numbers3SmoothSet :: [Integer]
numbers3SmoothSet :: [Integer]
numbers3SmoothSet =
   (Set Integer -> Maybe (Integer, Set Integer))
-> Set Integer -> [Integer]
forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr
      (((Integer, Set Integer) -> (Integer, Set Integer))
-> Maybe (Integer, Set Integer) -> Maybe (Integer, Set Integer)
forall a b. (a -> b) -> Maybe a -> Maybe b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\(Integer
m,Set Integer
rest) -> (Integer
m, Set Integer -> Set Integer -> Set Integer
forall a. Ord a => Set a -> Set a -> Set a
Set.union Set Integer
rest (Set Integer -> Set Integer) -> Set Integer -> Set Integer
forall a b. (a -> b) -> a -> b
$ [Integer] -> Set Integer
forall a. Eq a => [a] -> Set a
Set.fromAscList [Integer
2Integer -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
m,Integer
3Integer -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
m])) (Maybe (Integer, Set Integer) -> Maybe (Integer, Set Integer))
-> (Set Integer -> Maybe (Integer, Set Integer))
-> Set Integer
-> Maybe (Integer, Set Integer)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       Set Integer -> Maybe (Integer, Set Integer)
forall a. Set a -> Maybe (a, Set a)
Set.minView) (Set Integer -> [Integer]) -> Set Integer -> [Integer]
forall a b. (a -> b) -> a -> b
$
   Integer -> Set Integer
forall a. a -> Set a
Set.singleton Integer
1


{-
Hamming sequence
<http://oeis.org/A051037>
-}
numbers5Smooth :: [Integer]
numbers5Smooth :: [Integer]
numbers5Smooth = [Integer]
numbers5SmoothCorec

numbers5SmoothCorec :: [Integer]
numbers5SmoothCorec :: [Integer]
numbers5SmoothCorec =
   if Bool
False
     then -- causes permanent storage of numbers3SmoothCorec
          Integer -> [Integer] -> [Integer]
forall a. (Ord a, C a) => a -> [a] -> [a]
mergePowers Integer
5 ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$ [Integer]
numbers3SmoothCorec
     else Integer -> [Integer] -> [Integer]
forall a. (Ord a, C a) => a -> [a] -> [a]
mergePowers Integer
5 ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$ Integer -> [Integer] -> [Integer]
forall a. (Ord a, C a) => a -> [a] -> [a]
mergePowers Integer
3 ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer
2Integer -> Integer -> Integer
forall a. C a => a -> a -> a
*) Integer
1

numbers5SmoothFoldr :: [Integer]
numbers5SmoothFoldr :: [Integer]
numbers5SmoothFoldr =
   ([Integer] -> [Integer] -> [Integer])
-> [Integer] -> [[Integer]] -> [Integer]
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr
      (\(Integer
x0:Integer
x1:Integer
x2:[Integer]
xs) [Integer]
ys -> Integer
x0 Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: Integer
x1 Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: Integer
x2 Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: (Integer -> Integer -> Bool) -> [Integer] -> [Integer] -> [Integer]
forall a. (a -> a -> Bool) -> [a] -> [a] -> [a]
ListHT.mergeBy Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
(<=) [Integer]
xs [Integer]
ys)
      ([Char] -> [Integer]
forall a. HasCallStack => [Char] -> a
error [Char]
"numbers5SmoothFoldr: infinite list should not have an end") ([[Integer]] -> [Integer]) -> [[Integer]] -> [Integer]
forall a b. (a -> b) -> a -> b
$
   ([Integer] -> [Integer]) -> [Integer] -> [[Integer]]
forall a. (a -> a) -> a -> [a]
iterate ((Integer -> Integer) -> [Integer] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Integer
5Integer -> Integer -> Integer
forall a. C a => a -> a -> a
*)) ([Integer] -> [[Integer]]) -> [Integer] -> [[Integer]]
forall a b. (a -> b) -> a -> b
$
   [Integer]
numbers3SmoothFoldr

numbers5SmoothSet :: [Integer]
numbers5SmoothSet :: [Integer]
numbers5SmoothSet =
   (Set Integer -> Maybe (Integer, Set Integer))
-> Set Integer -> [Integer]
forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr
      (((Integer, Set Integer) -> (Integer, Set Integer))
-> Maybe (Integer, Set Integer) -> Maybe (Integer, Set Integer)
forall a b. (a -> b) -> Maybe a -> Maybe b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\(Integer
m,Set Integer
rest) -> (Integer
m, Set Integer -> Set Integer -> Set Integer
forall a. Ord a => Set a -> Set a -> Set a
Set.union Set Integer
rest (Set Integer -> Set Integer) -> Set Integer -> Set Integer
forall a b. (a -> b) -> a -> b
$ [Integer] -> Set Integer
forall a. Eq a => [a] -> Set a
Set.fromAscList [Integer
2Integer -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
m,Integer
3Integer -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
m,Integer
5Integer -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
m])) (Maybe (Integer, Set Integer) -> Maybe (Integer, Set Integer))
-> (Set Integer -> Maybe (Integer, Set Integer))
-> Set Integer
-> Maybe (Integer, Set Integer)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       Set Integer -> Maybe (Integer, Set Integer)
forall a. Set a -> Maybe (a, Set a)
Set.minView) (Set Integer -> [Integer]) -> Set Integer -> [Integer]
forall a b. (a -> b) -> a -> b
$
   Integer -> Set Integer
forall a. a -> Set a
Set.singleton Integer
1

ceilingPowerOfTwo :: (Ring.C a, Bits a) => a -> a
ceilingPowerOfTwo :: forall a. (C a, Bits a) => a -> a
ceilingPowerOfTwo a
0 = a
1
ceilingPowerOfTwo a
n =
   (a
1a -> a -> a
forall a. C a => a -> a -> a
+) (a -> a) -> a -> a
forall a b. (a -> b) -> a -> b
$ (a, a) -> a
forall a b. (a, b) -> a
fst ((a, a) -> a) -> (a, a) -> a
forall a b. (a -> b) -> a -> b
$ [(a, a)] -> (a, a)
forall a. HasCallStack => [a] -> a
head ([(a, a)] -> (a, a)) -> [(a, a)] -> (a, a)
forall a b. (a -> b) -> a -> b
$
   ((a, a) -> Bool) -> [(a, a)] -> [(a, a)]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile ((a -> a -> Bool) -> (a, a) -> Bool
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry a -> a -> Bool
forall a. Eq a => a -> a -> Bool
(/=)) ([(a, a)] -> [(a, a)]) -> [(a, a)] -> [(a, a)]
forall a b. (a -> b) -> a -> b
$
   (a -> a -> (a, a)) -> [a] -> [(a, a)]
forall a b. (a -> a -> b) -> [a] -> [b]
ListHT.mapAdjacent (,) ([a] -> [(a, a)]) -> [a] -> [(a, a)]
forall a b. (a -> b) -> a -> b
$
   (a -> Int -> a) -> a -> [Int] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\a
m Int
d -> a -> Int -> a
forall a. Bits a => a -> Int -> a
shiftR a
m Int
d a -> a -> a
forall a. Bits a => a -> a -> a
.|. a
m) (a
na -> a -> a
forall a. C a => a -> a -> a
-a
1) ([Int] -> [a]) -> [Int] -> [a]
forall a b. (a -> b) -> a -> b
$
   (Int -> Int) -> Int -> [Int]
forall a. (a -> a) -> a -> [a]
iterate (Int
2Int -> Int -> Int
forall a. C a => a -> a -> a
*) Int
1

{- |
It's not awfully efficient, but ok for our uses.
-}
ceilingPower :: (Integral.C a, Ord a) => a -> a -> a
ceilingPower :: forall a. (C a, Ord a) => a -> a -> a
ceilingPower a
base a
n = a
base a -> Integer -> a
forall a. C a => a -> Integer -> a
^ Int -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral (a -> a -> Int
forall a. (C a, Ord a) => a -> a -> Int
ceilingLog a
base a
n)

ceilingLog :: (Integral.C a, Ord a) => a -> a -> Int
ceilingLog :: forall a. (C a, Ord a) => a -> a -> Int
ceilingLog a
base =
   [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([a] -> Int) -> (a -> [a]) -> a -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>a
0) ([a] -> [a]) -> (a -> [a]) -> a -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate ((a -> a -> a) -> a -> a -> a
forall a b c. (a -> b -> c) -> b -> a -> c
flip a -> a -> a
forall a. C a => a -> a -> a
div a
base) (a -> [a]) -> (a -> a) -> a -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> a -> a
forall a. C a => a -> a -> a
subtract a
1

divideByMaximumPower ::
   (Integral.C a, ZeroTestable.C a) => a -> a -> a
divideByMaximumPower :: forall a. (C a, C a) => a -> a -> a
divideByMaximumPower a
b a
n =
   [a] -> a
forall a. HasCallStack => [a] -> a
last ([a] -> a) -> [a] -> a
forall a b. (a -> b) -> a -> b
$
   a
n a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> Maybe (a, a)) -> a -> [a]
forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr (\a
m -> case a -> a -> (a, a)
forall a. C a => a -> a -> (a, a)
divMod a
m a
b of (a
q,a
r) -> Bool -> (a, a) -> Maybe (a, a)
forall a. Bool -> a -> Maybe a
toMaybe (a -> Bool
forall a. C a => a -> Bool
isZero a
r) (a
q,a
q)) a
n

divideByMaximumPowerRecursive ::
   (Integral.C a, Eq a, ZeroTestable.C a) => a -> a -> a
divideByMaximumPowerRecursive :: forall a. (C a, Eq a, C a) => a -> a -> a
divideByMaximumPowerRecursive a
b =
   let recourse :: a -> a
recourse a
n =
          case a -> a -> (a, a)
forall a. C a => a -> a -> (a, a)
divMod a
n a
b of
             (a
q,a
0) -> a -> a
recourse a
q
             (a, a)
_ -> a
n
   in  a -> a
recourse

getMaximumExponent ::
   (Integral.C a, ZeroTestable.C a) =>
   a -> a -> (Int,a)
getMaximumExponent :: forall a. (C a, C a) => a -> a -> (Int, a)
getMaximumExponent a
b a
n =
   [(Int, a)] -> (Int, a)
forall a. HasCallStack => [a] -> a
last ([(Int, a)] -> (Int, a)) -> [(Int, a)] -> (Int, a)
forall a b. (a -> b) -> a -> b
$ (Int
0,a
n) (Int, a) -> [(Int, a)] -> [(Int, a)]
forall a. a -> [a] -> [a]
:
   ((Int, a) -> Maybe ((Int, a), (Int, a))) -> (Int, a) -> [(Int, a)]
forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr
      (\(Int
e,a
m) ->
         let (a
q,a
r) = a -> a -> (a, a)
forall a. C a => a -> a -> (a, a)
divMod a
m a
b
             eq :: (Int, a)
eq = (Int
eInt -> Int -> Int
forall a. C a => a -> a -> a
+Int
1,a
q)
         in  Bool -> ((Int, a), (Int, a)) -> Maybe ((Int, a), (Int, a))
forall a. Bool -> a -> Maybe a
toMaybe (a -> Bool
forall a. C a => a -> Bool
isZero a
r) ((Int, a)
eq,(Int, a)
eq))
      (Int
0,a
n)

is3Smooth :: Integer -> Bool
is3Smooth :: Integer -> Bool
is3Smooth =
   (Integer
1Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==) (Integer -> Bool) -> (Integer -> Integer) -> Integer -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   Integer -> Integer -> Integer
forall a. (C a, C a) => a -> a -> a
divideByMaximumPower Integer
3 (Integer -> Integer) -> (Integer -> Integer) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   Integer -> Integer -> Integer
forall a. (C a, C a) => a -> a -> a
divideByMaximumPower Integer
2

is5Smooth :: Integer -> Bool
is5Smooth :: Integer -> Bool
is5Smooth =
   (Integer
1Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==) (Integer -> Bool) -> (Integer -> Integer) -> Integer -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   Integer -> Integer -> Integer
forall a. (C a, C a) => a -> a -> a
divideByMaximumPower Integer
5 (Integer -> Integer) -> (Integer -> Integer) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   Integer -> Integer -> Integer
forall a. (C a, C a) => a -> a -> a
divideByMaximumPower Integer
3 (Integer -> Integer) -> (Integer -> Integer) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   Integer -> Integer -> Integer
forall a. (C a, C a) => a -> a -> a
divideByMaximumPower Integer
2


ceiling3Smooth :: Integer -> Integer
ceiling3Smooth :: Integer -> Integer
ceiling3Smooth = Integer -> Integer
ceiling3SmoothTrace

ceiling5Smooth :: Integer -> Integer
ceiling5Smooth :: Integer -> Integer
ceiling5Smooth = Integer -> Integer
ceiling5SmoothTrace

{- |
Compute the smallest composite of 2 and 3 that is as least as large as the input.
This can be interpreted as solving an integer linear programming problem with
min (\(a,b) -> a * log 2 + b * log 3)
over the domain {(a,b) : a>=0, b>=0, a * log 2 + b * log 3 >= log n}
-}
{-
This implementation looks stupid,
but it is drastically faster for large numbers than ceiling3SmoothNaive.
The reason is that the smooth numbers are logarithmically equally distributed.
That is, from @n@ to the next smooth number
it may be only 1% deviation from @n@,
but for huge @n@ the absolute difference @0.01*n@ is still huge.

@ceiling3Smooth (10^400+1)@ can be computed in about 0.1 seconds.
(Surprisingly, @ceiling3Smooth (10^500+1)@ needs almost 30 seconds.)
-}
ceiling3SmoothScan :: Integer -> Integer
ceiling3SmoothScan :: Integer -> Integer
ceiling3SmoothScan Integer
n =
   [Integer] -> Integer
forall a. HasCallStack => [a] -> a
head ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$ (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<Integer
n) [Integer]
numbers3Smooth

ceiling5SmoothScan :: Integer -> Integer
ceiling5SmoothScan :: Integer -> Integer
ceiling5SmoothScan Integer
n =
   [Integer] -> Integer
forall a. HasCallStack => [a] -> a
head ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$ (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<Integer
n) [Integer]
numbers5Smooth

ceiling3SmoothNaive :: Integer -> Integer
ceiling3SmoothNaive :: Integer -> Integer
ceiling3SmoothNaive =
   [Integer] -> Integer
forall a. HasCallStack => [a] -> a
head ([Integer] -> Integer)
-> (Integer -> [Integer]) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (Bool -> Bool
not (Bool -> Bool) -> (Integer -> Bool) -> Integer -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Bool
is3Smooth) ([Integer] -> [Integer])
-> (Integer -> [Integer]) -> Integer -> [Integer]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer
1Integer -> Integer -> Integer
forall a. C a => a -> a -> a
+)

ceiling5SmoothNaive :: Integer -> Integer
ceiling5SmoothNaive :: Integer -> Integer
ceiling5SmoothNaive =
   [Integer] -> Integer
forall a. HasCallStack => [a] -> a
head ([Integer] -> Integer)
-> (Integer -> [Integer]) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (Bool -> Bool
not (Bool -> Bool) -> (Integer -> Bool) -> Integer -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Bool
is5Smooth) ([Integer] -> [Integer])
-> (Integer -> [Integer]) -> Integer -> [Integer]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer
1Integer -> Integer -> Integer
forall a. C a => a -> a -> a
+)


{-
Problem: We cannot just start with the ceilingPowerOfTwo
and then multiply with 3/4 until we fall below n,
since the 3/4 decreases too fast.
27/32 is closer to one,
and higher powers of 3 and 2 in the ratio make the ratio even closer to one.

This implementation is different:
It always moves and tests above @n@.
For every power of 3 it computes the least power of 2,
such that their product is above @n@.
-}
ceiling3SmoothTrace :: Integer -> Integer
ceiling3SmoothTrace :: Integer -> Integer
ceiling3SmoothTrace Integer
n =
   [Integer] -> Integer
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$ Integer -> Integer -> Integer -> Integer -> [Integer]
ceilingSmoothsTrace Integer
2 Integer
3 Integer
n (Integer -> [Integer]) -> Integer -> [Integer]
forall a b. (a -> b) -> a -> b
$ Integer -> Integer
forall a. (C a, Bits a) => a -> a
ceilingPowerOfTwo Integer
n

{-
We must be careful not to skip combinations that are optimal.

E.g.:
> _ceiling5SmoothTraceWrong (10^70+1)
10002658207445093206727527411583349735126415100956607165326185795158016
> ceiling5Smooth (10^70+1)
10001329015408448808646079907338649600000000000000000000000000000000000
-}
_ceiling5SmoothTraceWrong :: Integer -> Integer
_ceiling5SmoothTraceWrong :: Integer -> Integer
_ceiling5SmoothTraceWrong Integer
n =
   [Integer] -> Integer
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer) -> [Integer] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map ([Integer] -> Integer
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum ([Integer] -> Integer)
-> (Integer -> [Integer]) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Integer -> Integer -> Integer -> [Integer]
ceilingSmoothsTrace Integer
3 Integer
5 Integer
n) ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$
   Integer -> Integer -> Integer -> Integer -> [Integer]
ceilingSmoothsTrace Integer
2 Integer
3 Integer
n (Integer -> [Integer]) -> Integer -> [Integer]
forall a b. (a -> b) -> a -> b
$ Integer -> Integer
forall a. (C a, Bits a) => a -> a
ceilingPowerOfTwo Integer
n

{-
For every reasonable pair of powers of 3 and 5
it computes the least power of 2,
such that their product is above @n@.
-}
ceiling5SmoothTrace :: Integer -> Integer
ceiling5SmoothTrace :: Integer -> Integer
ceiling5SmoothTrace Integer
n =
   [Integer] -> Integer
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer) -> [Integer] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map ([Integer] -> Integer
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum ([Integer] -> Integer)
-> (Integer -> [Integer]) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Integer -> Integer -> Integer -> [Integer]
ceilingSmoothsTrace Integer
2 Integer
5 Integer
n) ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$
   Integer -> Integer -> Integer -> Integer -> [Integer]
ceilingSmoothsTrace Integer
2 Integer
3 Integer
n (Integer -> [Integer]) -> Integer -> [Integer]
forall a b. (a -> b) -> a -> b
$ Integer -> Integer
forall a. (C a, Bits a) => a -> a
ceilingPowerOfTwo Integer
n

{- |
@ceilingSmoothsTrace a b n m@
replaces successively @a@ factors in @m@ by @b@ factors
while keeping the product above @n@.
-}
ceilingSmoothsTrace :: Integer -> Integer -> Integer -> Integer -> [Integer]
ceilingSmoothsTrace :: Integer -> Integer -> Integer -> Integer -> [Integer]
ceilingSmoothsTrace Integer
a Integer
b Integer
n =
   let divMany :: Integer -> Integer
divMany Integer
k =
          case Integer -> Integer -> (Integer, Integer)
forall a. C a => a -> a -> (a, a)
divMod Integer
k Integer
a of
             (Integer
q,Integer
r) -> if Integer
rInteger -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==Integer
0 Bool -> Bool -> Bool
&& Integer
qInteger -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>=Integer
n then Integer -> Integer
divMany Integer
q else Integer
k
       go :: Integer -> [Integer]
go Integer
m  =  Integer
m Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: if Integer -> Integer -> Integer
forall a. C a => a -> a -> a
mod Integer
m Integer
a Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 then Integer -> [Integer]
go (Integer -> [Integer]) -> Integer -> [Integer]
forall a b. (a -> b) -> a -> b
$ Integer -> Integer
divMany (Integer -> Integer) -> Integer -> Integer
forall a b. (a -> b) -> a -> b
$ Integer
mInteger -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
b else []
   in  Integer -> [Integer]
go


{- |
Compute all primes that occur in the course of dividing
a Fourier transform into sub-transforms.
-}
partialPrimes :: Integer -> [Integer]
partialPrimes :: Integer -> [Integer]
partialPrimes =
   let primeFactorSet :: Integer -> Set Integer
primeFactorSet = [Integer] -> Set Integer
forall a. Eq a => [a] -> Set a
Set.fromAscList ([Integer] -> Set Integer)
-> (Integer -> [Integer]) -> Integer -> Set Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> [Integer]
forall a. (C a, Bits a, C a, Ord a) => a -> [a]
uniquePrimeFactors
   in  (Set Integer -> Maybe (Integer, Set Integer))
-> Set Integer -> [Integer]
forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr
         (((Integer, Set Integer) -> (Integer, Set Integer))
-> Maybe (Integer, Set Integer) -> Maybe (Integer, Set Integer)
forall a b. (a -> b) -> Maybe a -> Maybe b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap
             (\(Integer
p,Set Integer
set) ->
                (Integer
p, Set Integer -> Set Integer -> Set Integer
forall a. Ord a => Set a -> Set a -> Set a
Set.union (Integer -> Set Integer
primeFactorSet (Integer
pInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
1)) Set Integer
set)) (Maybe (Integer, Set Integer) -> Maybe (Integer, Set Integer))
-> (Set Integer -> Maybe (Integer, Set Integer))
-> Set Integer
-> Maybe (Integer, Set Integer)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
          Set Integer -> Maybe (Integer, Set Integer)
forall a. Set a -> Maybe (a, Set a)
Set.maxView)
       (Set Integer -> [Integer])
-> (Integer -> Set Integer) -> Integer -> [Integer]
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       Integer -> Set Integer
primeFactorSet

-- cf. htam:NumberTheory
uniquePrimeFactors ::
   (Integral.C a, Bits a, ZeroTestable.C a, Ord a) =>
   a -> [a]
--   map snd . primeFactors
uniquePrimeFactors :: forall a. (C a, Bits a, C a, Ord a) => a -> [a]
uniquePrimeFactors a
n =
   let oddFactors :: a -> [a]
oddFactors =
          (a -> (a -> [a]) -> a -> [a]) -> (a -> [a]) -> [a] -> a -> [a]
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr
             (\a
p a -> [a]
go a
m ->
                let (a
q,a
r) = a -> a -> (a, a)
forall a. C a => a -> a -> (a, a)
divMod a
m a
p
                in  if a
ra -> a -> Bool
forall a. Eq a => a -> a -> Bool
==a
0
                      then a
p a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a -> [a]
go (a -> a -> a
forall a. (C a, C a) => a -> a -> a
divideByMaximumPower a
p a
q)
                      else
                        if a
q a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
p
                          then a -> [a]
go a
m
                          else if a
ma -> a -> Bool
forall a. Eq a => a -> a -> Bool
==a
1 then [] else a
m a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [])
             ([Char] -> a -> [a]
forall a. HasCallStack => [Char] -> a
error [Char]
"uniquePrimeFactors: end of infinite list")
             ((a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a
2a -> a -> a
forall a. C a => a -> a -> a
+) a
3)
   in  case a -> (a, a)
forall a. (Bits a, C a) => a -> (a, a)
powerOfTwoFactors a
n of
          (a
1,a
m) -> a -> [a]
oddFactors a
m
          (a
_,a
m) -> a
2 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a -> [a]
oddFactors a
m

{- |
Prime factors and their exponents in ascending order.
-}
primeFactors ::
   (PrimitiveRoot a, Ord a) => a -> [(Int, a)]
primeFactors :: forall a. (PrimitiveRoot a, Ord a) => a -> [(Int, a)]
primeFactors a
n =
   let oddFactors :: a -> [(Int, a)]
oddFactors =
          (a -> (a -> [(Int, a)]) -> a -> [(Int, a)])
-> (a -> [(Int, a)]) -> [a] -> a -> [(Int, a)]
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr
             (\a
p a -> [(Int, a)]
go a
m ->
                let (a
q0,a
r) = a -> a -> (a, a)
forall a. C a => a -> a -> (a, a)
divMod a
m a
p
                in  if a
ra -> a -> Bool
forall a. Eq a => a -> a -> Bool
==a
0
                      then
                        case a -> a -> (Int, a)
forall a. (C a, C a) => a -> a -> (Int, a)
getMaximumExponent a
p a
q0 of
                          (Int
e,a
q1) -> (Int
eInt -> Int -> Int
forall a. C a => a -> a -> a
+Int
1,a
p) (Int, a) -> [(Int, a)] -> [(Int, a)]
forall a. a -> [a] -> [a]
: a -> [(Int, a)]
go a
q1
                      else
                        if a
q0 a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
p
                          then a -> [(Int, a)]
go a
m
                          else if a
ma -> a -> Bool
forall a. Eq a => a -> a -> Bool
==a
1 then [] else (Int
1,a
m) (Int, a) -> [(Int, a)] -> [(Int, a)]
forall a. a -> [a] -> [a]
: [])
             ([(Int, a)] -> a -> [(Int, a)]
forall a b. a -> b -> a
const [])
             ((a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
filter (Bool -> Bool
not (Bool -> Bool) -> (a -> Bool) -> a -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> Bool
forall a. C a => a -> Bool
Units.isUnit) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$
              a -> [a]
forall a. PrimitiveRoot a => a -> [a]
primitiveRootCandidates a
n)
   in  case a -> a -> (Int, a)
forall a. (C a, C a) => a -> a -> (Int, a)
getMaximumExponent a
2 a
n of
          (Int
0,a
m) -> a -> [(Int, a)]
oddFactors a
m
          (Int
e,a
m) -> (Int
e,a
2) (Int, a) -> [(Int, a)] -> [(Int, a)]
forall a. a -> [a] -> [a]
: a -> [(Int, a)]
oddFactors a
m

{-
cf. htam:NumberTheory

Shall this be moved to NumericPrelude?

It should be replaced by a more sophisticated prime test.
-}
isPrime :: Integer -> Bool
isPrime :: Integer -> Bool
isPrime Integer
n =
   case Integer -> [(Int, Integer)]
forall a. (PrimitiveRoot a, Ord a) => a -> [(Int, a)]
primeFactors Integer
n of
      [] -> Bool
False
      (Int
e,Integer
m):[(Int, Integer)]
_ -> Int
eInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
1 Bool -> Bool -> Bool
&& Integer
mInteger -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==Integer
n

{- |
Find lengths of signals that require many interim Rader transforms
and end with the given length.

> raderWorstCases 2  =  <http://oeis.org/A061092>
> raderWorstCases 5  =  tail <http://oeis.org/A059411>

Smallest raderWorstCase numbers are 2,5,13,17,19,31,37,41,43,61,...
This matches the definition of <http://oeis.org/A061303>.
-}
raderWorstCases :: Integer -> [Integer]
raderWorstCases :: Integer -> [Integer]
raderWorstCases =
   (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate
      (\Integer
n ->
         [Integer] -> Integer
forall a. HasCallStack => [a] -> a
head ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$ (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (Bool -> Bool
not (Bool -> Bool) -> (Integer -> Bool) -> Integer -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Bool
isPrime) ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$
         [Integer] -> [Integer]
forall a. HasCallStack => [a] -> [a]
tail ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer
nInteger -> Integer -> Integer
forall a. C a => a -> a -> a
+) Integer
1)

{- |
This is usually faster than 'fastFourierRing'
since it does not need to factor large numbers.
However, the generated modulus is usually much greater.
-}
{-
I see the following opportunities for optimization:

1. Speedup 'fastFourierRing' by
   faster primality test (Miller-Rabin) and
   faster prime factorization (Pollard-Rho-method).
   These are also needed for
   maximumOrderOfPrimitiveRootsOfUnityInteger
   that is used by Fourier.Element.primitiveRoot
   in order to compute a root with maximum order.

2. Reduce the moduli produced by '_fastFourierRingAlt'
   by merging some orders that are passed to
   ringWithPrimitiveRootsOfUnityAndUnits,
   such that an LCM of a group of orders can still be handled.
   This is a kind of knapsack problem.
   Maybe we could collect the factors in a way
   such that (lcm orderGroup + 1) is prime.

3. Avoid to compute factorizations of numbers
   where we already know the factors
   or where we do not need the factors at all.
   Use the factors returned by partialPrimes
   in order to compute a prime factorization
   of lcmMulti (map pred (partialPrimes n)).
   Call this (product ps).
   Now search for rings of moduli (1 + k * product ps),
   where there are (small) primitive roots of order (product ps).
   We only need to check whether there are small numbers
   such as 2, 3, 5, 6, 7 that have a (product ps)-th power that is 1,
   using fast exponentiation.
   If there is a power being 1 then check for primitivity
   by computing (k * product ps / p)-th powers
   for all prime factors p of (k * product ps).
   If there is no primitive root <= 7,
   there may still be a primitive root of wanted order,
   but it is then cheaper to seek for larger moduli.

   If we finally have a nice modulus
   it is still stupid to factorize (modulus-1)
   and search for a primitive root
   in each invocation of Fourier.Element.primitiveRoot.
   We could define a special datatype analogously to ResidueClass,
   that stores the primitive root and its order
   additional to the ResidueClass modulus.
-}
_fastFourierRingAlt :: Int -> Integer
_fastFourierRingAlt :: Int -> Integer
_fastFourierRingAlt Int
n =
   case Int
n of
      Int
0 -> Integer
2
      Int
1 -> Integer
2
      Int
_ ->
         let ni :: Integer
ni = Int -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral Int
n
             ps :: [Integer]
ps = (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
filter (Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>Integer
1) ((Integer -> Integer) -> [Integer] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Integer -> Integer -> Integer
forall a. C a => a -> a -> a
subtract Integer
1) (Integer -> [Integer]
partialPrimes Integer
ni))
         in  [Order] -> [Integer] -> Integer
ringWithPrimitiveRootsOfUnityAndUnits ((Integer -> Order) -> [Integer] -> [Order]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> Order
Order ([Integer] -> [Order]) -> [Integer] -> [Order]
forall a b. (a -> b) -> a -> b
$ Integer
ni Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: [Integer]
ps) [Integer]
ps

{- |
Determine an integer residue ring
in which a Fast Fourier transform of size n can be performed.
It must contain certain primitive roots.
If we choose a non-primitive root,
then some off-diagonal values in F^-1·F are non-zero.
-}
{-
When we need roots of orders o1,...,ok and according inverse elements
we can also ask for a ring, where there is a root of order lcm(o1,...,ok).
The answer to both questions is the same set of rings.
This can be proven using the statement,
that the order of any primitive root
divides the carmichael value of the modulus.

Since ringWithPrimitiveRootsOfUnityAndUnits
multiplies the moduli of rings for o1,...,ok,
it will produce large moduli.
-}
fastFourierRing :: Int -> Integer
fastFourierRing :: Int -> Integer
fastFourierRing Int
n =
   case Int
n of
      Int
0 -> Integer
2
      Int
1 -> Integer
2
      Int
_ ->
         let ni :: Integer
ni = Int -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral Int
n
         in  {-
             We cannot use ringsWithPrimitiveRootOfUnityAndUnit
             since for 359 we already get an Int overflow.
             For 719, 1439, 2879 we also get a very large value.
             -}
             [Integer] -> Integer
forall a. HasCallStack => [a] -> a
head ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$ (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
filter Integer -> Bool
isPrime ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$
             (\Integer
order -> (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer
order Integer -> Integer -> Integer
forall a. C a => a -> a -> a
+) Integer
1) (Integer -> [Integer]) -> Integer -> [Integer]
forall a b. (a -> b) -> a -> b
$
             [Integer] -> Integer
forall a. C a => [a] -> a
lcmMulti ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$
             Integer
ni Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: (Integer -> Integer) -> [Integer] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Integer -> Integer -> Integer
forall a. C a => a -> a -> a
subtract Integer
1) (Integer -> [Integer]
partialPrimes Integer
ni)