--* From postmaster%watson.vnet.ibm.com@yktvmv.watson.ibm.com  Wed Jan 26 09:37:21 1994
--* Received: from yktvmv.watson.ibm.com by leonardo.watson.ibm.com (AIX 3.2/UCB 5.64/4.03)
--*           id AA26305; Wed, 26 Jan 1994 09:37:21 -0500
--* X-External-Networks: yes
--* Received: from watson.vnet.ibm.com by yktvmv.watson.ibm.com (IBM VM SMTP V2R3)
--*    with BSMTP id 9613; Wed, 26 Jan 94 09:43:28 EST
--* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0"
--*           id <A.BRONSTEI.NOTE.YKTVMV.3435.Jan.26.09:43:26.-0500>
--*           for asbugs@watson; Wed, 26 Jan 94 09:43:27 -0500
--* Received: from bernina.ethz.ch by watson.ibm.com (IBM VM SMTP V2R3) with TCP;
--*    Wed, 26 Jan 94 09:43:24 EST
--* Received: from neptune (actually neptune.ethz.ch) by bernina.ethz.ch
--*           with SMTP inbound; Wed, 26 Jan 1994 15:43:10 +0100
--* From: Manuel Bronstein <bronstei@inf.ethz.ch>
--* Received: from rutishauser.inf.ethz.ch (rutishauser-gw.inf.ethz.ch) by neptune
--*           id AA23448; Wed, 26 Jan 94 15:43:03 +0100
--* Date: Wed, 26 Jan 94 15:43:01 +0100
--* Message-Id: <9401261443.AA03117@rutishauser.inf.ethz.ch>
--* Received: from ru7.inf.ethz.ch.rutishauser by rutishauser.inf.ethz.ch
--*           id AA03117; Wed, 26 Jan 94 15:43:01 +0100
--* To: asbugs@watson.ibm.com
--* Subject: [1] Getting compile-time seg.fault or runtime core dump
--*          [bugore.as][33.2]

--@ Fixed  by:  PAB   Mon Mar 21 16:50:50 EST 1994 
--@ Tested by:  ore0.as 
--@ Summary:    Improved A#'s hashing scheme 

----------------------------- bugore.as --------------------------------
-- This file illustrates 2 different bugs, depending on whether it is
-- compiled as a whole, or whether the 2 subfiles  fullore.as and testore.as
-- are compiled separately.
-- In the first case, I get a runtime bug, in the second a compile time bug.
--
-- Bug #1: compiling and running the complete file produces at runtime:
-- peano [/u/manuel/a#] 52> bugore
-- Illegal instruction (core dumped)
--
-- Bug #2: if fullore.as and testore.as are compiled separately (with
--         the first 3 lines of testore.as uncommented, then fullore.as
--         compiles ok, and I get:
-- peano [/u/manuel/a#] 49> asharp -v -L. testore.as
-- A# version 33.2 for AIX RS/6000 (debug version)
-- (Error) Program fault (segmentation violation).


----------------------------- fullore.as --------------------------------
#include "aslib.as"

-- NonNegativeInteger is not yet in aslib, so we use Integer instead for now
macro NonNegativeInteger == Integer

macro {
    N == NonNegativeInteger;
    Z == Integer;
}

--% RetractableTo
RetractableTo(T:Type):Category == with {
	coerce :      T -> %;
	retract:      % -> T;
	retractIfCan: % -> Partial T;
	default retract(x:%):T == {
				import Partial T;
				not failed?(u := retractIfCan x) => u::T;
				error "retract: value is not retractable"
			}
}

--% CommutativeRing
CommutativeRing:Category == Ring

--% Module
Module(R:CommutativeRing):Category == Ring with {
	*: (R, %) -> %;
	*: (%, R) -> %
}

--% Algebra
Algebra(R:CommutativeRing):Category == Module R with {
	coerce: R -> %;
	default coerce(r:R):% == r * 1
}

--% Applicable
Applicable(A:Type, B:Type):Category == with { apply: (%, A) -> B }


--% OreCoefficientFieldCategory
+++ Author: Manuel Bronstein
+++ Date Created: 19 October 1993
+++ Date Last Updated: 1 November 1993
+++ References:
+++ Description:
+++ OreCoefficientFieldCategory is the category of coefficient fields for
+++ skew polynomial rings. Such fields must export an endomorphism \sigma
+++ as well as a \sigma-derivation \delta.
+++ We make the implementation restriction to require \sigma to be an
+++ automorphism of the field.
OreCoefficientFieldCategory: Category == Field with {
	sigma: (%, ..:Z == 1) -> %;
		++ sigma(x,n) returns sigma applied n times to x, n defaults to 1;
		++ if \spad{n < 0} then it returns the inverse of sigma applied
		++ \spad{-n} times to x.
	delta: % -> %;
		++ delta is a sigma-derivation on the field, i.e. it satisfies
		++    \spad{\delta(a + b) = \delta a + \delta b}
		++    \spad{\delta(a b) = \sigma(a) \delta b + b \delta a}.
	constant?: % -> Boolean;
		++ constant?(x) returns \spad{true} if x is constant, i.e.
		++ if \spad{\sigma(x) = x} and \spad{\delta x = 0}.
	default constant?(x:%):Boolean == zero? delta x and x = sigma x
}

--% OreCoefficientField
+++ Author: Manuel Bronstein
+++ Date Created: 19 October 1993
+++ Date Last Updated: 1 November 1993
+++ References:
+++ Description:
+++ OreCoefficientField(R) makes an OreCoefficientFieldCategory of R
+++ using the given maps for \sigma and \delta.
OreCoefficientField(R:Field, sig:(R, Z) -> R, del:R -> R):
	OreCoefficientFieldCategory with {
		coerce: R -> %;
		coerce: % -> R
} == R add {
	macro Rep == R;
	import Rep;

	sigma(x:%, n:Z):%        == per sig(rep x, n);
	delta(x:%):%             == per del rep x;
	coerce(r:R):%            == per r;
	coerce(x:%):R            == rep x
}

--% MonogenicLinearOperator
+++ Author: Stephen Watt
+++ Date Created: 1986
+++ Date Last Updated: 26 January 1994
+++ References:
+++ Description:
+++ MonogenicLinearOperator(R) is a generic linear operator type.
MonogenicLinearOperator(R:Ring):Category == Join(Ring, RetractableTo R) with {
    degree:             % -> N;
        ++ degree(sum_{i=0}^n monomial(a_i,i)) = n if a_n \ne 0.
    minimumDegree:      % -> N;
        ++ minimumDegree(sum_{i=0}^n monomial(a_i,i)) is the smallest i
        ++ such that a_i \ne 0.
    leadingCoefficient: % -> R;
        ++ leadingCoefficient(sum_{i=0}^n monomial(a_i,i)) = a_n if a_n \ne 0.
    reductum:           % -> %;
        ++ reductum(sum_{i=0}^n monomial(a_i,i)) =
        ++          sum_{i=0}^{n-1} monomial(a_i,i)a_n if a_n \ne 0.
    coefficient:        (%, N) -> R;
        ++ coefficient(sum_{i=0}^n monomial(a_i,i), m) = a_m.
    monomial:           (R, N) -> %;
        ++ monomial(c, n) = c monomial(1,1)^n, where monomial(1,1) is the
        ++ generating operator.
    coefficients:       % -> List R;
        ++ coefficients(sum_{i=0}^n monomial(a_i,i)) = [a_n,a_{n-1},\dots,a_0]
        ++ with all the zero coefficients removed.
	default coefficients(l:%):List(R) == {
				import List(R);
				ans:List(R) := empty();
				while l ~= 0 repeat {
					ans := cons(leadingCoefficient l, ans);
					l   := reductum l;
				}
				ans
			}
}

--% UnivariateSkewPolynomialCategory
+++ Author: Manuel Bronstein
+++ Date Created: 19 October 1993
+++ Date Last Updated: 26 January 1994
+++ Description:
+++   This is the category of univariate skew polynomials over an Ore
+++   coefficient field.
+++   The multiplication is given by \spad{x a = \sigma(a) x + \delta a}.
+++   This category is an evolution of the types
+++     MonogenicLinearOperator, OppositeMonogenicLinearOperator, and
+++     NonCommutativeOperatorDivision
+++   developped in Axiom by Jean Della Dora and Stephen M. Watt.
UnivariateSkewPolynomialCategory(F:OreCoefficientFieldCategory):
	Category == Join(Algebra F,MonogenicLinearOperator F,Applicable(%,%)) with {
		apply: (%, F, F) -> F;
			++ apply(p, c, m) returns \spad{p(m)} where the action is
			++ given by \spad{x m = c sigma(m) + delta(m)}.
		exactQuotient: (%, F) -> Partial %;
			++ exactQuotient(l, a) returns the exact quotient of l by a.
		content: % -> F;
			++ content(l) returns the gcd of all the coefficients of l.
		primitivePart: % -> %;
			++ primitivePart(l) returns l0 such that \spad{l = a * l0}
			++ for some a in F, and \spad{content(l0) = 1}.
		leftDivide:   (%, %) -> (%, %);
			++ leftDivide(a,b) returns the pair \spad{(q,r)} such that
			++ \spad{a = b*q + r} and the degree of \spad{r} is
		leftQuotient:  (%, %) -> %;
			++ leftQuotient(a,b) computes the pair \spad{(q,r)} such that
			++ \spad{a = b*q + r} and the degree of \spad{r} is
			++ less than the degree of \spad{b}.
			++ The value \spad{q} is returned.
		leftRemainder:  (%, %) -> %;
			++ leftRemainder(a,b) computes the pair \spad{(q,r)} such that
			++ \spad{a = b*q + r} and the degree of \spad{r} is
			++ less than the degree of \spad{b}.
			++ The value \spad{r} is returned.
		leftExactQuotient:(%, %) -> Partial %;
			++ leftExactQuotient(a,b) computes the value \spad{q}, if it exists,
			++  such that \spad{a = b*q}.
		leftGcd:   (%, %) -> %;
			++ leftGcd(a,b) computes the value \spad{g} of highest degree
			++ such that
			++    \spad{a = g*aa}
			++    \spad{b = g*bb}
			++ for some values \spad{aa} and \spad{bb}.
			++ The value \spad{g} is computed using left-division.
		leftExtendedGcd:   (%, %) -> Record(coef1:%, coef2:%, generator:%);
			++ leftExtendedGcd(a,b) returns \spad{[c,d]} such that
			++ \spad{g = a * c + b * d = leftGcd(a, b)}.
		rightLcm:   (%, %) -> %;
			++ rightLcm(a,b) computes the value \spad{m} of lowest degree
			++ such that \spad{m = a*aa = b*bb} for some values
			++ \spad{aa} and \spad{bb}.  The value \spad{m} is
			++ computed using left-division.
		rightDivide:   (%, %) -> (%, %);
			++ rightDivide(a,b) returns the pair \spad{(q,r)} such that
			++ \spad{a = q*b + r} and the degree of \spad{r} is
			++ less than the degree of \spad{b}.
			++ This process is called ``right division''.
		rightQuotient:  (%, %) -> %;
			++ rightQuotient(a,b) computes the pair \spad{(q,r)} such that
			++ \spad{a = q*b + r} and the degree of \spad{r} is
			++ less than the degree of \spad{b}.
			++ The value \spad{q} is returned.
		rightRemainder:  (%, %) -> %;
			++ rightRemainder(a,b) computes the pair \spad{(q,r)} such that
			++ \spad{a = q*b + r} and the degree of \spad{r} is
			++ less than the degree of \spad{b}.
			++ The value \spad{r} is returned.
		rightExactQuotient:(%, %) -> Partial %;
			++ rightExactQuotient(a,b) computes the value \spad{q}, if it exists
			++ such that \spad{a = q*b}.
		rightGcd:   (%, %) -> %;
			++ rightGcd(a,b) computes the value \spad{g} of highest degree
			++ such that
			++    \spad{a = aa*g}
			++    \spad{b = bb*g}
			++ for some values \spad{aa} and \spad{bb}.
			++ The value \spad{g} is computed using right-division.
		rightExtendedGcd:   (%, %) -> Record(coef1:%, coef2:%, generator:%);
			++ rightExtendedGcd(a,b) returns \spad{[c,d]} such that
			++ \spad{g = c * a + d * b = rightGcd(a, b)}.
		leftLcm:   (%, %) -> %;
			++ leftLcm(a,b) computes the value \spad{m} of lowest degree
			++ such that \spad{m = aa*a = bb*b} for some values
			++ \spad{aa} and \spad{bb}.  The value \spad{m} is
			++ computed using right-division.
		default {
			rightLcm(a:%, b:%):%          == nclcm(a, b, leftEEA);
			leftLcm(a:%, b:%):%           == nclcm(a, b, rightEEA);
			rightGcd(a:%, b:%):%          == ncgcd(a, b, rightRemainder);
			leftGcd(a:%, b:%):%           == ncgcd(a, b, leftRemainder);
			leftQuotient(a:%, b:%):%      == { (q, r) := leftDivide(a,b); q }
			leftRemainder(a:%, b:%):%     == { (q, r) := leftDivide(a,b); r }
			rightQuotient(a:%, b:%):%     == { (q, r) := rightDivide(a,b); q }
			rightRemainder(a:%, b:%):%    == { (q, r) := rightDivide(a,b); r }
			coerce(x:F):%                 == { import N; monomial(x, 0) }
			exactQuot(q:%, r:%):Partial % == { zero? r => [q]; failed }
			leftExactQuotient(a:%, b:%):Partial % == exactQuot leftDivide(a, b);
			rightExactQuotient(a:%,b:%):Partial % == exactQuot rightDivide(a,b);

			leftExtendedGcd(a:%, b:%):Record(coef1:%, coef2:%, generator:%) ==
				extended(a, b, leftEEA);

			rightExtendedGcd(a:%, b:%):Record(coef1:%, coef2:%, generator:%) ==
				extended(a, b, rightEEA);

			retractIfCan(x:%):Partial F == {
				import N;
				zero? degree x => [leadingCoefficient x];
				failed
			}

			(x:%) * (y:%):% == {
				zero? y => 0;
				z:% := 0;
				while x ~= 0 repeat {
					z := z + termPoly(leadingCoefficient x, degree x, y);
					x := reductum x
				}
				z
			}

			(a:F) * (y:%):% == {
				z:% := 0;
				while y ~= 0 repeat {
					z := z + monomial(a * leadingCoefficient y, degree y);
					y := reductum y
				}
				z
			}

			termPoly(a:F, n:N, y:%):% == {
				zero? y => 0;
				(n1 := n - 1) < 0 => a * y;
				z:% := 0;
				while y ~= 0 repeat {
					m := degree y;
					b := leadingCoefficient y;
					z := z + termPoly(a, n1, monomial(sigma b, m + 1))
							+ termPoly(a, n1, monomial(delta b, m));
					y := reductum y;
				}
				z
			}

			apply(x:%, y:%):% == {
				import N; import Segment N;
				ans:% := 0;
				zero? x => ans;
				yn:%  := 1;
				for i in 0..degree x repeat {
					ans := ans + coefficient(x, i) * yn;
					yn  := y * yn;
				}
				ans
			}

			-- this is meant for when the parameter won't be a field any more.
			exactQuotient(l:%, a:F):Partial(%) == {
				import F;
				ans:% := 0;
				while l ~= 0 repeat {
					(q, r) := divide(leadingCoefficient l, a);
					r ~= 0 => return failed;
					ans := ans + monomial(q, degree l);
					l   := reductum l;
				}
				[ans]
			}

			-- this is meant for when the parameter won't be a field any more.
			content(l:%):F == {
				import List F;
				empty?(cf := coefficients l) => 0;
				g := first cf;
				while not empty?(cf := rest cf) repeat g := gcd(g, first cf);
				g;
			}

			-- this is meant for when the parameter won't be a field any more.
			primitivePart(l:%):% == {
				import Partial %;
				retract exactQuotient(l, content l)
			}

			apply(p:%, c:F, x:F):F == {
				import N; import Segment N;
				w:F  := 0;
				xn:F := x;
				for i in 0..degree p repeat {
					w  := w + coefficient(p, i) * xn;
					xn := c * sigma xn + delta xn;
				}
				w
			}

			-- leftDivide(a, b) returns [q, r] such that a = b q + r
			leftDivide(a:%, b:%):(%, %) == {
				import N; import F;
				-- must comment out the error check, due to a compiler bug
				-- in version 33.2,  MB 1/94
				--zero? b => error "leftDivide: division by 0";
				n := degree(a) - (m := degree b);
				zero? a or n < 0 => (0, a);
				q0 := monomial(
					sigma(leadingCoefficient a / leadingCoefficient b, -m), n);
				(q, r) := leftDivide(a - b * q0, b);
				(q + q0, r)
			}

			-- returns [g = leftGcd(a, b), c, d, l = rightLcm(a, b)]
			-- such that g := a c + b d
			leftEEA(a:%, b:%):(%, %, %, %) == {
				a0 := a;
				u0:% := v:% := 1;
				v0:% := u:% := 0;
				while b ~= 0 repeat {
					(q, r) := leftDivide(a, b);
					(a, b) := (b, r);
					(u0, u):= (u, u0 - u * q);
					(v0, v):= (v, v0 - v * q)
				}
				(a, u0, v0, a0 * u)
			}

			-- rightDivide(a, b) returns [q, r] such that a = q b + r
			rightDivide(a:%, b:%):(%, %) == {
				import N; import F;
				-- must comment out the error check, due to a compiler bug
				-- in version 33.2,  MB 1/94
				--zero? b => error "rightDivide: division by 0";
				n := degree(a) - (m := degree b);
				zero? a or n < 0 => (0, a);
				q0 := monomial(leadingCoefficient a /
											sigma(leadingCoefficient b, n), n);
				(q, r) := rightDivide(a - q0 * b, b);
				(q + q0, r)
			}

			ncgcd(a:%, b:%, ncrem:(%, %) -> %):% == {
				import N;
				zero? a => b;
				zero? b => a;
				degree a < degree b => ncgcd(b, a, ncrem);
				while b ~= 0 repeat (a, b) := (b, ncrem(a, b));
				a
			}

			extended(a:%, b:%, eea:(%,%) -> (%,%,%,%)):
				Record(coef1:%, coef2:%, generator:%) == {
				import N;
				zero? a => [0, 1, b];
				zero? b => [1, 0, a];
				degree a < degree b => {(g, c1, c2, l) := eea(b,a); [c2, c1, g]}
				(g, c1, c2, l) := eea(a, b);
				[c1, c2, g]
			}

			nclcm(a:%, b:%, eea:(%,%) -> (%,%,%,%)):% == {
				import N;
				zero? a or zero? b => 0;
				degree a < degree b => nclcm(b, a, eea);
				(g, c1, c2, l) := eea(b, a);
				l
			}

			-- returns [g = rightGcd(a, b), c, d, l = leftLcm(a, b)]
			-- such that g := a c + b d
			rightEEA(a:%, b:%):(%, %, %, %) == {
				a0 := a;
				u0:% := v:% := 1;
				v0:% := u:% := 0;
				while b ~= 0 repeat {
					(q, r) := rightDivide(a, b);
					(a, b) := (b, r);
					(u0, u):= (u, u0 - q * u);
					(v0, v):= (v, v0 - q * v)
				}
				(a, u0, v0, u * a0)
			}
		}   -- end from default definitions
}

--% SparseUnivariateSkewPolynomial
+++ Author: Manuel Bronstein
+++ Date Created: 2 November 1993
+++ Date Last Updated: 26 January 1994
+++ Description:
+++ SparseUnivariateSkewPolynomial is the domain of sparse univariate skew
+++ polynomials over an Ore coefficient field.
+++ The multiplication is given by \spad{x a = \sigma(a) x + \delta a}.
SparseUnivariateSkewPolynomial(F:OreCoefficientFieldCategory):
 UnivariateSkewPolynomialCategory F with {
	apply: (OutPort, %, String) -> OutPort;
		++ apply(p, q, x) prints q on p using x for the anonymous variable.
  	} == add {
		Term: BasicType with {
			term: (F, N) -> %;    -- term(c, n) creates the monomial c x^n
			coef: % -> F;         -- coef(c x^n) returns c
			expt: % -> N;         -- expt(c x^n) returns n
		} == add {
			macro Rep == Record(coef:F, expt:N);
			import Rep;

			term(c:F, n:N):% == per [c, n];
			coef(t:%):F      == rep(t).coef;
			expt(t:%):N      == rep(t).expt;
			sample:%         == term(1, 1);
			apply(p:OutPort, x:%):OutPort == p;
            (u:%) = (v:%):Boolean == coef u = coef v and expt u = expt v;
		}
		macro Rep == List Term;   -- sorted, highest term first
		import Term;
		import Rep;
		import F;
		import N;

		0:%                       == per empty();
		1:%                       == per list term(1, 0);
		degree(p:%):N             == {zero? p => 0; expt first rep p}
		minimumDegree(p:%):N      == {zero? p => 0; expt last rep p}
		leadingCoefficient(p:%):F == {zero? p => 0; coef first rep p}
    	reductum(p:%):%           == {zero? p => 0; per rest rep p}
		monomial(c:F, n:N):%      == {zero? c => 0; per list term(c, n)}
  		coerce(c:F):%             == monomial(c, 0);
		(x:%) = (y:%):Boolean     == rep x = rep y;
		-(p:%):%                  == per map(cterm(-1), rep p);
  		(c:F) * (p:%):%           == per map(cterm c, rep p);
  		(p:%) * (c:F):%           == p * (c::%);
		cterm(c:F):(Term->Term)   == (t:Term):Term +-> term(c * coef t,expt t);
		apply(p:OutPort, x:%):OutPort == {import String; apply(p, x, " ?")}

		(x:%) + (y:%):% == {
			zero? x => y; zero? y => x;
			nx := degree x; ny := degree y;
			nx > ny => per cons(first rep x, rep(reductum x + y));
			nx < ny => per cons(first rep y, rep(x + reductum y));
			z := reductum x + reductum y;
			zero?(c := leadingCoefficient x + leadingCoefficient y) => z;
			per cons(term(c, nx), rep z);
		}

		applyTerm(p:OutPort, c:F, n:N, v:String):OutPort == {
			p := p(c);
			if n > 0 then {p := p v; if n > 1 then p := p("^")(n);}
			p;
		}

		apply(p:OutPort, x:%, v:String):OutPort == {
			zero? x => p(0$F);
			while x ~= 0 repeat {
				p := applyTerm(p, leadingCoefficient x, degree x, v);
				x := reductum x;
				if (x ~= 0) then p := p(" + ");
			}
			p;
		}

		retractIfCan(p:%):Partial F == {
			zero? p or zero? degree p => [leadingCoefficient p];
			failed
		}

		coefficient(p:%, n:N):F == {
			for t in rep p repeat { n = expt t => return coef t }
			0
		}

}

------------------------------ testore.as ---------------------------------
--#include "aslib.as"

--#library ore       "fullore.aso"
--import ore

-- NonNegativeInteger is not yet in aslib, so we use Integer instead for now
macro NonNegativeInteger == Integer

macro {
	N == NonNegativeInteger;
	Z == Integer;
	Q == Ratio Z;
	ORE == OreCoefficientField(Q, ident, zero);
	ORESUP == SparseUnivariateSkewPolynomial ORE;
}

ident(x:Q, n:Z):Q == x;
zero(x:Q):Q == 0;

main():Z == {
	import N; import ORE; import ORESUP;
	a:ORE := 1;
	print("a = ")(a)();
	p:ORESUP := monomial(a, 2) + monomial(a, 0);
	print("x^2 + 1 = ")(p)();
	0
}

main()

 
--+ Runtime hashing now in the system (actually coming soon). The second bug
--+ regarding splitting the file is still open, and I have resubmitted it.
--+ 
--+ Pete.
