{-# LANGUAGE FlexibleInstances, UndecidableInstances #-}
-- | Prufer domains are non-Noetherian analogues of Dedekind domains. That is
-- integral domains in which every finitely generated ideal is invertible. This 
-- implementation is mainly based on:
--
-- http:\/\/hlombardi.free.fr\/liens\/salouThesis.pdf
--
module Algebra.Structures.PruferDomain 
  ( PruferDomain(..)
  , propCalcUVW, propPruferDomain
  , calcUVW_B, calcUVWT, propCalcUVWT, fromUVWTtoUVW
  , computePLM_PD
  , invertIdeal
  , intersectionPD, intersectionPDWitness, solvePD
  ) where

import Test.QuickCheck
import Data.List (nub, (\\))

import Algebra.Structures.IntegralDomain
import Algebra.Structures.BezoutDomain
import Algebra.Structures.Coherent
import Algebra.Ideal
import Algebra.Matrix


-------------------------------------------------------------------------------
-- | Given a and b it computes u, v and we such that:
-- 
--  (1) au = bv
--
--  (2) b(1-u) = aw
--
class IntegralDomain a => PruferDomain a where
  --         a    b     u v w
  calcUVW :: a -> a -> (a,a,a)

-- | Property specifying that:
-- au = bv and b(1-u) = aw
propCalcUVW :: (PruferDomain a, Eq a) => a -> a -> Bool
propCalcUVW a b = a <*> u == b <*> v && b <*> (one <-> u) == a <*> w
  where (u,v,w) = calcUVW a b

propPruferDomain :: (PruferDomain a, Eq a) => a -> a -> a -> Property
propPruferDomain a b c | propCalcUVW a b = propIntegralDomain a b c
                       | otherwise       = whenFail (print "propCalcUVW") False


-- | Alternative characterization of Prufer domains, given a and b compute u, v, 
-- w, t such that:
-- 
-- ua = vb && wa  = tb && u+t = 1
calcUVWT :: PruferDomain a => a -> a -> (a,a,a,a)
calcUVWT a b = (x,y,z,one <-> x)
  where (x,y,z) = calcUVW a b

propCalcUVWT :: (PruferDomain a, Eq a) => a -> a -> Bool
propCalcUVWT a b = u <*> a == v <*> b && w <*> a == t <*> b && u <+> t == one
  where (u,v,w,t) = calcUVWT a b

-- | Go back to the original definition.
fromUVWTtoUVW :: (a,a,a,a) -> (a,a,a)
fromUVWTtoUVW (u,v,w,t) = (u,v,w) 


--------------------------------------------------------------------------------
-- | Bezout domain -> Prufer domain
--
{-
Prufer: forall a b exists u v w t.  u+t = 1 &  ua = vb & wa = tb

We consider only domain.
We assume we have the Bezout condition: given a, b we can find g,a1,b1,c,d s.t.

a = g a1
b = g b1
1 = c a1 + d b1

We try then 

u = d b1
t = c a1

We should find v such that
a d b1 = b v
this simplifies to 
g a1 d b1 = g b1 v
and we can take 
v = a1 d
Similarly we can take 
w = b1 c

We have shown that Bezout domain -> Prufer domain.

instance (BezoutDomain a, Eq a) => PruferDomain a where
-}

-- | Proof that all Bezout domains are Prufer domains.
calcUVW_B :: (BezoutDomain a, Eq a) => a -> a -> (a,a,a)
calcUVW_B a b | a == zero = (one,zero,zero)
              | b == zero = (zero,zero,zero)
              | otherwise = fromUVWTtoUVW (u,v,w,t)
    where
    -- Compute g, a1 and b1 such that:
    -- a = g*a1
    -- b = g*b1
    (g,[_,_],[a1,b1])  = toPrincipal (Id [a,b])
    
    -- Compute c and d such that:
    -- 1 = a1*c + a2*d
    (_,[c,d],_) = toPrincipal (Id [a1,b1])

    u = d <*> b1
    t = c <*> a1
    v = d <*> a1
    w = c <*> b1

-------------------------------------------------------------------------------
-- Coherence

-- | Compute a principal localization matrix for an ideal in a Prufer domain.
computePLM_PD :: (PruferDomain a, Eq a) => Ideal a -> Matrix a
computePLM_PD (Id [_])   = matrix [[one]]
computePLM_PD (Id [a,b]) = let (u,v,w,t) = calcUVWT b a 
                           in M [ Vec [u,v], Vec [w,t]]
computePLM_PD (Id xs)    = matrix a
  where
  -- Use induction hypothesis to construct a matrix for n-1:
  x_is = init xs
  b    = unMVec $ computePLM_PD (Id x_is)
  m    = length b - 1
  
  -- Let s_i be b_ii:
  s_is = [ (b !! i) !! i | i <- [0..m]]

  -- Take out x_n:
  x_n  = last xs

  -- Compute (u_i, v_i, w_i, t_i) for <x_n,x_i>:
  uvwt_i = [ calcUVWT x_n x_i | x_i <- x_is ]
    
  -- Take out all u, v, w, and t:
  u_is = [ u_i | (u_i,_,_,_) <- uvwt_i ]
  v_is = [ v_i | (_,v_i,_,_) <- uvwt_i ]
  w_js = [ w_i | (_,_,w_i,_) <- uvwt_i ]
  t_is = [ t_i | (_,_,_,t_i) <- uvwt_i ]
  
  -- COMPUTE a_ij for 1 <= i,j < n
  -- i = row
  -- j = column
  a_ij = [ [ if i == j 
                then (s_is !! i) <*> (u_is !! i)
                else (u_is !! i) <*> (b !! i !! j)
           | j <- [0..m] ]
         | i <- [0..m] ]

  -- COMPUTE a_nn
  a_nn = sumRing $ zipWith (<*>) s_is t_is

  -- COMPUTE a_ni for 1 <= i < n
  -- THIS IS THE LAST ROW
  a_ni = [ sumRing [ (b !! j !! i) <*> (w_js !! j)
                   | j <- [0..m] ]
         | i <- [0..m] ]

  -- COMPUTE a_in for 1 <= i < n
  -- THIS IS THE LAST COLUMN
  a_in = [ (s_is !! i) <*> (v_is !! i)
         | i <- [0..m] ]

  -- ASSEMBLE EVERYTHING
  a = [ x ++ [y] | (x,y) <- zip a_ij a_in ] ++ [a_ni ++ [a_nn]]


-- | Ideal inversion. Given I compute J such that IJ is principal.
-- Uses the principal localization matrix for the ideal.
invertIdeal :: (PruferDomain a, Eq a) => Ideal a -> Ideal a
invertIdeal xs = 
  let a = unMVec $ computePLM_PD xs

      -- Pick out the first column
      a_njs = [ head (a !! j) | j <- [0..length a - 1]]
  in Id a_njs

-- | Compute the intersection of I and J by:
--       
--       (I \\cap J)(I + J) = IJ  => (I \\cap J)(I + J)(I + J)' = IJ(I + J)'
--
intersectionPDWitness :: (PruferDomain a, Eq a) => Ideal a -> Ideal a -> (Ideal a,[[a]],[[a]])
intersectionPDWitness (Id is) (Id js) = (int,wis,wjs)
  where
  lj  = length js 
  li  = length is

  ij  = Id (is ++ js)

  plm = computePLM_PD ij

  as  = take li $ unMVec $ transpose plm
  as' = drop li $ unMVec $ transpose plm

  int = Id $ concat [ map (j <*>) a | j <- js , a <- as ]

  wis = concat [ [ addZ i li a | a <- as ] | as <- as', i <- [0..li-1] ]
  wjs = [ addZ i lj a | i <- [0..lj-1], a <- concat as ]

  addZ n l x = replicate n zero ++ x : replicate (l-n-1) zero

intersectionPD :: (PruferDomain a, Eq a) => Ideal a -> Ideal a -> Ideal a
intersectionPD i j = fst3 (intersectionPDWitness i j)
  where fst3 (x,_,_) = x

-- | Coherence of Prufer domains.
solvePD :: (PruferDomain a, Eq a) => Vector a -> Matrix a
solvePD x = solveWithIntersection x intersectionPDWitness

-- instance (PruferDomain a, Eq a) => Coherent a where
--   solve = solvePD