--------------------------------------------------------------------
--------------------------------------------------------------------
--  The Algebraic Domain Constructor  DoCon,   version 2.09
--
--  Copyright  Serge Mechveliani,    2005
--------------------------------------------------------------------
--------------------------------------------------------------------





module Polrel_ 

  -- Finding linear relation generators  for  (e-)polynomials
  -- (syzygy module generators).
  -- and  Algebraic relation generators for the polynomials.
  --
  -- All needed from here is  reexported by  GBasis.

  (polRelGens, polRelGens_e, algRelsPols)

where
import qualified Data.Map as Map (empty)

import Maybe (fromMaybe)
import List  (partition, genericLength, genericTake, genericDrop)

import DPrelude   (PropValue(..))
import SetGroup   (Set(..), zeroS, isZero, unity)
import RingModule (EuclideanRing(..), isEucRing, upRing) 
import VecMatr    (Vector(..), Matrix(..), vecRepr, scalarMt)
import PP_        (ppComp_blockwise)
import UPol_      (PPOrdTerm, PolVar, ppoComp)
import Pol_       (Pol(..), polMons, polPPComp)
import EPol0_     (EPol(..), epolPol, epolLCoord, polToEPol)
import EPol1_     (sEPVecP)
import PolNF_     (polNF_e)
import GBas_      (gBasis, gBasis_e)

import qualified GBas_ (initVecs_)





--------------------------------------------------------------------
polRelGens :: EuclideanRing a => String -> [Pol a] -> [[Pol a]]
                                 -- mode   fs         rs


{-  The generator list for the linear relation (syzygy) module 
    for a polynomial list  fs  over an Euclidean GCD-ring  a.

fs  should be a non-empty list of the elements of  P = a[x1..xn].

rs :: [[P]]  is a list of rows,   ri <- P^n.

each row  r = [r1..rn]  is a linear relation for  fs: 
                                                f1*r1+...+fn*rn = 0.

mode =  "g"   means  fs  is a weak reduced Groebner basis
                            - here the computation is more direct,
        ""           generic case.
For the weak Groebner basis  fs  over a *field*  and ecpTOP_weight
ordering [MoM] on P^n generated by the weights  lpp(fi), fi <- fs,   
the result is also a weak Groebner basis of the submodule [MoM].
(Does this also hold for the Euclidean coefficients?)

METHOD. [Bu]

1. Case  mode = "g"

It is known that in this case the "s-relations" generate the 
relation module.
Here the local function gCase  reduces the s-polynomials sPol(fi,fj)
by  fs to zeroes.  The coefficient vectors of this reduction summed 
with the complementary monomial factors for  fi,fj  form a basis for
syzygy generators for  fs.  
To build the relation for (fi,fj)  gCase  creates a copy of  fi-s
supplied with the vectors:
fvs = [fv1..fvn],  fvi = (fi,vi)   where  vi  is the  i-th  vector 
from the unity matrix over  P.
Let  m1,m2  are the factors for the s-polynomial of fi,fj,

(zeroPol,qs) = moduloBasis "g" fs (m1*fi-m2*fj)

Then,  m1*vi - m2*vj - qs   is the relation for (i,j).
 

2. Case  mode = "" 

reduces to the "g" case as follows.

1) (gs,mt) = gBasis "" fs   produces the matrix  mt  such that  
                            mt*transp(fs) = gs  is a Groebner basis.
2) map (moduloBasis "g" gs) fs   
                       produces  mt1  such that  mt1*transp(gs) = fs
3) rels1  = mt1*mt - unityMatrix                  
4) relsG  = gCase gs             makes the basic relations for  gs
5) result = rels1 ++ (relsG*mt)   
                   - the second matrix is set "under" the first one.

  The matrices are multiplied over  Pol a. 

This method founds on the following consideration (hint):

   relsG*mt   are the ported relations for the Groebner basis.
   rels1      are the relations caused by the fact that   mt   and
              mt1  may be not precisely  mutually  inverse  module 
              homomorphisms.   

This all is valid for the  e-polynomials -  and  polRelGens  calls 
polRelGens_e.
Note also that in the "g" case for the  e-polynomials   the 
s-relations have to form only for  fvi, fvj  having the same leading
coordinate. 
-}
--------------------------------------------------------------------   


polRelGens  mode  []       = error ("polRelGens "++mode++" []\n")
polRelGens  mode  fs@(f:_) =
                         let cp                = polPPComp f
                             ecp (_, p) (_, q) = cp p q
                             eo                = (ecp, "a", [], cp)
                         in  
                         polRelGens_e mode $ map (polToEPol 1 eo) fs



--------------------------------------------------------------------
polRelGens_e :: EuclideanRing a => String -> [EPol a] -> [[Pol a]]
polRelGens_e                       mode      fs       =  
  let
    pol@(Pol _ c _ _ dC) = epolPol $ head fs
    (_, eC)              = baseEucRing c dC
    (zeroPol, unPol)     = (zeroS pol, unity pol)
    isZeroV              = all (== zeroPol)
    dP                   = upRing pol Map.empty
                                   -- for matrices over polynomials,
                                   -- - somewhat dummy 

    msgFs = (("polRelGens_e "++mode++" fs,")++) . 
            ("\nlength fs =  "++) . shows (genericLength fs) .
            ("\nhead   fs =  "++) . shows (head fs) . 
            ("\n <-  "        ++) . showsDomOf (head fs) . 
            ("\n\n"++)
    ----------------------------------------------------------------
    -- forming the relation on  gs  from the pair  (fv,gv),  
    -- (f,u) = fv,  (g,v) = gv,  gs  a non-empty weak Groebner basis

    relForPair gs (fv, gv) =       
      let 
        ((s, w), _, _) =                              -- (_, m1, m2)
             fromMaybe (error $ msgFs $ msg 
                                  ("... sEPVecP fv gv -> Nothing "++
                                   "  for some  fv, gv    -?\n"
                                  ) 
                       ) $ sEPVecP fv gv
                                   -- sPol for vector-ed polynomials
                                   -- s = m1*f-m2*g
                                   -- w = m1*u-m2*v
        msg = ("... gCase gs,"++) . 
              ("\nlength gs =  "++) . shows (genericLength gs) . 
              ("\nhead   gs =  "++) . shows (head gs) . ("\n\n"++)
 
        (r, qs) = polNF_e "" gs s
                          -- vi is the canonical unit vector for fi.
                          -- We have:  m1*f-m2*g = qs.gs;   hence
                          -- (m1*u-m2*v-qs).gs = 0,  and  (w-qs)  is
                          -- the needed (i,j) relation for  gs.
      in
      if  isZero r  then  zipWith (-) w qs
      else
      error $ msgFs $ msg
                       ("Some s-polynomial does not reduce to zero." 
                        ++"\nMaybe, the mode was set wrongly\n"
                       )
    ----------------------------------------------------------------
    gCase fs =         -- case fs is a non-empty weak Groebner basis   
      let
        fvs = GBas_.initVecs_ unPol fs
                    -- we need a list  pairs  consisting of (fvi,fvj)

        splitByCoord []            = []
        splitByCoord ((f, v): fvs) = 
            let
              o          = epolLCoord f
              (gvs, hvs) = partition ((== o) . epolLCoord . fst) fvs
            in  
            ((f, v): gvs): (splitByCoord hvs)

        blocks = splitByCoord fvs

        critPairs []       = []
        critPairs (fv:fvs) = [(fv, x)| x <- fvs] ++ (critPairs fvs)
      in
      map (relForPair fs) $ concat $ map critPairs blocks
    ----------------------------------------------------------------
                       -- Reduction to the Groebner basis case. 
    genericCase fs =   -- Here  fs  is non-empty and free of zeroes.
      let
        {(gs, mt) = gBasis_e fs;  fBy_g_s = map (polNF_e "" gs) fs}
      in 
      if  any (not . isZero . fst) fBy_g_s  then
                  error $
                  msgFs ("genericCase:   gs = ..gBasis_e fs,  ..."++
                         "map (polNF_e \"\" gs) fs:  "++
                         "\n\nNon-zero remainder appeared  - ? \n"
                        )
      else
      let mt'             = Mt mt dP
          mt1'            = Mt (map snd fBy_g_s) dP
          m1m'@(Mt m1m _) = mt1'*mt'
          unityMt'        = Mt (scalarMt m1m unPol zeroPol) dP
          Mt rels1 _      = m1m' - unityMt'
          relsG'          = Mt (gCase gs) dP
          Mt rsG_m _      = relsG' * mt'
      in
      filter (not . isZeroV) (rels1 ++ rsG_m)
    ----------------------------------------------------------------
  in
  case (fs, all isZero fs, isEucRing eC == Yes, mode)
  of
  ([], _   , _    , _  ) -> error ("polRelGens_e "++mode++" []\n")
  (_ , True, _    , _  ) -> scalarMt fs unPol zeroPol
  (_ , _   , False, _  ) -> 
                          error $ msgFs
                          "Euclidean ring needed for coefficients\n"

  (_ , _   , _    , "g") -> gCase fs 
  _                      ->
    let         
      -- Reduce to the case of  fs  free of zeroes.
      -- First, numerate fs and part zeroes ...

      m               = genericLength fs
      nfs             = zip [1 .. m] fs
      (zNfs, nonzNfs) = partition (isZero . snd) nfs
      zInds           = map fst zNfs
      (nonzInds, fs') = unzip nonzNfs    -- fs' are non-zeroes
      relsForNonZero  = genericCase fs'
      relsForZero     =  
                    if  null zNfs  then  []
                    else                 scalarMt zNfs unPol zeroPol
                                                     -- unity matrix
           
              -- Now  relsForNonZero,relsForZero  extend with zeroes 
              -- in each vector (in the appropriate positions), and 
              -- the two obtained matrices concatenate.

                         -- complete each vector from relsForNonZero 
                         -- with zeroes set in the rest positions

      nRelsNz = map (zip nonzInds) relsForNonZero
      zNfs_p  = [(i, zeroPol) | (i, _) <- zNfs]
      relsForNonZero_completed =  map (merge zNfs_p) nRelsNz

      merge []                nfs               = map snd nfs
      merge nzs               []                = map snd nzs
      merge nzs@((i,z): nzs') nfs@((j,f): nfs') =  
                                 if i < j then  z: (merge nzs' nfs )
                                 else           f: (merge nzs  nfs')
      nRelsZ  = map (zip zInds) relsForZero
      nzNfs_p = [(i, zeroPol) | (i, _) <- nonzNfs]
      relsForZero_completed = map (merge nzNfs_p) nRelsZ
    in
    relsForNonZero_completed ++ relsForZero_completed





{- RESERVE *******************************************************
This old version yields the smaller list of the syzygy generators. 
It uses for the "g" case the *criteria* for the critical pairs.
But we doubt whether in this case it will be a Groebner basis
relations for  ecpTOP_weights  ordering.
    repeatUpdatePairs_field = repeatUpdatePairs_ev     False
    repeatUpdatePairs_euc   = repeatUpdatePairs_ev_euc False
        -- here `False' is for skipping the  ppMutPrime  criterion
    --------------------------------------------------------------
    gCase fs =       -- case fs is a non-empty weak Groebner basis   
      let  r    = genericLength fs
        fvs  = GBas_.initVecs_ unityPol fs
                  -- we need a list  pairs  consisting of (fvi,fvj)
        fvsM = listToFM  (zipWith  (\i f->(i,(f,False)))  [1..r] fvs)
        bringIn (fv,_) = [(fst (f_No fvsM i), fv)| (i,_) <- ? ]
      in if  prop=="f"  then  
          let indBlocks = repeatUpdatePairs_field [] fvsM [2..r]
              iRows     = concat (map fst indBlocks)
              rs      = [bringIn (f_No fvsM j) row| (row,j)<-iRows]
          in  map (relForPair fs) (concat rs)
        else
          let indBlocks = repeatUpdatePairs_euc [] fvsM [2..r]
              iRows     = concat (map fst indBlocks)
              rs = [bringIn (f_No fvsM j) row | (row,j)<- iRows]
          in map (relForPair fs) (concat rs)
END RESERVE ********************************************************
-}







--------------------------------------------------------------------
algRelsPols :: 
              EuclideanRing k
              => 
              [Pol k] -> [PolVar] -> PPOrdTerm -> [Pol k] -> [Pol k]
              -- hs      ys          oY           fs         rels

  -- Ideal generators for the algebraic relations for polynomials fs
  -- considered modulo Ideal(hs).
  -- f(i),h(i) are from  A = k[x1,...,xn].
  -- k  a field   (does this work for any Euclidean k too?)
  --
  -- Any  rel <- rels   is from  B = k[y1..ym],
  --
  -- yi    form the list  ys  of the variables corresponding 
  --       bijectively to fi,  so that  rel(f1,...,fm) = 0  in A,
  -- rels  is the reduced Groebner basis in B for the ideal of all 
  --       such algebraic relations for  fs.
  -- oY    is the power product ordering description chosen for B.
  --
  -- No  fi  should belong to  Ideal(hs).
  --
  -- METHOD
  -- is taken from the paper [GTZ]:
  --
  -- 1. Embed  fi,hi  to  C = k[x1..xn,y1..ym] = k[X,Y],   
  --
  -- 2. Let  pi = fi - yi  in  C,   
  --    C  is viewed under the direct-sum power product comparison
  --    (ppComp_blockwise |xs| cpX cpY)  - so   
  --                            X-power-products > Y-power-products,
  --    gs = GroebnerBasis ([p1..pm]++hs)  in  C.
  -- 3. rels' =  [g | g <- gs  and  g  does not depend on X]
  -- 4. rels  =  resl'  embedded in  k[Y].                 


algRelsPols hs ys oY fs = 
  let
    [hsN, fsN]         = map genericLength [hs, fs]
    (Pol _ c oX xs dC) = head fs
                            -- the case [] is processed after `msg0'

    [ysN, xsN]         = map genericLength [ys, xs]
    (xZeroes, yZeroes) = (map (const 0) xs, map (const 0) ys)
    (un, cpX, cpY)     = (unity c, ppoComp oX, ppoComp oY)
    xsys               = xs ++ ys  
    cp                 = ppComp_blockwise xsN cpX cpY
    oXY                = (("", xsN+fsN), cp, [])

                                           -- A = k[X] -> C = k[X,Y]
    toC f = let extPP (Vec js) = Vec (js ++ yZeroes)
                (cs, exps)     = unzip $ polMons f
                exps'          = map extPP exps
            in  Pol (zip cs exps') un oXY xsys dC

                                           -- C = k[X,Y] -> B = k[Y]
    fromCtoB f = let dropX (Vec js) = Vec $ genericDrop xsN js
                     (cs, exps)     = unzip $ polMons f
                     exps'          = map dropX exps
                 in  Pol (zip cs exps') un oY ys dC

    (hCs, fCs) = (map toC hs, map toC fs)

                                    -- forming  pi = fi - yi  in C
    rows  = scalarMt ys 1 0   
    y_pps = map (xZeroes ++) rows
    yMons = [(un, Vec js) | js <- y_pps]   -- yi  as monomial
    ps    = 
       zipWith (\ fi mon -> fi-(Pol [mon] un oXY xsys dC)) fCs yMons

    dependsNotOnX f = 
                   all (all (== 0)) $
                   map (genericTake xsN . vecRepr . snd) $ polMons f

    (gs, _) = gBasis (ps ++ hCs)
    ----------------------------------------------------------------
    msg  = msg0 . msgDom . msgHd hs "hs" . msgHd fs "fs" 
    msg0 = ("algRelsPols hs ys ordY fs, "++) . 
           ("\n(length hs, length fs) =  "++) . shows (hsN,fsN) .
           ("\nys                     =  "++) . shows ys
    msgDom = case  hs ++ fs  
             of 
             []  -> id
             g:_ -> ("\nh_i, f_i <-  "++) . showsDomOf g

    msgHd gs str = case gs of  
                       []  -> id
                       g:_ -> (("\nhead "++str++"  = ")++) . shows g
    ----------------------------------------------------------------
  in
  case (fs, ysN == fsN, any isZero fs) 
                       -- (any isZero fs, fsN == (genericLength ys)) 
  of                                            
  ([],   _  , _   ) -> error$ msg "\n\nfs = []\n"
  (_ , False, _   ) -> error$ msg "\n\nlength ys /= length fs \n"
  (_ , _    , True) -> error$ msg "\n\nZeroes among  fs\n" 
  _                 -> map fromCtoB $ filter dependsNotOnX gs
