--Copyright (C) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--Copyright (C) 2007-2010, Gabriel Dos Reis.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

)abbrev category IDPC IndexedDirectProductCategory
++ Author: James Davenport, Gabriel Dos Reis
++ Date Created:
++ Date Last Updated: June 28, 2010
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This category represents the direct product of some set with
++ respect to an ordered indexing set.

IndexedDirectProductCategory(A:BasicType,S:OrderedType): Category ==
  BasicType with
    if A has SetCategory and S has SetCategory then SetCategory
    map:           (A -> A, %) -> %
       ++ map(f,z) returns the new element created by applying the
       ++ function f to each component of the direct product element z.
    monomial:         (A, S) -> %
       ++ monomial(a,s) constructs a direct product element with the s
       ++ component set to \spad{a}
    leadingCoefficient:   % -> A
       ++ leadingCoefficient(z) returns the coefficient of the leading
       ++ (with respect to the ordering on the indexing set)
       ++ monomial of z.
       ++ Error: if z has no support.
    leadingSupport:   % -> S
       ++ leadingSupport(z) returns the index of leading
       ++ (with respect to the ordering on the indexing set) monomial of z.
       ++ Error: if z has no support.
    reductum:      % -> %
       ++ reductum(z) returns a new element created by removing the
       ++ leading coefficient/support pair from the element z.
       ++ Error: if z has no support.
    terms: % -> List Pair(S,A)
      ++ \spad{terms x} returns the list of terms in \spad{x}.
      ++ Each term is a pair of a support (the first component)
      ++ and the corresponding value (the second component).

)abbrev domain IDPO IndexedDirectProductObject
++ Author: James Davenport, Gabriel Dos Reis
++ Date Created:
++ Date Last Updated: June 28, 2010
++ Description:
++   Indexed direct products of objects over a set \spad{A}
++   of generators indexed by an ordered set S. All items have finite support.
IndexedDirectProductObject(A,S): Public == Private where
  A: BasicType
  S: OrderedType
  Public == IndexedDirectProductCategory(A,S)
  Private == add
    Term == Pair(S,A)
    Rep == List Term
    -- Return the index of a term
    termIndex(t: Term): S == first t
    -- Return the value of a term
    termValue(t: Term): A == second t

    x = y ==
      x' := rep x
      y' := rep y
      while not null x' and not null y' repeat
        termIndex first x' ~= termIndex first y' => return false
        termValue first x' ~= termValue first y' => return false
        x' := rest x'
        y' := rest y'
      null x' and null y'

    if A has CoercibleTo OutputForm and S has CoercibleTo OutputForm then
      coerce(x:%):OutputForm ==
         bracket [rarrow(termIndex(t)::OutputForm, termValue(t)::OutputForm)
                   for t in rep x]

    -- sample():% == [[sample()$S,sample()$A]$Term]$Rep

    monomial(r,s) == per [[s,r]]
    map(f,x) == per [[termIndex tm,f termValue tm] for tm in rep x]

    reductum x ==
       per rest rep x
    leadingCoefficient x  ==
       null rep x =>
         error "Can't take leadingCoefficient of empty product element"
       termValue first rep x
    leadingSupport x  ==
       null rep x =>
         error "Can't take leadingCoefficient of empty product element"
       termIndex first rep x

    terms x == rep x

)abbrev domain IDPAM IndexedDirectProductAbelianMonoid
++ Indexed direct products of abelian monoids over an abelian monoid \spad{A} of
++ generators indexed by the ordered set S. All items have finite support.
++ Only non-zero terms are stored.
IndexedDirectProductAbelianMonoid(A:AbelianMonoid,S:OrderedType):
    Join(AbelianMonoid,IndexedDirectProductCategory(A,S))
 ==  IndexedDirectProductObject(A,S) add
    --representations
       Term ==  Pair(S,A)
       import Term
       import List Term
       termIndex(t: Term): S == first t
       termValue(t: Term): A == second t

       r: A
       n: NonNegativeInteger
       f: A -> A
       s: S
       0  == nil$List(Term) pretend %
       zero? x ==  null terms x

       import %tail: List Term -> List Term from Foreign Builtin

       qsetrest!: (List Term, List Term) -> List Term
       qsetrest!(l, e) ==
         %store(%tail l,e)$Foreign(Builtin)

       -- PERFORMANCE CRITICAL; Should build list up
       --  by merging 2 sorted lists.   Doing this will
       -- avoid the recursive calls (very useful if there is a
       -- large number of vars in a polynomial.
       x + y ==
         x' := terms x
         y' := terms y
         null x' => y
         null y' => x
         endcell: List Term := nil
         res: List Term := nil
         while not empty? x' and not empty? y' repeat 
           newcell: List Term := nil
           if termIndex x'.first = termIndex y'.first then
             r := termValue x'.first + termValue y'.first
             if not zero? r then 
               newcell := cons([termIndex x'.first, r], empty())
             x' := rest x'
             y' := rest y'
           else if termIndex x'.first > termIndex y'.first then
             newcell := cons(x'.first, empty())
             x' := rest x'
           else
             newcell := cons(y'.first, empty())
             y' := rest y'
           if not empty? newcell then 
             if not empty? endcell then
               qsetrest!(endcell, newcell)
               endcell := newcell
             else
               res     := newcell;
               endcell := res
         end := 
           empty? x' => y'
           x'
         if empty? res then res := end
         else qsetrest!(endcell, end)
         res pretend %

       n * x  ==
         n = 0 => 0
         n = 1 => x
         [[termIndex u,a] for u in terms x
            | not zero?(a:=n * termValue u)] pretend %

       monomial(r,s) ==
         zero? r => 0
         [[s,r]] pretend %

       map(f,x) ==
         [[termIndex tm,a] for tm in terms x
            | not zero?(a:=f termValue tm)] pretend %

       reductum x ==
         null terms x => 0
         rest(terms x) pretend %

       leadingCoefficient x  ==
         null terms x => 0
         termValue terms(x).first

       pair2Term(t: Pair(A,S)): Term ==
         [second t, first t]

)abbrev domain IDPOAM IndexedDirectProductOrderedAbelianMonoid
++ Indexed direct products of ordered abelian monoids \spad{A} of
++ generators indexed by the ordered set S.
++ The inherited order is lexicographical.
++ All items have finite support: only non-zero terms are stored.
IndexedDirectProductOrderedAbelianMonoid(A:OrderedAbelianMonoid,S:OrderedType):
    Join(OrderedAbelianMonoid,IndexedDirectProductCategory(A,S))
 ==  IndexedDirectProductAbelianMonoid(A,S) add
    --representations
       Term:=  Record(k:S,c:A)
       Rep:=  List Term
       x,y: %
       x<y ==
         empty? y => false
         empty? x => true   -- note careful order of these two lines
         y.first.k > x.first.k => true
         y.first.k < x.first.k => false
         y.first.c > x.first.c => true
         y.first.c < x.first.c => false
         x.rest < y.rest

)abbrev domain IDPOAMS IndexedDirectProductOrderedAbelianMonoidSup
++ Indexed direct products of ordered abelian monoid sups \spad{A},
++ generators indexed by the ordered set S.
++ All items have finite support: only non-zero terms are stored.
IndexedDirectProductOrderedAbelianMonoidSup(A:OrderedAbelianMonoidSup,S:OrderedSet):
    Join(OrderedAbelianMonoidSup,IndexedDirectProductCategory(A,S))
 ==  IndexedDirectProductOrderedAbelianMonoid(A,S) add
    --representations
       Term:=  Record(k:S,c:A)
       Rep:=  List Term
       x,y: %
       r: A
       s: S

       subtractIfCan(x,y) ==
         empty? y => x
         empty? x => "failed"
         x.first.k < y.first.k => "failed"
         x.first.k > y.first.k =>
             t:= subtractIfCan(x.rest, y)
             t case "failed" => "failed"
             cons( x.first, t)
         u:=subtractIfCan(x.first.c, y.first.c)
         u case "failed" => "failed"
         zero? u => subtractIfCan(x.rest, y.rest)
         t:= subtractIfCan(x.rest, y.rest)
         t case "failed" => "failed"
         cons([x.first.k,u],t)

       sup(x,y) ==
         empty? y => x
         empty? x => y
         x.first.k < y.first.k => cons(y.first,sup(x,y.rest))
         x.first.k > y.first.k => cons(x.first,sup(x.rest,y))
         u:=sup(x.first.c, y.first.c)
         cons([x.first.k,u],sup(x.rest,y.rest))

)abbrev domain IDPAG IndexedDirectProductAbelianGroup
++ Indexed direct products of abelian groups over an abelian group \spad{A} of
++ generators indexed by the ordered set S.
++ All items have finite support: only non-zero terms are stored.
IndexedDirectProductAbelianGroup(A:AbelianGroup,S:OrderedType):
    Join(AbelianGroup,IndexedDirectProductCategory(A,S))
 ==  IndexedDirectProductAbelianMonoid(A,S) add
    --representations
       Term == Pair(S,A)
       termIndex(t: Term): S == first t
       termValue(t: Term): A == second t

       -x == [[termIndex u,-termValue u] for u in terms x] pretend %
       n:Integer * x:% ==
         n = 0 => 0
         n = 1 => x
         [[termIndex u,a] for u in terms x
            | not zero?(a := n * termValue u)] pretend %

       import %tail: List Term -> List Term from Foreign Builtin

       qsetrest!: (List Term, List Term) -> List Term
       qsetrest!(l, e) ==
         %store(%tail l,e)$Foreign(Builtin)

       x - y ==
         x' := terms x
         y' := terms y
         null x' => -y
         null y' => x
         endcell: List Term := nil
         res: List Term := nil
         while not empty? x' and not empty? y' repeat
           newcell: List Term := nil
           if termIndex x'.first = termIndex y'.first then
             r := termValue x'.first - termValue y'.first
             if not zero? r then
               newcell := cons([termIndex x'.first, r], empty())
             x' := rest x'
             y' := rest y'
           else if termIndex x'.first > termIndex y'.first then
             newcell := cons(x'.first, empty())
             x' := rest x'
           else
             newcell := cons([termIndex y'.first,-termValue y'.first], empty())
             y' := rest y'
           if not empty? newcell then
             if not empty? endcell then
               qsetrest!(endcell, newcell)
               endcell := newcell
             else
               res     := newcell;
               endcell := res
         end := 
           empty? x' => terms(-(y' pretend %))
           x'
         if empty? res then res := end
         else qsetrest!(endcell, end)
         res pretend %

