--* From BMT%WATSON.vnet.ibm.com@yktvmh.watson.ibm.com  Wed Oct 27 18:24:49 1993
--* Received: from yktvmh.watson.ibm.com by radical.watson.ibm.com (AIX 3.2/UCB 5.64/900524)
--*           id AA20581; Wed, 27 Oct 1993 18:24:49 -0400
--* Received: from watson.vnet.ibm.com by yktvmh.watson.ibm.com (IBM VM SMTP V2R3)
--*    with BSMTP id 2558; Wed, 27 Oct 93 18:31:46 EDT
--* Received: from YKTVMH by watson.vnet.ibm.com with "VAGENT.V1.0"
--*           id <A.BMT.NOTE.VAGENT2.4708.Oct.27.18:31:45.-0400>
--*           for asbugs@watson; Wed, 27 Oct 93 18:31:46 -0400
--* Received: from YKTVMH by watson.vnet.ibm.com with "VAGENT.V1.0"
--*           id 4705; Wed, 27 Oct 1993 18:31:44 EDT
--* Received: from cyst.watson.ibm.com by yktvmh.watson.ibm.com (IBM VM SMTP V2R3)
--*    with TCP; Wed, 27 Oct 93 18:31:43 EDT
--* Received: from spadserv.watson.ibm.com by cyst.watson.ibm.com (AIX 3.2/UCB 5.64/900528)
--*   id AA84728; Wed, 27 Oct 1993 18:31:19 -0400
--* Received: by spadserv.watson.ibm.com (AIX 3.2/UCB 5.64/900524)
--*           id AA15273; Wed, 27 Oct 1993 18:33:24 -0400
--* Date: Wed, 27 Oct 1993 18:33:24 -0400
--* From: bmt@spadserv.watson.ibm.com
--* X-External-Networks: yes
--* Message-Id: <9310272233.AA15273@spadserv.watson.ibm.com>
--* To: asbugs@watson.ibm.com
--* Subject: compiler complains about missing gcd, quo, rem but these are defaulted by Field [algext.as][32.1 (current)]

--@ Fixed  by:  SSD   Tue Nov 8 09:33:18 EST 1994 
--@ Tested by:  none 
--@ Summary:    Bug no longer appears in r1-0-2. 

-- simple algebraic extension fields

#include "aslib.as"

macro {
	SI    == SingleInteger;
	B     == Boolean;
}

SimpleAlgebraicExtensionField(F: Field, n:SingleInteger, minp:PrimitiveArray(F)):  with {
        Field;
        *: (F, %) -> %;
--        norm: % -> F;
        trace: % -> F;
        coordinates: % -> PrimitiveArray(F);
        represents: PrimitiveArray(F) -> %;
        generator: %;
        definingPolynomial: PrimitiveArray(F);
        map:      (F->F, %) -> %;
	map:	  ((F,F)->F, %, %) -> %;
}
== add {
	Rep ==> PrimitiveArray F;
	import from Rep;
	import from SI;
	import from F;
        -- minp.(n+1) ~= 1 => error "defining polynomial must be monic";
	_^ (x:%, i:Integer):% == x;
        coordinates(x:%):Rep == rep x;
        represents(v:Rep):% == per v;
        0:% == per new(n, 0);
        1:% == (v:=new(n,0); v.1 := 1; per v);
        generator:% == (v:=new(n,0); v.2 := 1; per v);
        definingPolynomial:Rep == minp;
	map(f:F-> F, v: %):% == {
		 vv:Rep := new(n);
		 for i in 1..n repeat
		 	vv.i := f((rep v).i);
		 per vv
	}
	map(f:(F,F) -> F, v1:%, v2:%) :% == {
		 vv:Rep := new(n);
		 for i in 1..n repeat
		 	vv.i := f(rep(v1).i, rep(v2).i);
		 per vv
	}
	-(v:%):%          == map(-,v);
	(v1:%) + (v2:%):% == map(+,v1,v2);
	(v1:%) - (v2:%):% == map(-,v1,v2);
	(c: F) * (x: %):% == map((cc:F):F+-> c*cc, x);

        fmecg!(v1:Rep,e:SI,r:F,v2:Rep):Rep == {      -- p1 - r * X**e * p2
           for i in 1..n repeat
              (ii:= e+i; v1.ii:=v1.ii - r * v2.i);
           v1
        }

        pmult(v1:Rep, v2:Rep):Rep == { -- polynomial style product
          prodv := new(2*n-1,0);
          for i in 1..n repeat {
             if not zero?(c := v1.i) then
                for j in 1..n for k in i.. repeat
                    prodv.k := prodv.k + c * v2.j
          }
          prodv
        }

        (x: %) * (y: %): % == {
          prodv := pmult(rep x, rep y);
          for i in (2*n-1)..(n+1) by -1 repeat
             if not zero?(c:=prodv.i) then {
                 prodv.i := 0;
                 for j in 1..n for k in i-n+1 .. repeat
                    prodv.k := prodv.k - c * minp.j
             }
          ans := new(n);
          for i in 1..n repeat ans.i := prodv.i;
          per ans
        }

        inv(x:%): % == {
           zero? x => error "inverse of zero";
           x; -- temporary
           --per(extendedEuclidean(rep x, minp, 1).coef1);
        }

        trace(x:%): F == {
           xn := x;
           ans:F := rep(xn).1;
           for i in 2..n repeat {
              xn := generator * xn; -- this could be special cased
              ans := rep(xn).i + ans;
           }
           ans;
           }
--       norm(x:%):F          == resultant(definingPolynomial(), lift x)

	(v1:%) = (v2:%):B == {
		 for i: SI in 1..n repeat
		 	rep(v1).i ~= rep(v2).i => return false;
		 true
	}

	apply(p: OutPort, v: %): OutPort == {
		  p "[";
                  p(rep(v).1);
		  for i in 2..n repeat p(", ")(rep(v).i);
		  p "]";
	}
}