module MPS.MathExt where

import Data.List
import Data.Ratio
import MPS.Hack.Dot
import MPS.Math.CombinatoricsGeneration
import MPS.Math.NumberTheoryFundamentals
import qualified MPS.Math.PermutationGroups as Group
import MPS.Math.Primes
import MPS.Snippets

import Prelude hiding ((.), (^))

phi           = (sqrt(5)+1)/2
fib n         = (( phi^n - (1-phi)^n ) / sqrt(5)).round
fibs          = 1 : 1 : zipWith (+) fibs (fibs.tail)
primes        = primesTo10000 ++ [10007, 10009 ..].filter(isPrime)

m_indices m n = [ (i,j) | i <- [0..m-1], j <- [0..n-1] ]
at2 (i,j) m   = m.at(i).at(j)

factors       = primePowerFactors
divisors n    = n.factors.map(\(p,q) -> [ p ^ i | i <- [0..q] ]).cartProd.map(product).sort.init


perm_divisors 1 = [[]]
perm_divisors n = [i:xs | 
  i <- (n: n.divisors.tail.reverse),
  let j = n `div` i,
  let xss = perm_divisors j,
  xs <- xss 
  ].map sort .unique .rsort

totient m = m.factors.map phi_token.product
  where phi_token (p, l) = (p-1)*p^(l-1)

combinations_for 0 xs = [[]]
combinations_for n xs = [y:ys | y <- xs, ys <- combinations_for (n - 1) xs]

is_prime              = isPrime

power_permutations xs = xs.powerset.concatMap(permutationsOf)

is_proper_integer x   = x.floor .is ( x.ceiling )
is_square x           = x.fromIntegral.sqrt.is_proper_integer
is_pythagorean a b c  = a^2 + b^2 == c^2



cf_iter x y = x + 1 / y
cf xs = xs.map(%1) .foldr1 cf_iter

cf_e = 2: ([[1,n,1] | n <- [2,4..]] .join')

-- farey sequence
median_q x y = ( x.numerator + y.numerator ) % ( x.denominator + y.denominator )
farey_iter xs = xs
  .zip (xs.tail)
  .map (splash median_q)
  .zip (xs.tail)
  .map swap
  .map list2
  .join'
  .(xs.first : )

fareys' xs 1 = xs
fareys' xs n = fareys' (farey_iter xs) (n-1)
fareys n     = fareys' [0,1] n

bounded_fareys' xs 1 _ = xs
bounded_fareys' xs n c = bounded_fareys' filtered (n-1) (c+1) where
  filtered             = farey_iter xs
    .reject (\q -> q.denominator > c)


bounded_fareys n = bounded_fareys' [0,1] n 2 .tail .init

-- pythagoreans generators
pythagoreans = [ (s, (a, b, c), (m, n) ) |
  m <- [2..],
  n <- [1..m-1],
  (m + n).odd,
  gcd m n == 1, 
  let a = m^2 - n ^2,
  let b = 2 * m * n,
  let c = m^2 + n ^ 2,
  let s = a + b + c
  ]

pythagoreans' x = pythagoreans
  .ub(\(_, _, (m,_)) -> m <= bound)
  .filter(\(s, _, _) -> s <= x)
  where bound = x.from_i.sqrt.floor


turns xs            = xs.iterate Group.rotateL .take(xs.length)