module WWN 

where 

import Data.List
import Data.Ratio
import Data.Char
import Data.Complex
import Nats

binary :: Integral a => a -> [Int]
binary x = reverse (bits x)
  where bits 0 = [0]
        bits 1 = [1]
        bits n = fromIntegral (rem n 2) : bits (quot n 2)

showDigits :: [Int] -> String
showDigits = map intToDigit

bin ::  Integral a => a -> String
bin = showDigits . binary

coprime :: Integer -> Integer -> Bool
coprime m n = (gcd m n) == 1

gcde :: Integral a => a -> a -> (a,a,a) 
gcde a b | a == 0 && b == 0 = error "gcd undefined"
         | a < 0  || b < 0  = error "gcd undefined"
         | a == 0           = (b,0,1)
         |           b == 0 = (a,1,0)
         | a > b            = let 
                               (k,r)   = quotRem a b
                               (d,x,y) = gcde b r
                              in (d, y, x-k*y)
         | otherwise        = let
                               (k,r)   = quotRem b a
                               (d,x,y) = gcde a r
                              in (d, x-k*y, y)

data Sgn = P | N deriving (Eq,Show)
type MyInt = (Sgn,Natural)

myplus :: MyInt -> MyInt -> MyInt
myplus (s1,m) (s2,n) | s1 == s2           = (s1,m+n)
                     | s1 == P  && n <= m = (P,m-n)
                     | s1 == P  && n > m  = (N,n-m)
                     | otherwise          = myplus (s2,n) (s1,m) 

type NatPair = (Natural,Natural)

plus1 :: NatPair -> NatPair -> NatPair 
plus1 (m1, m2) (n1, n2) = (m1+n1, m2+n2)

subtr1 :: NatPair -> NatPair -> NatPair 
subtr1 (m1, m2) (n1, n2) = plus1 (m1, m2) (n2, n1)

mult1 :: NatPair -> NatPair -> NatPair 
mult1 (m1, m2) (n1, n2) = (m1*n1 + m2*n2, m1*n2 + m2*n1)

eq1 :: NatPair -> NatPair -> Bool
eq1 (m1, m2) (n1, n2) = (m1+n2) == (m2+n1)

reduce1 :: NatPair -> NatPair
reduce1 (m1,Z) = (m1,Z)
reduce1 (Z,m2) = (Z,m2) 
reduce1 (S m1, S m2) = reduce1 (m1, m2)

decExpand :: Rational -> [Integer]
decExpand x | x < 0     = error "negative argument"
            | r == 0    = [q]
            | otherwise = q : decExpand ((r*10) % d)
   where 
   (q,r) = quotRem n d
   n     = numerator x 
   d     = denominator x 

decForm :: Rational -> (Integer,[Int],[Int]) 
decForm x | x < 0     = error "negative argument"
          | otherwise = (q,ys,zs) 
  where 
  (q,r)   = quotRem n d 
  n       = numerator x 
  d       = denominator x 
  (ys,zs) = decF (r*10) d []

decF :: Integer -> Integer -> [(Int,Integer)] -> ([Int],[Int])
decF n d xs | r == 0        = (reverse (q: (map fst xs)),[])
            | elem (q,r) xs = (ys,zs) 
            | otherwise     =  decF (r*10) d ((q,r):xs)
     where 
     (q',r)  = quotRem n d
     q       = fromIntegral q'
     xs'     = reverse xs
     Just k  = elemIndex (q,r) xs'
     (ys,zs) = splitAt k (map fst xs')

mechanicsRule :: Rational -> Rational -> Rational 
mechanicsRule p x = (1 % 2) * (x + (p * (recip x)))

mechanics :: Rational -> Rational -> [Rational]
mechanics p x = iterate (mechanicsRule p) x

sqrtM :: Rational -> [Rational]
sqrtM p | p < 0     = error "negative argument"
        | otherwise = mechanics p s
   where 
   s = if xs == [] then 1 else last xs
   xs = takeWhile  (\ m -> m^2 <= p) [1..] 

approximate :: Rational -> [Rational] -> Rational 
approximate eps (x:y:zs) 
  | abs (y-x) < eps = y
  | otherwise       = approximate eps (y:zs) 

apprx :: [Rational] -> Rational
apprx = approximate (1/10^6)

mySqrt :: Rational -> Rational
mySqrt p = apprx (sqrtM p)

solveQ :: (Complex Float, Complex Float, Complex Float)
                               -> (Complex Float, Complex Float) 
solveQ =  \ (a,b,c) -> if a == 0 then error "not quadratic"
                                 else let d = b^2 - 4*a*c in 
                                           ((- b + sqrt d) / 2*a,
                                            (- b - sqrt d) / 2*a)

