{-# 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 forall a. Bits a => a -> a -> a
.&. (-a
n)
   in  (a
powerOfTwo, 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  forall a b. (a -> b) -> [a] -> [b]
map (\(Integer
a,Integer
b) -> (Integer
bforall a. C a => a -> a -> a
-Integer
a,Integer
bforall a. C a => a -> a -> a
+Integer
a)) forall a b. (a -> b) -> a -> b
$
       forall a b c. Ord a => [(a, b)] -> [(a, c)] -> [(b, c)]
mergeAndFilter
          (forall a b. [a] -> [b] -> [(a, b)]
zip (forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl forall a. C a => a -> a -> a
(+) Integer
n [Integer
1,Integer
3..]) [Integer
0 .. forall a. C a => a -> a -> a
div (Integer
nforall a. C a => a -> a -> a
-Integer
1) Integer
2])
          (forall a b. [a] -> [b] -> [(a, b)]
zip (forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl forall a. C a => a -> a -> a
(+) (Integer
rootforall a. C a => a -> a -> a
*Integer
root) forall a b. (a -> b) -> a -> b
$ forall a. (a -> a) -> a -> [a]
iterate (Integer
2forall a. C a => a -> a -> a
+) (Integer
2forall a. C a => a -> a -> a
*Integer
rootforall 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 forall a. Ord a => a -> a -> Ordering
compare a
a0 a
a1 of
      Ordering
LT -> forall a b c. Ord a => [(a, b)] -> [(a, c)] -> [(b, c)]
mergeAndFilter [(a, b)]
a0s ((a
a1,c
c)forall a. a -> [a] -> [a]
:[(a, c)]
a1s)
      Ordering
GT -> forall a b c. Ord a => [(a, b)] -> [(a, c)] -> [(b, c)]
mergeAndFilter ((a
a0,b
b)forall a. a -> [a] -> [a]
:[(a, b)]
a0s) [(a, c)]
a1s
      Ordering
EQ -> (b
b,c
c) forall a. a -> [a] -> [a]
: 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 forall a. Set a -> Maybe (a, Set a)
Set.minView Set Integer
candidates of
             Maybe (Integer, Set Integer)
Nothing -> forall a. HasCallStack => [Char] -> a
error forall a b. (a -> b) -> a -> b
$ forall a. Show a => a -> [Char]
show Integer
p forall a. [a] -> [a] -> [a]
++ [Char]
" is not a prime"
             Just (Integer
x,Set Integer
rest) ->
                case forall a. Ord a => T a -> Set a
orbitSet forall a b. (a -> b) -> a -> b
$ 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 forall a. Eq a => a -> a -> Bool
== forall a. Ord a => [a] -> Set a
Set.fromList [Integer
1..Integer
pforall a. C a => a -> a -> a
-Integer
1]
                        then Integer
x
                        else Set Integer -> Integer
search (forall a. Ord a => Set a -> Set a -> Set a
Set.difference Set Integer
rest Set Integer
new)
   in  Set Integer -> Integer
search forall a b. (a -> b) -> a -> b
$ forall a. Ord a => [a] -> Set a
Set.fromList [Integer
1..Integer
pforall a. C a => a -> a -> a
-Integer
1]

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


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

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

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

countOrder :: [a] -> Order
countOrder :: forall a. [a] -> Order
countOrder = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (\Order
o a
_ -> 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) =
   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
moduforall 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 =
   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 = forall a b. (a -> b) -> [a] -> [b]
map (forall a. C a => a -> a -> a
div Integer
order) forall a b. (a -> b) -> a -> b
$ forall a. (C a, Bits a, C a, Ord a) => a -> [a]
uniquePrimeFactors Integer
order
   in  forall a. (a -> Bool) -> [a] -> [a]
filter
          (\a
n ->
             let pow :: Integer -> a
pow Integer
y = forall a. T a -> a
RC.representative forall a b. (a -> b) -> a -> b
$ (a
n forall a. C a => a -> a -> T a
/: a
modu) forall a. C a => a -> Integer -> a
^ Integer
y
             in  forall a. C a => a -> a -> Bool
PID.coprime a
n a
modu
                 Bool -> Bool -> Bool
&&
                 Integer -> a
pow Integer
order forall a. Eq a => a -> a -> Bool
== forall a. C a => a
one
                 Bool -> Bool -> Bool
&&
                 forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Integer
y -> Integer -> a
pow Integer
y forall a. Eq a => a -> a -> Bool
/= forall a. C a => a
one) [Integer]
greatDivisors) forall a b. (a -> b) -> a -> b
$
       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) =
   forall a. (a -> Bool) -> [a] -> [a]
filter
      (\a
n ->
         let ([a]
prefix,a
end:[a]
_) =
                forall i a. Integral i => i -> [a] -> ([a], [a])
genericSplitAt (Integer
expoforall a. C a => a -> a -> a
-Integer
1) forall a b. (a -> b) -> a -> b
$ forall y. T y -> [y]
SigS.toList forall a b. (a -> b) -> a -> b
$ forall a. C a => a -> a -> T a
orbit a
modu a
n
         in  forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (a
1forall a. Eq a => a -> a -> Bool
/=) [a]
prefix Bool -> Bool -> Bool
&& a
endforall a. Eq a => a -> a -> Bool
==a
1) forall a b. (a -> b) -> a -> b
$
   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 =
   forall x acc. (x -> acc -> acc) -> acc -> T x -> acc
SigS.foldR
      (\a
new Set a -> Set a
cont Set a
seen ->
         if 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 (forall a. Ord a => a -> Set a -> Set a
Set.insert a
new Set a
seen))
      forall a. a -> a
id T a
list 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 = forall a. (a -> a) -> a -> T a
SigS.iterate (\a
y -> forall a. C a => a -> a -> a
mod (a
xforall 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 =
          forall a b c. (a -> b -> c) -> b -> a -> c
flip forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a. Set a -> Maybe (a, Set a)
Set.minView Set a
candidates) forall a b. (a -> b) -> a -> b
$ \(a
x,Set a
rest) ->
          forall b c a. (b -> c) -> (a, b) -> (a, c)
mapSnd (forall a. Ord a => Set a -> Set a -> Set a
Set.difference Set a
rest forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Ord a => [a] -> Set a
Set.fromList) forall a b. (a -> b) -> a -> b
$
          forall a. (PrimitiveRoot a, Ord a) => a -> Order -> a -> ([a], [a])
primitiveRootsOfOrbit a
modu Order
expo a
x
   in  forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat forall a b. (a -> b) -> a -> b
$ forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr Set a -> Maybe ([a], Set a)
search forall a b. (a -> b) -> a -> b
$ forall a. Ord a => [a] -> Set a
Set.fromList forall a b. (a -> b) -> a -> b
$
       -- needed for modules with many divisors
       forall a. (a -> Bool) -> [a] -> [a]
filter (forall a. C a => a -> a -> Bool
PID.coprime a
modu) forall a b. (a -> b) -> a -> b
$
       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 =
          forall a b c. (a -> b -> c) -> b -> a -> c
flip forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a. Set a -> Maybe (a, Set a)
Set.minView Set a
candidates) forall a b. (a -> b) -> a -> b
$ \(a
x,Set a
rest) ->
          forall a c b d. (a -> c, b -> d) -> (a, b) -> (c, d)
mapPair ((,) a
x,
                   forall a. Ord a => Set a -> Set a -> Set a
Set.difference Set a
rest forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Ord a => [a] -> Set a
Set.fromList) forall a b. (a -> b) -> a -> b
$
          forall a. (PrimitiveRoot a, Ord a) => a -> Order -> a -> ([a], [a])
primitiveRootsOfOrbit a
modu Order
expo a
x
   in  forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr Set a -> Maybe ((a, [a]), Set a)
search forall a b. (a -> b) -> a -> b
$ forall a. Ord a => [a] -> Set a
Set.fromList forall a b. (a -> b) -> a -> b
$
       forall a. (a -> Bool) -> [a] -> [a]
filter (forall a. C a => a -> a -> Bool
PID.coprime a
modu) forall a b. (a -> b) -> a -> b
$
       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
1forall a. a -> [a] -> [a]
:) forall a b. (a -> b) -> a -> b
$ forall a. (a -> Bool) -> [a] -> [a]
takeWhile (a
1forall a. Eq a => a -> a -> Bool
/=) forall a b. (a -> b) -> a -> b
$ forall a. (a -> a) -> a -> [a]
iterate (\a
y -> forall a. C a => a -> a -> a
mod (a
xforall a. C a => a -> a -> a
*a
y) a
modu) a
x
       (Order Integer
orbitSize) = forall a. [a] -> Order
countOrder [a]
orb
   in  (if Integer
expoforall 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 forall a. C a => a -> a -> (a, a)
divMod Integer
orbitSize Integer
expo of
               (Integer
s,Integer
0) ->
                  forall a b. (a -> b) -> [a] -> [b]
map forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall a. (a -> Bool) -> [a] -> [a]
filter (forall a. C a => a -> a -> Bool
PID.coprime Integer
expo forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a, b) -> a
fst) forall a b. (a -> b) -> a -> b
$
                  forall a b. [a] -> [b] -> [(a, b)]
zip
                     [Integer
0 .. Integer
expoforall a. C a => a -> a -> a
-Integer
1]
                     -- (ListHT.sieve s $ orb)
                     (forall a b. (a -> b) -> [a] -> [b]
map forall a. [a] -> a
head forall a b. (a -> b) -> a -> b
$ forall a. (a -> a) -> a -> [a]
iterate (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 =
   forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Order -> Order -> Bool
dividesOrder Order
expo forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a, b) -> b
snd) forall a b. (a -> b) -> a -> b
$
   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 =
   forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl forall a. Ord a => a -> a -> a
max (Integer -> Order
Order Integer
1) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> b) -> [a] -> [b]
map forall a b. (a, b) -> b
snd forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 =
          forall a b c. (a -> b -> c) -> b -> a -> c
flip forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a. Set a -> Maybe (a, Set a)
Set.minView Set a
candidates) forall a b. (a -> b) -> a -> b
$ \(a
x,Set a
rest) ->
          forall a c b d. (a -> c, b -> d) -> (a, b) -> (c, d)
mapPair ((,) a
x,
                   forall a. Ord a => Set a -> Set a -> Set a
Set.difference Set a
rest forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Ord a => [a] -> Set a
Set.fromList) forall a b. (a -> b) -> a -> b
$
          forall a. (PrimitiveRoot a, Ord a) => a -> a -> (Order, [a])
orderOfOrbit a
modu a
x
   in  forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr Set a -> Maybe ((a, Order), Set a)
search forall a b. (a -> b) -> a -> b
$ forall a. Ord a => [a] -> Set a
Set.fromList forall a b. (a -> b) -> a -> b
$
       forall a. (a -> Bool) -> [a] -> [a]
filter (forall a. C a => a -> a -> Bool
PID.coprime a
modu) forall a b. (a -> b) -> a -> b
$
       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 = forall a. (a -> Bool) -> [a] -> [a]
takeWhile (forall a. C a => a
oneforall a. Eq a => a -> a -> Bool
/=) forall a b. (a -> b) -> a -> b
$ forall y. T y -> [y]
SigS.toList forall a b. (a -> b) -> a -> b
$ forall a. C a => a -> a -> T a
orbit a
modu a
x
   in  (forall a. Enum a => a -> a
succ forall a b. (a -> b) -> a -> b
$ 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 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 forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   forall a. C a => [a] -> a
lcmMulti forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   forall a b. (a -> b) -> [a] -> [b]
map
      (\(Integer
e,Integer
p) ->
         if Integer
p forall a. Eq a => a -> a -> Bool
== Integer
2 Bool -> Bool -> Bool
&& Integer
e forall a. Ord a => a -> a -> Bool
> Integer
2
           then Integer
pforall a. C a => a -> Integer -> a
^(Integer
eforall a. C a => a -> a -> a
-Integer
2)
           else Integer
pforall a. C a => a -> Integer -> a
^(Integer
eforall a. C a => a -> a -> a
-Integer
1) forall a. C a => a -> a -> a
* (Integer
pforall a. C a => a -> a -> a
-Integer
1)) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   forall a b. (a -> b) -> [a] -> [b]
map (forall a c b. (a -> c) -> (a, b) -> (c, b)
mapFst forall a b. (C a, C b) => a -> b
fromIntegral) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   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 =
   forall a b c. (a -> b -> c) -> b -> a -> c
flip forall a b. (a -> b) -> [a] -> [b]
map [Integer
1..] forall a b. (a -> b) -> a -> b
$ \Integer
modu ->
   let maxOrder :: Order
maxOrder = forall a. PrimitiveRoot a => a -> Order
maximumOrderOfPrimitiveRootsOfUnity (Integer
modu::Integer)
   in  forall a b. (a -> b) -> [a] -> [b]
map (forall (t :: * -> *) a. Foldable t => t a -> Int
length forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (PrimitiveRoot a, Eq a) => a -> Order -> [a]
primitiveRootsOfUnityPower Integer
modu) forall a b. (a -> b) -> a -> b
$
--       filter (flip divides maxOrder) $
       [Integer -> Order
Order Integer
1 .. Order
maxOrder]

ordersOfRootsOfUnityInteger :: [[Int]]
ordersOfRootsOfUnityInteger :: [[Int]]
ordersOfRootsOfUnityInteger =
   forall a b c. (a -> b -> c) -> b -> a -> c
flip forall a b. (a -> b) -> [a] -> [b]
map [Integer
1..] forall a b. (a -> b) -> a -> b
$ \Integer
modu ->
   forall a b. (a -> b) -> [a] -> [b]
map (forall (t :: * -> *) a. Foldable t => t a -> Int
length forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (PrimitiveRoot a, Eq a) => a -> Order -> [a]
rootsOfUnityPower (Integer
modu::Integer)) 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 =
   forall a b c. (a -> b -> c) -> b -> a -> c
flip forall a b. (a -> b) -> [a] -> [b]
map [Integer
1..] forall a b. (a -> b) -> a -> b
$ \Integer
modu ->
   let maxOrder :: Order
maxOrder = forall a. PrimitiveRoot a => a -> Order
maximumOrderOfPrimitiveRootsOfUnity (Integer
modu::Integer)
   in  forall a b. (a -> b) -> [a] -> [b]
map (forall (t :: * -> *) a. Foldable t => t a -> Int
length forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (PrimitiveRoot a, Eq a) => a -> Order -> [a]
rootsOfUnityPower Integer
modu) 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) =
   forall a. (a -> Bool) -> [a] -> [a]
filter
      (\a
n ->
         forall a. C a => a -> a -> Bool
PID.coprime a
n a
modu
         Bool -> Bool -> Bool
&&
         forall a. T a -> a
RC.representative ((a
n forall a. C a => a -> a -> T a
/: a
modu) forall a. C a => a -> Integer -> a
^ Integer
expo) forall a. Eq a => a -> a -> Bool
== forall a. C a => a
one) forall a b. (a -> b) -> a -> b
$
   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) =
   forall a. (a -> Bool) -> [a] -> [a]
filter (forall a b c. (a -> b -> c) -> b -> a -> c
flip Integer -> Order -> Bool
hasPrimitiveRootOfUnityInteger Order
order) forall a b. (a -> b) -> a -> b
$
   forall a. (a -> a) -> a -> [a]
iterate (Integer
kforall a. C a => a -> a -> a
+) Integer
1


ringsWithPrimitiveRootsOfUnityAndUnitsNaive :: [Order] -> [Integer] -> [Integer]
ringsWithPrimitiveRootsOfUnityAndUnitsNaive :: [Order] -> [Integer] -> [Integer]
ringsWithPrimitiveRootsOfUnityAndUnitsNaive [Order]
rootOrders [Integer]
units =
   forall a. (a -> Bool) -> [a] -> [a]
filter
      (\Integer
n ->
         forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Integer -> Order -> Bool
hasPrimitiveRootOfUnityInteger Integer
n) [Order]
rootOrders Bool -> Bool -> Bool
&&
         forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (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 = forall a. C a => [a] -> a
lcmMulti [Integer]
units
   in  forall a. C a => [a] -> a
lcmMulti forall a b. (a -> b) -> a -> b
$
       forall a b. (a -> b) -> [a] -> [b]
map (forall a. [a] -> a
head forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (a -> Bool) -> [a] -> [a]
filter (forall a. C a => a -> a -> Bool
PID.coprime Integer
p) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
            Order -> [Integer]
ringsWithPrimitiveRootOfUnityAndUnit) 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
_ ->
         forall a. C a => [a] -> a
product forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry forall a b. (C a, C b) => b -> a -> a
ringPower) forall a b. (a -> b) -> a -> b
$ forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$
         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 forall k a. Ord k => a -> k -> Map k a -> a
Map.findWithDefault Int
0 Integer
p Map Integer Int
factors forall a. Ord a => a -> a -> Bool
>= Int
e
                 then (Map Integer Int
factors, (Int
0,Integer
p))
                 else
                   if Integer
pforall 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
eforall a. C a => a -> a -> a
+Int
2, Integer
2))
                     else
                       (forall k a. Ord k => (a -> a -> a) -> Map k a -> Map k a -> Map k a
Map.unionWith forall a. Ord a => a -> a -> a
max Map Integer Int
factors forall a b. (a -> b) -> a -> b
$
                           forall k a. Ord k => [(k, a)] -> Map k a
Map.fromList forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall a b. (a, b) -> (b, a)
swap forall a b. (a -> b) -> a -> b
$ forall a. (PrimitiveRoot a, Ord a) => a -> [(Int, a)]
primeFactors forall a b. (a -> b) -> a -> b
$ Integer
pforall a. C a => a -> a -> a
-Integer
1,
                        (Int
eforall a. C a => a -> a -> a
+Int
1, Integer
p)))
            forall k a. Map k a
Map.empty forall a b. (a -> b) -> a -> b
$
         forall a. [a] -> [a]
reverse forall a b. (a -> b) -> a -> b
$ forall a. (PrimitiveRoot a, Ord a) => a -> [(Int, a)]
primeFactors forall a b. (a -> b) -> a -> b
$ forall a. C a => [a] -> a
lcmMulti forall a b. (a -> b) -> a -> b
$
         Integer
n forall a. a -> [a] -> [a]
: forall a b. (a -> b) -> [a] -> [b]
map (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 = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl forall a. C a => a -> a -> a
lcm 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 = forall a. (Ord a, C a) => a -> [a] -> [a]
mergePowers Integer
3 forall a b. (a -> b) -> a -> b
$ forall a. (a -> a) -> a -> [a]
iterate (Integer
2forall 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 forall a. a -> [a] -> [a]
: forall a. (a -> a -> Bool) -> [a] -> [a] -> [a]
ListHT.mergeBy forall a. Ord a => a -> a -> Bool
(<=) [a]
xs (forall a b. (a -> b) -> [a] -> [b]
map (a
pforall a. C a => a -> a -> a
*) [a]
ys)
   in  [a]
ys

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

numbers3SmoothSet :: [Integer]
numbers3SmoothSet :: [Integer]
numbers3SmoothSet =
   forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr
      (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\(Integer
m,Set Integer
rest) -> (Integer
m, forall a. Ord a => Set a -> Set a -> Set a
Set.union Set Integer
rest forall a b. (a -> b) -> a -> b
$ forall a. Eq a => [a] -> Set a
Set.fromAscList [Integer
2forall a. C a => a -> a -> a
*Integer
m,Integer
3forall a. C a => a -> a -> a
*Integer
m])) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       forall a. Set a -> Maybe (a, Set a)
Set.minView) forall a b. (a -> b) -> a -> b
$
   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
          forall a. (Ord a, C a) => a -> [a] -> [a]
mergePowers Integer
5 forall a b. (a -> b) -> a -> b
$ [Integer]
numbers3SmoothCorec
     else forall a. (Ord a, C a) => a -> [a] -> [a]
mergePowers Integer
5 forall a b. (a -> b) -> a -> b
$ forall a. (Ord a, C a) => a -> [a] -> [a]
mergePowers Integer
3 forall a b. (a -> b) -> a -> b
$ forall a. (a -> a) -> a -> [a]
iterate (Integer
2forall a. C a => a -> a -> a
*) Integer
1

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

numbers5SmoothSet :: [Integer]
numbers5SmoothSet :: [Integer]
numbers5SmoothSet =
   forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr
      (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\(Integer
m,Set Integer
rest) -> (Integer
m, forall a. Ord a => Set a -> Set a -> Set a
Set.union Set Integer
rest forall a b. (a -> b) -> a -> b
$ forall a. Eq a => [a] -> Set a
Set.fromAscList [Integer
2forall a. C a => a -> a -> a
*Integer
m,Integer
3forall a. C a => a -> a -> a
*Integer
m,Integer
5forall a. C a => a -> a -> a
*Integer
m])) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       forall a. Set a -> Maybe (a, Set a)
Set.minView) forall a b. (a -> b) -> a -> b
$
   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
1forall a. C a => a -> a -> a
+) forall a b. (a -> b) -> a -> b
$ forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall a. [a] -> a
head forall a b. (a -> b) -> a -> b
$
   forall a. (a -> Bool) -> [a] -> [a]
dropWhile (forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry forall a. Eq a => a -> a -> Bool
(/=)) forall a b. (a -> b) -> a -> b
$
   forall a b. (a -> a -> b) -> [a] -> [b]
ListHT.mapAdjacent (,) forall a b. (a -> b) -> a -> b
$
   forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\a
m Int
d -> forall a. Bits a => a -> Int -> a
shiftR a
m Int
d forall a. Bits a => a -> a -> a
.|. a
m) (a
nforall a. C a => a -> a -> a
-a
1) forall a b. (a -> b) -> a -> b
$
   forall a. (a -> a) -> a -> [a]
iterate (Int
2forall 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 forall a. C a => a -> Integer -> a
^ forall a b. (C a, C b) => a -> b
fromIntegral (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 =
   forall (t :: * -> *) a. Foldable t => t a -> Int
length forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (a -> Bool) -> [a] -> [a]
takeWhile (forall a. Ord a => a -> a -> Bool
>a
0) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (a -> a) -> a -> [a]
iterate (forall a b c. (a -> b -> c) -> b -> a -> c
flip forall a. C a => a -> a -> a
div a
base) forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 =
   forall a. [a] -> a
last forall a b. (a -> b) -> a -> b
$
   a
n forall a. a -> [a] -> [a]
: forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr (\a
m -> case forall a. C a => a -> a -> (a, a)
divMod a
m a
b of (a
q,a
r) -> forall a. Bool -> a -> Maybe a
toMaybe (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 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 =
   forall a. [a] -> a
last forall a b. (a -> b) -> a -> b
$ (Int
0,a
n) forall a. a -> [a] -> [a]
:
   forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr
      (\(Int
e,a
m) ->
         let (a
q,a
r) = forall a. C a => a -> a -> (a, a)
divMod a
m a
b
             eq :: (Int, a)
eq = (Int
eforall a. C a => a -> a -> a
+Int
1,a
q)
         in  forall a. Bool -> a -> Maybe a
toMaybe (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
1forall a. Eq a => a -> a -> Bool
==) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   forall a. (C a, C a) => a -> a -> a
divideByMaximumPower Integer
3 forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   forall a. (C a, C a) => a -> a -> a
divideByMaximumPower Integer
2

is5Smooth :: Integer -> Bool
is5Smooth :: Integer -> Bool
is5Smooth =
   (Integer
1forall a. Eq a => a -> a -> Bool
==) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   forall a. (C a, C a) => a -> a -> a
divideByMaximumPower Integer
5 forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   forall a. (C a, C a) => a -> a -> a
divideByMaximumPower Integer
3 forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   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 =
   forall a. [a] -> a
head forall a b. (a -> b) -> a -> b
$ forall a. (a -> Bool) -> [a] -> [a]
dropWhile (forall a. Ord a => a -> a -> Bool
<Integer
n) [Integer]
numbers3Smooth

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

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

ceiling5SmoothNaive :: Integer -> Integer
ceiling5SmoothNaive :: Integer -> Integer
ceiling5SmoothNaive =
   forall a. [a] -> a
head forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (a -> Bool) -> [a] -> [a]
dropWhile (Bool -> Bool
not forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Bool
is5Smooth) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (a -> a) -> a -> [a]
iterate (Integer
1forall 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 =
   forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum forall a b. (a -> b) -> a -> b
$ Integer -> Integer -> Integer -> Integer -> [Integer]
ceilingSmoothsTrace Integer
2 Integer
3 Integer
n forall a b. (a -> b) -> a -> b
$ 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 =
   forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Integer -> Integer -> Integer -> [Integer]
ceilingSmoothsTrace Integer
3 Integer
5 Integer
n) forall a b. (a -> b) -> a -> b
$
   Integer -> Integer -> Integer -> Integer -> [Integer]
ceilingSmoothsTrace Integer
2 Integer
3 Integer
n forall a b. (a -> b) -> a -> b
$ 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 =
   forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Integer -> Integer -> Integer -> [Integer]
ceilingSmoothsTrace Integer
2 Integer
5 Integer
n) forall a b. (a -> b) -> a -> b
$
   Integer -> Integer -> Integer -> Integer -> [Integer]
ceilingSmoothsTrace Integer
2 Integer
3 Integer
n forall a b. (a -> b) -> a -> b
$ 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 forall a. C a => a -> a -> (a, a)
divMod Integer
k Integer
a of
             (Integer
q,Integer
r) -> if Integer
rforall a. Eq a => a -> a -> Bool
==Integer
0 Bool -> Bool -> Bool
&& Integer
qforall 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 forall a. a -> [a] -> [a]
: if forall a. C a => a -> a -> a
mod Integer
m Integer
a forall a. Eq a => a -> a -> Bool
== Integer
0 then Integer -> [Integer]
go forall a b. (a -> b) -> a -> b
$ Integer -> Integer
divMany forall a b. (a -> b) -> a -> b
$ Integer
mforall 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 = forall a. Eq a => [a] -> Set a
Set.fromAscList forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (C a, Bits a, C a, Ord a) => a -> [a]
uniquePrimeFactors
   in  forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr
         (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap
             (\(Integer
p,Set Integer
set) ->
                (Integer
p, forall a. Ord a => Set a -> Set a -> Set a
Set.union (Integer -> Set Integer
primeFactorSet (Integer
pforall a. C a => a -> a -> a
-Integer
1)) Set Integer
set)) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
          forall a. Set a -> Maybe (a, Set a)
Set.maxView)
       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 =
          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) = forall a. C a => a -> a -> (a, a)
divMod a
m a
p
                in  if a
rforall a. Eq a => a -> a -> Bool
==a
0
                      then a
p forall a. a -> [a] -> [a]
: a -> [a]
go (forall a. (C a, C a) => a -> a -> a
divideByMaximumPower a
p a
q)
                      else
                        if a
q forall a. Ord a => a -> a -> Bool
>= a
p
                          then a -> [a]
go a
m
                          else if a
mforall a. Eq a => a -> a -> Bool
==a
1 then [] else a
m forall a. a -> [a] -> [a]
: [])
             (forall a. HasCallStack => [Char] -> a
error [Char]
"uniquePrimeFactors: end of infinite list")
             (forall a. (a -> a) -> a -> [a]
iterate (a
2forall a. C a => a -> a -> a
+) a
3)
   in  case 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 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 =
          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) = forall a. C a => a -> a -> (a, a)
divMod a
m a
p
                in  if a
rforall a. Eq a => a -> a -> Bool
==a
0
                      then
                        case forall a. (C a, C a) => a -> a -> (Int, a)
getMaximumExponent a
p a
q0 of
                          (Int
e,a
q1) -> (Int
eforall a. C a => a -> a -> a
+Int
1,a
p) forall a. a -> [a] -> [a]
: a -> [(Int, a)]
go a
q1
                      else
                        if a
q0 forall a. Ord a => a -> a -> Bool
>= a
p
                          then a -> [(Int, a)]
go a
m
                          else if a
mforall a. Eq a => a -> a -> Bool
==a
1 then [] else (Int
1,a
m) forall a. a -> [a] -> [a]
: [])
             (forall a b. a -> b -> a
const [])
             (forall a. (a -> Bool) -> [a] -> [a]
filter (Bool -> Bool
not forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. C a => a -> Bool
Units.isUnit) forall a b. (a -> b) -> a -> b
$
              forall a. PrimitiveRoot a => a -> [a]
primitiveRootCandidates a
n)
   in  case 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) 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 forall a. (PrimitiveRoot a, Ord a) => a -> [(Int, a)]
primeFactors Integer
n of
      [] -> Bool
False
      (Int
e,Integer
m):[(Int, Integer)]
_ -> Int
eforall a. Eq a => a -> a -> Bool
==Int
1 Bool -> Bool -> Bool
&& Integer
mforall 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 =
   forall a. (a -> a) -> a -> [a]
iterate
      (\Integer
n ->
         forall a. [a] -> a
head forall a b. (a -> b) -> a -> b
$ forall a. (a -> Bool) -> [a] -> [a]
dropWhile (Bool -> Bool
not forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Bool
isPrime) forall a b. (a -> b) -> a -> b
$
         forall a. [a] -> [a]
tail forall a b. (a -> b) -> a -> b
$ forall a. (a -> a) -> a -> [a]
iterate (Integer
nforall 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 = forall a b. (C a, C b) => a -> b
fromIntegral Int
n
             ps :: [Integer]
ps = forall a. (a -> Bool) -> [a] -> [a]
filter (forall a. Ord a => a -> a -> Bool
>Integer
1) (forall a b. (a -> b) -> [a] -> [b]
map (forall a. C a => a -> a -> a
subtract Integer
1) (Integer -> [Integer]
partialPrimes Integer
ni))
         in  [Order] -> [Integer] -> Integer
ringWithPrimitiveRootsOfUnityAndUnits (forall a b. (a -> b) -> [a] -> [b]
map Integer -> Order
Order forall a b. (a -> b) -> a -> b
$ Integer
ni 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 = 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.
             -}
             forall a. [a] -> a
head forall a b. (a -> b) -> a -> b
$ forall a. (a -> Bool) -> [a] -> [a]
filter Integer -> Bool
isPrime forall a b. (a -> b) -> a -> b
$
             (\Integer
order -> forall a. (a -> a) -> a -> [a]
iterate (Integer
order forall a. C a => a -> a -> a
+) Integer
1) forall a b. (a -> b) -> a -> b
$
             forall a. C a => [a] -> a
lcmMulti forall a b. (a -> b) -> a -> b
$
             Integer
ni forall a. a -> [a] -> [a]
: forall a b. (a -> b) -> [a] -> [b]
map (forall a. C a => a -> a -> a
subtract Integer
1) (Integer -> [Integer]
partialPrimes Integer
ni)