root/trunk/Library/FortressLibrary.fss @ 2739

(*******************************************************************************
    Copyright 2008 Sun Microsystems, Inc.,
    4150 Network Circle, Santa Clara, California 95054, U.S.A.
    All rights reserved.

    U.S. Government Rights - Commercial software.
    Government users are subject to the Sun Microsystems, Inc. standard
    license agreement and applicable provisions of the FAR and its supplements.

    Use is subject to license terms.

    This distribution may include materials developed by third parties.

    Sun, Sun Microsystems, the Sun logo and Java are trademarks or registered
    trademarks of Sun Microsystems, Inc. in the U.S. and other countries.
 ******************************************************************************)

component FortressLibrary
import NativeArray.{...}
import NatReflect.{...}
import JavaString.{...}
import List.{AnyList}

export FortressLibrary

(************************************************************
* Simple Combinators
************************************************************)

(** Casting *)

cast[\T extends Any\](x:Any):T =
  typecase x of
    T => x
    else => throw CastError
  end

instanceOf[\T extends Any\](x:Any):Boolean =
  typecase x of
    T => true
    else => false
  end

(** Useful functions *)

ignore(_:Any):() = ()

identity[\T extends Any\](x:T):T = x

(* Should we depracate tuple and use identity instead?  Decision: no. *)
tuple[\T\](x:T):T = x

(* Function composition *)
opr COMPOSE[\A,B,C\](f: B->C, g: A->B): A->C = fn (a:A): C => f(g(a))

(** Reflection of static type parameters for overloading purposes.
    Works around shortcomings in the story on parametric overloading
    (all overloadings must have the same parameters with the same
    bounds).  Allows us to overload a function based on the parametric
    type of an output. *)

object __Proxy[\T extends (* U *) Any\]
    (* extends { __Proxy[\U\], Object } where { U extends Any } *)
  getter toString() = "__Proxy object, should not be user-visible."
end

fail[\T\](s:String):T = do
   println("FAIL: " s)
   (* throwError(s) *)
   throw FailCalled(s)
end

(************************************************************
* Control over locality and location
************************************************************)

(* At the moment all Fortress objects are immediately shared by default. *)

shared[\T extends Any\](x:T): T = x

isShared(x:Any): Boolean = true

localize[\T extends Any\](x:T): T = x

(* copy is presently unimplemented.
copy[\T extends Any\](x:T): T =
*)

trait Region extends Equality[\Region\]
    getter toString(): String
    isLocalTo(r: Region): Boolean = false
    opr =(self, other:Region): Boolean = self SEQV other
end

object Global extends Region
    getter toString(): String = "Global"
    isLocalTo(_: Region): Boolean = true
end

region(a:Any): Region = Global

here(): Region = Global

(************************************************************
* Equality and ordering
************************************************************)

opr =(a:Any, b:Any):Boolean = a SEQV b

opr =/=(a:Any, b:Any):Boolean = NOT (a=b)

trait Equality[\Self extends Equality[\Self\]\]
    opr =(self, other:Self): Boolean
end

(** Total ordering *)

object LexicographicPartialReduction
        extends { MonoidReduction[\Comparison\],
                  ReductionWithZeroes[\Comparison, Comparison\] }
    getter toString():String = "Lexicographic reduction on partial order"
    empty(): Comparison = EqualTo
    join(a:Comparison, b:Comparison):Comparison = a LEXICO b
    isZero(_:Unordered): Boolean = true
    isZero(_:Comparison): Boolean = false
end

object LexicographicReduction
        extends { MonoidReduction[\TotalComparison\],
                  ReductionWithZeroes[\TotalComparison, TotalComparison\] }
    getter toString():String = "Lexicographic reduction on total order"
    empty(): TotalComparison = EqualTo
    join(a:TotalComparison, b:TotalComparison):TotalComparison = a LEXICO b
    isLeftZero(_:EqualTo): Boolean = false
    isLeftZero(_:Comparison): Boolean = true
end

opr BIG LEXICO(): BigReduction[\TotalComparison, TotalComparison\] =
  BigReduction[\TotalComparison,TotalComparison\](LexicographicReduction)
  

trait Comparison
        extends { StandardPartialOrder[\Comparison\] }
        comprises { Unordered, TotalComparison }
    getter toString(): String
    opr =(self, other:Comparison): Boolean = false
    opr LEXICO(self, other:Comparison): Comparison = Unordered
    opr INVERSE(self): Comparison
end

(** Unordered is the outcome of a CMP b when a and b are partially
    ordered and no ordering relationship exists between them. **)
object Unordered extends Comparison
    getter toString(): String = "Unordered"
    opr =(self, other:Unordered): Boolean = true
    opr <(self, other:Comparison): Boolean = false
    opr INVERSE(self): Comparison = Unordered
end

trait TotalComparison
        extends { Comparison, StandardTotalOrder[\TotalComparison\] }
        comprises { LessThan, EqualTo, GreaterThan }
    (* We're both a partial order (including Unordered) and a total
       order (TotalComparison alone).  Avoid ambiguity between the
       default definitions of CMP and >=. *)
    opr =(self, other:TotalComparison): Boolean = false
    opr =(self, other:Unordered): Boolean = false
    opr CMP(self, other:Unordered): Comparison = Unordered
    opr <(self, other:Unordered): Boolean = false
    opr >=(self, other:Unordered): Boolean = false
    opr >=(self, other:Comparison): Boolean = NOT (other < self)
    opr LEXICO(self, other:TotalComparison): TotalComparison = self
    opr LEXICO(self, other:()->TotalComparison): TotalComparison = self
    opr INVERSE(self): TotalComparison
end

object LessThan extends TotalComparison
    getter toString(): String = "LessThan"
    opr =(self, other:LessThan): Boolean = true
    opr CMP(self, other:LessThan): TotalComparison = EqualTo
    opr CMP(self, other:TotalComparison): TotalComparison = GreaterThan
    opr <(self, other:LessThan): Boolean = false
    opr <(self, other:TotalComparison): Boolean = true
    opr INVERSE(self): TotalComparison = GreaterThan
end

object GreaterThan extends TotalComparison
    getter toString(): String = "GreaterThan"
    opr =(self, other:GreaterThan): Boolean = true
    opr CMP(self, other:GreaterThan): TotalComparison = EqualTo
    opr CMP(self, other:TotalComparison): TotalComparison = LessThan
    opr <(self, other:TotalComparison): Boolean = false
    opr INVERSE(self): TotalComparison = LessThan
end

object EqualTo extends TotalComparison
    getter toString(): String = "EqualTo"
    opr =(self, other:EqualTo): Boolean = true
    opr CMP(self, other:TotalComparison): TotalComparison = INVERSE other
    opr <(self, other:GreaterThan): Boolean = true
    opr <(self, other:TotalComparison): Boolean = false
    opr LEXICO(self, other:TotalComparison): TotalComparison = other
    opr LEXICO(self, other:()->TotalComparison): TotalComparison = other()
    opr INVERSE(self): TotalComparison = EqualTo
end

(** StandardPartialOrder is partial ordering using <,>,<=,>=,=, and CMP.
    This is primarily for floating-point values.  Minimal complete
    definition: CMP or { <, = }. **)
trait StandardPartialOrder[\Self extends StandardPartialOrder[\Self\]\]
        extends { Equality[\Self\] }
    opr CMP(self, other:Self): Comparison =
        if self < other then LessThan
        elif other > self then GreaterThan
        elif self = other then EqualTo
        else Unordered
        end
    opr <(self, other:Self): Boolean = LessThan = (self CMP other)
    opr >(self, other:Self): Boolean = other < self
    opr =(self, other:Self): Boolean = EqualTo = (self CMP other)
    opr <=(self, other:Self): Boolean = other >= self
    opr >=(self, other:Self): Boolean = (self = other OR: self > other)
end

(** %StandardMin% is a MIN operator; most types that implement %MIN%
    will implement a corresponding total order.  It's a separate type
    to account for the existence of floating point numbers, for which
    NaN counts as a bottom that is less than anything else but doesn't
    actually participate in the standard total ordering.  It is
    otherwise the case that %a MIN b = a% when %a <= b% and that
    %a MIN b = b MIN a%. **)
trait StandardMin[\T extends StandardMin[\T\]\]
    opr MIN(self, other:T): T
end

(** %StandardMax% is a MAX operator; most types that implement %MAX%
    will implement a corresponding total order.  It's a separate type
    to account for the existence of floating point numbers, for which
    NaN counts as a bottom that is less than anything else but doesn't
    actually participate in the standard total ordering.  It is
    otherwise the case that %a MAX b = a% when %a <= b% and that
    %a MAX b = b MAX a%. **)
trait StandardMax[\T extends StandardMax[\T\]\]
    opr MAX(self, other:T): T
end

(** StandardTotalOrder is the usual total order using <,>,<=,>=,=, and
    CMP.  Most values that define a comparison should do so using
    this.  Minimal complete definition: either CMP or < (it's
    advisable to define = in the latter case).  As noted above, %MIN%
    and %MAX% respect the total order and are defined in the obvious
    way. **)
trait StandardTotalOrder[\Self extends StandardTotalOrder[\Self\]\]
        extends { StandardPartialOrder[\Self\], StandardMin[\Self\], StandardMax[\Self\] }
    opr CMP(self, other:Self): TotalComparison =
        if self < other then LessThan
        elif other < self then GreaterThan
        else EqualTo
        end
    opr >=(self, other:Self): Boolean = NOT (self < other)
    opr <=(self, other:Self): Boolean = NOT (other < self)
    opr MIN(self, other:Self): Self = if other < self then other else self end
    opr MAX(self, other:Self): Self = if self < other then other else self end
end

(** Assertion *)
assert(flag:Boolean): () =
    if NOT flag then
        fail("Assertion failed!")
    end

assert(flag: Boolean, failMsg: String): () =
    if NOT flag then
        fail(failMsg)
    end

assert(x:Any, y:Any, failMsg: Any...): () =
    if x =/= y then
        msg = x " =/= " y "; " (BIG || failMsg)
        fail(msg)
    end
    
deny(flag:Boolean): () = assert(NOT flag)

deny(flag: Boolean, failMsg: String): () = assert(NOT flag, failMsg)

deny(x:Any, y:Any, failMsg: Any...): () =
    if x = y then
        msg = x " = " y "; " (BIG || failMsg)
        fail(msg)
    end
    
shouldRaise⟦Ex extends Exception⟧ (expr: ()→()): () = do
    try
        expr()
        throw ForbiddenException
    catch x
        Ex => assert(true) (* This is what we hope for *) 
    forbid Exception
    end
end
        

(************************************************************
* Numeric hierarchy
************************************************************)

(** Additive group making use of +.  Must define + and either unary or binary -. **)
trait AdditiveGroup[\Self extends AdditiveGroup[\Self\]\]
    getter zero(): Self = self - self
    opr +(self, other: Self): Self
    opr -(self, other: Self): Self = self + (-other)
    opr -(self) : Self = self.zero - self
end

(** Place holder for exclusions of MultiplicativeRing **)
trait AnyMultiplicativeRing end

(** Multiplicative ring using TIMES and juxtaposition.
    Define opr TIMES; juxtaposition is defined in terms of TIMES. **)
trait MultiplicativeRing[\Self extends MultiplicativeRing[\Self\]\]
        extends { AdditiveGroup[\Self\], AnyMultiplicativeRing }
    getter one(): Self
    opr TIMES(self, other:Self): Self
    opr juxtaposition(self, other:Self): Self = self TIMES other
    (** Exponentiation need only deal with natural exponents. **)
    opr ^(self, other:AnyIntegral): Self
end

(** The use of StandardTotalOrder here is a lie; it will have to suffice
    until we can deal gracefully with NaNs disobeying the ordering rules
    (and even the equivalence rules for = and =/=). **)
trait Number
        extends { StandardPartialOrder[\Number\], StandardMin[\Number\], StandardMax[\Number\],
                  AdditiveGroup[\Number\], MultiplicativeRing[\Number\] }
        comprises { RR64 }
    getter zero(): Number = 0.0
    getter one(): Number = 1.0
    getter toString(): String
    opr |self| : Number =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Abs")
    opr =(self, b:Number):Boolean =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Eq")
    opr =/=(self, b:Number):Boolean =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$NEq")
    opr <(self, b:Number):Boolean =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Less")
    opr <=(self, b:Number):Boolean =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$LessEq")
    opr >(self, b:Number):Boolean =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Greater")
    opr >=(self, b:Number):Boolean =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$GreaterEq")
    opr CMP(self, b:Number):Comparison =
        if self<b then LessThan
        elif self>b then GreaterThan
        elif self=b then EqualTo
        else Unordered
        end
    (** In case of NaN, %MIN% and %MAX% return a NaN, otherwise it respects the
        total order. **)
    opr MIN(self, b:Number):Number =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Min")
    opr MAX(self, b:Number):Number =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Max")

    opr -(self):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Negate")
    opr +(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Add")
    opr -(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Sub")
    opr DOT(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Mul")
    opr TIMES(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Mul")
    opr juxtaposition
         (self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Mul")
    opr /(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Div")
    opr SQRT(self):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Sqrt")
    opr PLUS_UP(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$AddUpNoNaN")
    opr MINUS_UP(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$SubUpNoNaN")
    opr DOT_UP(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$MulUpNoNaN")
    opr SLASH_UP(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$DivUpNoNaN")
    opr SQRT_UP(self):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$SqrtUp")
    opr PLUS_DOWN(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$AddDownNoNaN")
    opr MINUS_DOWN(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$SubDownNoNaN")
    opr DOT_DOWN(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$MulDownNoNaN")
    opr SLASH_DOWN(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$DivDownNoNaN")
    opr SQRT_DOWN(self):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$SqrtDown")
    opr IEEE_PLUS_UP(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$AddUp")
    opr IEEE_MINUS_UP(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$SubUp")
    opr IEEE_DOT_UP(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$MulUp")
    opr IEEE_SLASH_UP(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$DivUp")
    opr IEEE_PLUS_DOWN(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$AddDown")
    opr IEEE_MINUS_DOWN(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$SubDown")
    opr IEEE_DOT_DOWN(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$MulDown")
    opr IEEE_SLASH_DOWN(self,b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$DivDown")
    opr ^(self, b:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Pow")
    (** Shouldn't need this extra declaration. **)
    opr ^(self, b:AnyIntegral):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Pow")
    sin(self):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Sin")
    cos(self):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Cos")
    tan(self):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Tan")
    asin(self):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$ASin")
    acos(self):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$ACos")
    atan(self):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$ATan")
    atan2(self,x:Number):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$ATan2")
    log(self):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Log")
    exp(self):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Exp")
    floor(self):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Floor")
    opr |\self/| : ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$IFloor")
    ceiling(self):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Ceiling")
    opr |/self\| : ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$ICeiling")
    truncate(self):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Truncate")
end

trait RR64 extends Number comprises { Float, AnyIntegral, FloatLiteral }
    (** returns true if the value is an IEEE NaN **)
    getter isNaN(): Boolean =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$isNaN")
    (** returns true if the value is an IEEE infinity **)
    getter isInfinite(): Boolean =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$isInfinite")
    (** returns true if the value is a valid number (not NaN) **)
    getter isNumber(): Boolean = NOT self.isNaN
    (** returns true if the value is finite **)
    getter isFinite(): Boolean = NOT (self.isInfinite OR self.isNaN)
    (** check returns Just(its argument) if it is finite, otherwise Nothing. **)
    getter check(): Maybe[\RR64\] =
        if self.isFinite then Just[\RR64\](self) else Nothing[\RR64\] end
    (** check_star returns Just(its argument) if it is non-NaN, otherwise Nothing. **)
    getter check_star(): Maybe[\RR64\] =
        if self.isNaN then Nothing[\RR64\] else Just[\RR64\](self) end
    (** obtain the raw bits of the IEEE floating-point representation of this value. **)
    getter rawBits():ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$RawBits")
    getter signBit():ZZ32 = if rawBits() < 0 then 1 else 0 end
    (** next higher IEEE float **)
    getter nextUp():RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$NextUp")
    (** next lower IEEE float **)
    getter nextDown():RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$NextDown")
    (** %MINNUM% and %MAXNUM% return a numeric result where possible (avoiding NaN).
        Note that %MINNUM% and %MAX% form a lattice with NaN at the top, and
        that %MAXNUM% and %MIN% form a lattice with NaN at the bottom.  **)
    opr MINNUM(self, b:RR64):RR64 = do
        r = self MIN b
        if r.isNaN then
            if self.isNumber then self
            elif b.isNumber then b
            else r end
        else
            r
        end
      end
    opr MAXNUM(self, b:RR64):RR64 = do
        r = self MAX b
        if r.isNaN then
            if self.isNumber then self
            elif b.isNumber then b
            else r end
        else
            r
        end
      end
    (** returns a value of type RR32 **)
    narrow(self): RR32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Narrow")
end

trait QQ extends { RR64, StandardPartialOrder[\QQ\] }
    getter zero(): QQ = 0
    getter one(): QQ = 1
    getter isNaN(): Boolean
    getter isInfinite(): Boolean
    getter isNumber(): Boolean = NOT self.isNaN
    getter isFinite(): Boolean = (denominator(self) =/= 0)
    numerator(self): ZZ
    denominator(self): ZZ  (* Ideally would be NN *)
    opr |self| : QQ = |numerator(self)| / denominator(self)
    (* Note that 0/0 = 0/0 is supposed to be true (//unlike// the behavior of floating-point NaN) *)
    opr =(self, other:QQ):Boolean =
        numerator(self) = numerator(other) AND: denominator(self) = denominator(other)
    opr =/=(self, other:QQ):Boolean = NOT (self = other)
    opr <(self, other:QQ):Boolean = do
        (a, b) = (numerator(self), denominator(self))
        (c, d) = (numerator(other), denominator(other))
        (b=0 AND d=0 AND a=-1 AND c=1) OR (a d < b c)
    end
    opr <=(self, other:QQ):Boolean = (self < other) OR (self = other)
    opr >(self, other:QQ):Boolean = (other < self)
    opr >=(self, other:QQ):Boolean = (other <= self)
    (** In case of 0/0, %MIN% and %MAX% return 0/0, otherwise it respects the
        total order. **)
    opr MIN(self, other:QQ):QQ =
        if self < other then self elif other < self then other elif self = 0/0 then self else other end
    opr MAX(self, other:QQ):QQ =
        if self > other then self elif other > self then other elif self = 0/0 then self else other end

    opr -(self):QQ = (-numerator(self)) / denominator(self)
    opr +(self,other:QQ):QQ =
        if self = other AND denominator(self) = 0 then self else
            (a, b) = (numerator(self), denominator(self))
            (c, d) = (numerator(other), denominator(other))
            (a d + b c) / b d
        end
    opr -(self,other:QQ):QQ = self + (-other)
    opr DOT(self,other:QQ):QQ = do
        (a, b) = (numerator(self), denominator(self))
        (c, d) = (numerator(other), denominator(other))
        a c / b d
    end
    opr TIMES(self,other:QQ):QQ = self DOT other
    opr juxtaposition (self,other:QQ):QQ = self DOT other
    opr /(self,other:QQ):QQ = do
        (a, b) = (numerator(self), denominator(self))
        (c, d) = (numerator(other), denominator(other))
        a d / b c
    end
    opr ^(self, other:AnyIntegral):QQ =
        if -other > 0 then   (* Way too clever by half (to avoid infinite recursion on 2^(-63)) *)
            (denominator(self)/numerator(self))^(-other)
        else
            (numerator(self))^other / (denominator(self))^other
        end
    floor(self):ZZ = if self < 0 then -ceiling(-self) else truncate(self) end
    opr |\self/| : ZZ = floor(self)
    ceiling(self):ZZ = if denominator(self) = 0 then self else
        (a, b) = (numerator(self), denominator(self))
        if a < 0 then -floor(-self) else (a + b - 1) DIV b end
    end
    opr |/self\| : ZZ = ceiling(self)
    truncate(self):ZZ = if denominator(self) = 0 then self else numerator(self) DIV denominator(self) end
    round(self): ZZ = do x = self + 1/2; z = floor(x); if z = x AND odd z then z-1 else z end end
    opr MINNUM(self, other:QQ):QQ =
        if self < other then self elif other < self then other elif self = 0/0 then other else self end
    opr MAXNUM(self, other:QQ):QQ =
        if self > other then self elif other > self then other elif self = 0/0 then other else self end
end

object Ratio(num: ZZ, den: ZZ) extends { QQ }
    getter toString(): String = num.toString || "/" || den.toString
    getter isNaN(): Boolean = (numerator(self) = 0 AND denominator(self) = 0)
    getter isInfinite(): Boolean = (numerator(self) =/= 0 AND denominator(self) = 0)
    numerator(self):ZZ = num
    denominator(self):ZZ = den
end

trait AnyIntegral extends { QQ } end

trait Integral[\I extends Integral[\I\]\] extends { StandardTotalOrder[\I\], AnyIntegral }
    getter zero(): I
    getter one(): I
    getter isNaN(): Boolean = false
    getter isInfinite(): Boolean = false  (* For now, no infinite integers. *)
    denominator(self):ZZ = 1
    opr -(self):I
    opr +(self,b:I):I
    opr -(self,b:I):I
    opr DOT(self,b:I):I
    opr TIMES(self,b:I):I
    opr juxtaposition(self,b:I):I
    opr DIV(self,b:I):I
    opr REM(self,b:I):I
    opr MOD(self,b:I):I
    opr GCD(self,b:I):I
    opr LCM(self,b:I):I
    opr DIVIDES(self,b:I):Boolean = b =/= zero() AND: (b MOD self) = zero()
    opr CHOOSE(self,b:I):I
    opr BITAND(self,b:I):I
    opr BITOR(self,b:I):I
    opr BITXOR(self,b:I):I
    opr LSHIFT(self,b:AnyIntegral):I
    opr RSHIFT(self,b:AnyIntegral):I
    opr BITNOT(self):I
    opr ^(self, b:AnyIntegral):RR64
end

trait ZZ32 extends { ZZ64, Integral[\ZZ32\] } comprises { Int, IntLiteral }
    getter zero(): ZZ32 = 0
    getter one(): ZZ32 = 1
    getter minimum(): ZZ32 = -2147483647 - 1
    getter maximum(): ZZ32 = 2147483647
        
    opr |self| : ZZ32 = if self>=0 then self else -self end
    opr =(self, b:ZZ32):Boolean =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$Eq")
    opr <(self, b:ZZ32):Boolean =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$Less")

    opr -(self):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$Negate")
    opr +(self,b:ZZ32):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$Add")
    opr -(self,b:ZZ32):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$Sub")
    opr DOT(self,b:ZZ32):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$Mul")
    opr TIMES(self,b:ZZ32):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$Mul")
    opr juxtaposition(self,b:ZZ32):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$Mul")
    opr DIV(self,b:ZZ32):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$Div")
    opr REM(self,b:ZZ32):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$Rem")
    opr MOD(self,b:ZZ32):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$Mod")
    opr GCD(self,b:ZZ32):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$Gcd")
    opr LCM(self,b:ZZ32):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$Lcm")
    opr CHOOSE(self,b:ZZ32):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$Choose")
    opr BITAND(self,b:ZZ32):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$BitAnd")
    opr BITOR(self,b:ZZ32):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$BitOr")
    opr BITXOR(self,b:ZZ32):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$BitXor")
    opr LSHIFT(self,b:ZZ64):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$LShift")
    opr RSHIFT(self,b:ZZ64):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$RShift")
    opr BITNOT(self):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$BitNot")
    opr ^(self, b:AnyIntegral):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$Pow")
    widen(self):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$ToLong")
    partitionL(self):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Int$Partition")
    unsigned(self):NN32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.ZZ32$ToNN32")
end

trait ZZ64 extends { ZZ, Integral[\ZZ64\] } comprises { Long, ZZ32 }
    getter zero(): ZZ64 = widen(0)
    getter one(): ZZ64 = widen(1)
    getter minimum(): ZZ64 = -9223372036854775807 - 1
    getter maximum(): ZZ64 = 9223372036854775807

    opr |x:ZZ64| : ZZ64 = if x>=0 then x else -x end

    opr =(self, b:ZZ64):Boolean =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$Eq")
    opr <(self, b:ZZ64):Boolean =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$Less")
    (* Argh!  Due to method ambiguities with ZZ, these definitions must be given
       explicitly here. *)
    opr >(self, b:ZZ64):Boolean = b < self
    opr >=(self, b:ZZ64):Boolean = NOT (self < b)
    opr <=(self, b:ZZ64):Boolean = NOT (b < self)
    opr CMP(self, b:ZZ64): TotalComparison =
        if self < b then LessThan
        elif b < self then GreaterThan
        else EqualTo end

    opr -(self):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$Negate")
    opr +(self,b:ZZ64):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$Add")
    opr -(self,b:ZZ64):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$Sub")
    opr DOT(self,b:ZZ64):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$Mul")
    opr TIMES(self,b:ZZ64):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$Mul")
    opr juxtaposition(self,b:ZZ64):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$Mul")
    opr DIV(self,b:ZZ64):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$Div")
    opr REM(self,b:ZZ64):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$Rem")
    opr MOD(self,b:ZZ64):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$Mod")
    opr GCD(self,b:ZZ64):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$Gcd")
    opr LCM(self,b:ZZ64):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$Lcm")
    opr CHOOSE(self,b:ZZ64):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$Choose")
    opr BITAND(self,b:ZZ64):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$BitAnd")
    opr BITOR(self,b:ZZ64):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$BitOr")
    opr BITXOR(self,b:ZZ64):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$BitXor")
    opr LSHIFT(self,b:AnyIntegral):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$LShift")
    opr RSHIFT(self,b:AnyIntegral):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$RShift")
    opr BITNOT(self):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$BitNot")
    opr ^(self, b:AnyIntegral):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$Pow")
    partitionL(self):ZZ64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$Partition")
    narrow(self):ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$FromLong")
    unsigned(self):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$ToUnsignedLong")
    big(self):ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$ToBigNum")
end

trait NN64 extends { ZZ, Integral[\NN64\] } comprises { UnsignedLong, NN32 }
    opr |self| : NN64 = self
    opr =(self, b:NN64):Boolean =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$Eq")
    opr <(self, b:NN64):Boolean =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$Less")
    opr -(self):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$Negate")
    opr +(self,b:NN64):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$Add")
    opr -(self,b:NN64):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$Sub")
    opr DOT(self,b:NN64):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$Mul")
    opr TIMES(self,b:NN64):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$Mul")
    opr juxtaposition(self,b:NN64):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$Mul")
    opr DIV(self,b:NN64):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$Div")
    opr REM(self,b:NN64):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$Rem")
    opr MOD(self,b:NN64):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$Mod")
    opr GCD(self,b:NN64):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$Gcd")
    opr LCM(self,b:NN64):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$Lcm")
    opr CHOOSE(self,b:NN64):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$Choose")
    opr BITAND(self,b:NN64):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$BitAnd")
    opr BITOR(self,b:NN64):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$BitOr")
    opr BITXOR(self,b:NN64):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$BitXor")
    opr LSHIFT(self,b:AnyIntegral):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$LShift")
    opr RSHIFT(self,b:AnyIntegral):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$RShift")
    opr BITNOT(self):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$BitNot")
    opr ^(self, b:AnyIntegral):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$Pow")
    partitionL(self):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$Partition")
    narrow(self):NN32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$FromLong")
    signed(self):NN64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.UnsignedLong$ToLong")
end

trait ZZ extends Integral[\ZZ\]
        comprises { BigNum, ZZ64 }
    getter zero(): ZZ =
        1934791870947204798109283471902037419 - 1934791870947204798109283471902037419
    getter one(): ZZ = zero() + 1
    numerator(self):ZZ = self
    opr |self| : ZZ32 = if self>=0 then self else -self end
    opr =(self, b: ZZ):Boolean =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$Eq")
    opr <(self, other:ZZ): Boolean = self.cmp(other) < 0
    opr <=(self, other:ZZ): Boolean = self.cmp(other) <= 0
    opr >(self, other:ZZ): Boolean = self.cmp(other) > 0
    opr >=(self, other:ZZ): Boolean = self.cmp(other) >= 0
    opr CMP(self, other:ZZ): TotalComparison = self.cmp(other) CMP 0
    cmp(b:ZZ): ZZ32 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$Cmp")

    opr -(self): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$Negate")
    opr +(self, b: ZZ): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$Add")
    opr -(self, b: ZZ): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$Sub")
    opr DOT(self, b: ZZ): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$Mul")
    opr juxtaposition(self, b: ZZ): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$Mul")
    opr TIMES(self, b: ZZ): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$Mul")
    opr /(self,other:ZZ):QQ = if other < 0 then (-self)/(-other) else
        g = self GCD other
        if g=0 then Ratio(0,0) else
            a = self DIV g
            b = other DIV g
            if b = 1 then a else Ratio(a,b) end
        end
    end
    opr DIV(self, b: ZZ): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$Div")
    opr REM(self, b: ZZ): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$Rem")
    opr MOD(self, b: ZZ): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$Mod")
    opr GCD(self, b: ZZ): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$Gcd")
    opr LCM(self, b: ZZ): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$Lcm")
    opr CHOOSE(self, b: ZZ): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$Choose")
    opr BITAND(self, b: ZZ): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$BitAnd")
    opr BITOR(self, b: ZZ): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$BitOr")
    opr BITXOR(self, b: ZZ): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$BitXor")
    opr LSHIFT(self, b:AnyIntegral): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$LShift")
    opr RSHIFT(self, b:AnyIntegral): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$RShift")
    opr BITNOT(self): ZZ =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$BitNot")
    opr ^(self, b:AnyIntegral):RR64 =
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.BigNum$Pow")
    odd(self): Boolean = NOT (even self)
    even(self): Boolean = (self MOD 2) = 0
end

(************************************************************
* Generator support
************************************************************)

(** Generator
 *
 * We say an object which extends Generator[\T\] "generates objects of
 * type T."
 *
 * Generators are used to express iteration in Fortress.  Every
 * generated expression in Fortress (eg for loop, comprehension) is
 * desugared into calls to methods of Generator, chiefly the generate
 * method.
 *
 * Every generator has a notion of a "natural order" (which by default is
 * unspecified), which describes the ordering of reduction operations,
 * and also describes the order in which elements are produced by the
 * sequential version of the same generator (given by the seq(self)
 * method).  The default implementation of seq(self) guarantees that these
 * orders will match.
 *
 * Note in particular that the natural order of a Generator must be
 * consistent; that is, if a SEQV b then a and b must yield SEQV
 * elements in the same natural order.  However, note that unless a type
 * specifically documents otherwise, no particular element ordering is
 * guaranteed, nor is it necessary to guarantee that a=b have the same
 * natural order when equality is defined.
 *
 * Note that more complex derived generators are specified further
 * down in the definition of Generator.  These have the same notions
 * of natural order and by default are defined in terms of the
 * generate() method.
 *
 * Minimal complete definition of a Generator is the generate(...) method.
 *)
trait Generator[\E\]
    excludes { Number }
    (** generate is the core of Generator.  It generates elements of
        type E and passes them to the body function.  This generation
        can occur using any mixture of serial and parallel execution
        desired by the author of the generator; by default uses of a
        generator must assume every call to the body occurs in
        parallel.

        The results of generation are combined using the reduction
        object R, which specifies a monoidal operation (associative
        and with an identity).  Body results must be combined together
        following the natural order of the generator.  The author of
        the generator is free to use the identity element anywhere
        desired in this computation, and to group reductions in any
        way desired; if no elements are generated the identity must be
        returned. *)
    generate[\R\](r: Reduction[\R\], body: E->R): R

    (** Transforming generators into new generators *)
    (** map applies a function f to each element generated and yields
        the result.  The resulting generator must have the same
        ordering and cross product properties as the generator from
        which it is derived. *)
    map[\G\](f: E->G): Generator[\G\] = SimpleMappedGenerator[\E,G\](self,f)
    (** seq produces a sequential version of the same generator, in
        which elements are produced in the generator's natural order. *)
    seq(self): SequentialGenerator[\E\] = NaiveSeqGenerator[\E\](self)

    (** Nesting of two generators; the innermost is data-dependent
        upon the outer one.  This is specifically designed to be
        overloaded so that the combined generators have properties
        appropriate to the pairing.  Because of the data dependency
        the natural order of the nesting is the natural order of the
        inner generators, in the natural order the outer nesting
        produces them.  So for example, if we write:
          (0#3).nest[\ZZ32\](\(n:ZZ32):Generator[\ZZ32\] => (n*100#4))
        then the natural order is 0,1,2,3,100,101,102,103,200,201,202,203
     **)
    nest[\G\](f: E -> Generator[\G\]): Generator[\G\] =
        SimpleNestedGenerator[\E,G\](self,f)

    (** Filtering data from a generator.  Only elements that satisfy
        the predicate p are retained.  Natural order and cross product
        properties are otherwise preserved. **)
    filter(f: E -> Condition[\()\]): Generator[\E\] =
        SimpleFilterGenerator[\E\](self,f)

    (** Cross product of two generators.  This is specifically
        designed to be overloaded, such that pairs of independent
        generators can be combined to produce a generator which
        possibly interleaves the iteration spaces of the two
        generators.  For example, we might combine (0#16).cross(0#32)
        such that it first splits the second range in half, then the
        first range in half, then the second, and so forth.

        Consider a grid for which the rows are labeled from top to
        bottom with the elements of a in natural order, and the
        columns are labeled from left to right with the elements of g
        in natural order.  Each point in the grid corresponds to a
        pair (a,b) that must be generated by self.cross(g).  In the
        natural order of the cross product, an element must occur
        after those that lie above and to the left of it in the grid.
        By default the natural order of the cross product is
        left-to-right, top to bottom.  Programmers must not rely on
        the default order, except that cross products involving one or
        more sequential generators are always performed in the default
        order.  Note that this means that the following have the same
        natural order:
           seq(a).cross(b)
           a.cross(seq(b))
           seq(a).cross(seq(b))
        But seq(a.cross(b)) may have a different natural order. *)
    cross[\G\](g: Generator[\G\]): Generator[\(E,G)\] =
        SimplePairGenerator[\E,G\](self, g)

    (** Derived generation operations *)
    (** mapReduce is equivalent to generate, but takes an explicit join
        and zero which can have any type.  It still assumes join is
        associative and that zero is the identity of join. **)
    mapReduce[\R\](body: E->R, join:(R,R)->R, id:R): R =
        generate[\R\](MapReduceReduction[\R\](join,id), body)
    (** reduce works much like generate or mapReduce,
        but has no body expression **)
    reduce(j:(E,E)->E, z:E):E = mapReduce[\E\](fn (e:E)=>e, j, z)
    reduce(r: Reduction[\E\]):E = generate[\E\](r, fn(e:E)=>e)
    (** loop is a version of generate which discards the void results
        of the body computation.  It is used to translated
        reduction-variable-free for loops. **)
    loop(f:E->()): () = generate[\()\](VoidReduction, f)
    (** Is x generated by this generator? **)
    opr IN(x:E, self): Boolean = BIG OR [y <- self] x=y
    (** The MATCH operator is used (as a temporary hack) by CASE when the
        match clause is a Generator[\E\].  We must choose either =, SEQV, or
        IN as appropriate, depending on the type of the LHS. **)
    opr MATCH(x:Any, self): Boolean =
        typecase x of
            E => (x IN self)
            Generator[\E\] => (x = self)
            else => fail("Type error in MATCH of Generator")
        end
end

(** The following top-level functions are used as glue in the
    desugaring of generated expressions (comprehensions, BIG
    operators, and loops).  They exist at the top level in part to
    work around shortcomings in the ablity of the Fortress interpreter
    to check types dynamically for anything except top-level
    non-overloaded function calls. **)
__generate[\E,R\](g:Generator[\E\], r: Reduction[\R\], b:E->R): R =
    g.generate[\R\](r,b)

__nest[\E1,E2\](g:Generator[\E1\], f:E1->Generator[\E2\]):Generator[\E2\] =
    g.nest[\E2\](f)

__singleton[\E\](body:E): Generator[\E\] = Just[\E\](body)

__loop[\E,R\](g:Generator[\E\], body:E->R): () = g.loop(body)

__cond[\E,R\](c:Condition[\E\], t:E->R, e:()->R): R = c.cond[\R\](t,e)
__cond[\E\](c:Condition[\E\], t:E->()): () = c.cond[\()\](t, fn () => ())

(** Used in desugaring binding %while% **)
__whileCond[\E\](c:Condition[\E\], b:E->()): () =
    c.cond[\Boolean\](fn (e): Boolean => do b(e); true end, fn ():Boolean => false)

(** Unlike the user-visible dotted map method, this method is used for
    the internals of comprehension desugaring [except that at present
    it is not actually in use, but we're planning to use it].  As a
    result we know that its result will immediately be consumed, and
    we should *not* create an intermediate collection if we're passed
    a collection. **)
__map[\E,R\](g:Generator[\E\], f:E->R): Generator[\R\] =
    typecase g of
        SomeMappedGenerator[\E\] => g.map[\R\](f)
        SequentialGenerator[\E\] => SimpleMappedSeqGenerator[\E,R\](g,f)
        else => SimpleMappedGenerator[\E,R\](g,f)
    end

(** Unlike the user-visible dotted filter method, but like the __map
    method, this method is intended for use in the internals of
    comprehension desugaring.  Again its result will immediately be
    consumed, and we do not create an intermediate collection but
    perform the filtering lazily on demand. **)
__filter[\E\](g:Generator[\E\], p:E->Condition[\()\]): Generator[\E\] =
    typecase g of
        FilterGenerator[\E\] => g.filter(p)
        SequentialGenerator[\E\] => SimpleSeqFilterGenerator[\E\](g,p)
        else => SimpleFilterGenerator[\E\](g,p)
    end

(** Application of a BIG operator, the topmost level of comprehension desugaring. **)
__bigOperator[\I,O,R,L\](o:BigOperator[\I,O,R,L\],desugaredClauses:(Reduction[\L\],I->L)->L): O = do
    r = o.reduction()
    body(i): L = r.lift((o.body())(i))
    (o.unwrap())(r.unlift(desugaredClauses(r,body)))
  end

(** Application of two nested BIG operators, possibly with fusion.  This only covers
    a comprehension of the form:
       BIG outer [ xs <- expr_o ] (BIG inner [x <- xs] expr_i)
    The desugarer extracts comprehensions of this form from more complex nests of
    comprehensions, using a combination of splitting:
        BIG OP [ gs_1, gs_2 ] expr = BIG OP [ gs_1 ] (BIG OP [ gs_2 ] expr)
    and filter squeezing:
        BIG OP [ xs <- g, p(xs) ] expr = BIG OP [ xs <- g.filter(p) ] expr
    There are some big caveats to this explanation in practice.  Most important is that
    we don't unlift and lift or do input/output conversion except where neccessary, so
    splitting skips these operations in between the inner and outer comprehension.
    **)
__bigOperator2[\I,M,O,R1,L1,R2,L2,E\](outer:BigOperator[\M,O,R1,L1\],
                                    inner:BigOperator[\I,M,R2,L2\],
                                    gg: Generator[\Generator[\E\]\],
                                    innerBody:E->L2): L1 = do
    o = outer.reduction()
    i = inner.reduction()
    gg.generate(o,fn (xs) =>
                   o.lift((outer.wrap())((inner.unwrap())(i.unlift(
                                                xs.generate(i,innerBody))))))
  end

trait SequentialGenerator[\E\] extends { Generator[\E\] }
    seq(self) = self
    map[\G\](f: E->G): SequentialGenerator[\G\] =
        SimpleMappedSeqGenerator[\E,G\](self,f)
    (* TODO: make overloaded *)
    nest[\G\](f: E -> Generator[\G\]): Generator[\G\] =
        typecase f of
            E -> SequentialGenerator[\G\] =>
                SimpleNestedSeqGenerator[\E,G\](self,f)
            else => SimpleNestedGenerator[\E,G\](self,f)
        end
    (* TODO: make overloaded *)
    cross[\F\](g:Generator[\F\]): Generator[\(E,F)\] =
        typecase g of
            SequentialGenerator[\F\] =>
                    SimplePairSeqGenerator[\E,F\](self,g)
            else => SimplePairGenerator[\E,F\](self,g)
        end
end

(** A Condition is a Generator that generates 0 or 1 element.
    Conditions can be used as nullary comprehension generators or
    as predicates in an if expression. **)
trait Condition[\E\] extends { ZeroIndexed[\E\], SequentialGenerator[\E\] }
    getter isEmpty(): Boolean = NOT holds()
    getter nonEmpty(): Boolean = NOT isEmpty()
    getter holds(): Boolean = cond[\Boolean\](fn (_:E):Boolean => true, fn () => false)
    getter size(): ZZ32 = if holds() then 1 else 0 end
    getter get(): E throws NotFound = cond[\E\](identity[\E\], fn () => throw NotFound)
    getter bounds(): FullRange[\ZZ32,true\] = 0 # |self|
    getter indices(): FullRange[\ZZ32,true\] = bounds()
    getter indexValuePairs(): Condition[\(ZZ32,E)\] =
        cond[\E\](fn (e:E) => Just[\(ZZ32,E)\](0,e), fn () => Nothing[\(ZZ32,E)\])
    opr |self|: ZZ32 = if holds() then 1 else 0 end
    opr [i:ZZ32]:E throws NotFound = if i=0 then get() else throw NotFound end
    opr [r:Range[\ZZ32\]]: Condition[\E\] =
        cond[\E\](fn (e:E): Condition[\E\] => do
                    r1 = (0#1)[r]
                    if r1.isEmpty then Nothing[\E\] else self end
                  end,
                fn (): Condition[\E\] => do r1 = (0#0)[r]; self end)
    getDefault(e:E): E = cond[\E\](identity[\E\], fn () => e)
    cond[\G\](t: E -> G, e: () -> G): G
    generate[\G\](r:Reduction[\G\], body: E -> G): G = cond[\G\](body,fn () => r.empty())

    (** For a Condition, these methods run eagerly. **)
    map[\G\](f: E->G): Condition[\G\] =
        cond[\G\](fn (e:E) => Just[\G\](f(e)), fn () => Nothing[\G\])
    ivmap[\G\](f: (ZZ32,E)->G): Condition[\G\] =
        cond[\G\](fn (e:E) => Just[\G\](f(0,e)), fn () => Nothing[\G\])
    nest[\G\](f: E -> Generator[\G\]): Generator[\G\] =
        cond[\G\](f, fn () => Nothing[\G\])
    cross[\G\](g: Generator[\G\]): Generator[\(E,G)\] =
        cond[\G\](fn (e:E):Generator[\(E,G)\] => g.map[\(E,G)\](fn (g':G):(E,G) => (e,g')),
               fn () => Nothing[\(E,G)\])
    mapReduce[\R\](body: E->R, _:(R,R)->R, id:R): R =
        cond[\R\](body, fn () => id)
    reduce(_:(E,E)->E, z:E):E = cond[\E\](identity[\E\], fn ():E => z)
    reduce(r: Reduction[\E\]):E = cond[\E\](identity[\E\], fn ():E => r.empty())
    loop(f:E->()): () = cond[\()\](f,fn ():() => ())
    opr IN(x:E, self):Boolean = cond[\Boolean\](fn (e:E):Boolean => x=e, fn () => false)
end

(** Conjunction and disjunction of condition functions **)
opr ANDCOND[\I\](p1:I->Condition[\()\], p2:I->Condition[\()\]): I->Condition[\()\] =
    fn (i:I):Boolean => p1(i).holds AND: p2(i).holds

(** Conjunction and disjunction of condition functions **)
opr ORCOND[\I\](p1:I->Condition[\()\], p2:I->Condition[\()\]): I->Condition[\()\] =
    fn (i:I):Boolean => p1(i).holds OR: p2(i).holds

(** Operations which use generation internally.  Should be functional
    methods of Generator, but that didn't work when last tested. **)

(** IN returns true if any element generated by its second argument is
    = to its first argument.  x NOTIN g is simply NOT (x IN g).
    Unless documented otherwise, this is O(n) for an n-element
    Generator (it simply performs naive matching).  **)
opr NOTIN[\E\](x: E, this: Generator[\E\]): Boolean = NOT (x IN this)

(* seq[\T\](g:Generator[\T\]):SequentialGenerator[\T\] = g.seq() *)
sequential[\T\](g:Generator[\T\]):SequentialGenerator[\T\] = seq(g)

(************************************************************
* Meta-generators: BIG GENERATOR
************************************************************)

(** The problem: Your code is sprinkled with the same generator code,
    over and over again.  Maybe it looks something like this:
    % for (i,v) <- a.indexValuePairs, v =/= 0, j <- 0:i do (code using i and j) %
    You goal: abstract this into a single generator.  You can do this as follows:
    % usefulGen = BIG GENERATOR[(i,v) <- a.indexValuePairs, v =/= 0, j <- 0:i] (i,j) %
    Not much savings there, but the important part is that you can just
    cut and paste the generators from the original code.  And all occurrences
    can be cleaned up:
    % for (i,j) <- usefulGen do (code using i and j) %
    This is possible using %map%, %filter%, and %nest%, but you end up
    effectively desugaring the code by hand.  Using BIG GENERATOR is much easier.
 **)

(** TODO

The big problem here is that the new desugaring story actually gets in
the way.  We can't simply package up a function passed in, because
that's being done by the __bigOperator functionality directly.

How do we deal with reduction properties?  It's all very sticky.  And
we're already running into difficulties with co- and contra-variance
of the type parameters floating around that force us to use Any in
awkward places.

In particular, we don't want to generate the data when creating the
generator, as that misses the whole point.  We want to defer the
object traversals until the moment the resulting generator is used.
**)

(************************************************************
* The Maybe type, used instead of null
************************************************************)

(* This makes excludes work without where clauses, and allows opr =()
   to remain non-parametric. *)
value trait AnyMaybe extends Equality[\AnyMaybe\] excludes Number
        (** \vspace{-4ex} NOT YET: %comprises Maybe[\T\] where [\T\]% *)
    getter holds() : Boolean
    opr =(self, other:AnyMaybe): Boolean = false
end

just(t:Any):AnyMaybe = Just(t)

(** %Maybe% represents either %Nothing% or a single element of type
    %T% (%Just[\T\]%), which may be retrieved by calling %get%.  An
    object of type %Maybe[\T\]% can be used as a generator; it is either
    empty (%Nothing%) or generates the single element yielded by
    %get%, so there is no issue of canonical order or parallelism.

    Thus, %Just[\T\]% can be used as a single-element generator, and
    %Nothing% can be used as an empty generator. *)
value trait Maybe[\T\]
        extends { AnyMaybe, Condition[\T\] }
        comprises { Nothing[\T\], Just[\T\] }
end

value object Just[\T\](x:T) extends Maybe[\T\]
    getter size()=1
    getter toString():String = "Just(" x ")"
    getter holds() = true
    getter get() = x
    opr |self| : ZZ32 = 1
    getDefault(_:T): T = x
    cond[\R\](t:T->R, _:()->R): R = t(x)
    generate[\R\](_:Reduction[\R\],m:T->R): R = m(x)
    opr[i:ZZ32]:T = if i=0 then x else fail("Maybe[" i "] nonzero index") end

    map[\G\](f: T->G): Just[\G\] = Just[\G\](f(x))
    cross[\G\](g: Generator[\G\]): Generator[\(T,G)\] =
        g.map[\(T,G)\](fn (e:G):(T,G) => (x,e))

    mapReduce[\R\](m: T->R, _:(R,R)->R, _:R): R = m(x)
    reduce(_:(T,T)->T, _:T):T = x
    reduce(_: Reduction[\T\]):T = x
    loop(f:T->()): () = f(x)
    opr IN(y:T, self): Boolean = x=y
    opr =(self,o:Just[\T\]): Boolean = x=o.get
end

(* Obviously ought to be a non-parametric singleton when we get where
   clauses working. *)
value object Nothing[\T\] extends Maybe[\T\]
    getter size()=0
    getter holds() = false
    getter get() = throw NotFound
    getter toString():String = "Nothing"
    opr |self| : ZZ32 = 0
    getDefault(t:T):T = t
    cond[\R\](_:T->R, e:()->R): R = e()
    generate[\R\](r:Reduction[\R\],_:T->R): R = r.empty()
    opr[_:ZZ32]: T = fail("Cannot index Nothing")

    map[\G\](f: T->G): Nothing[\G\] = Nothing[\G\]
    cross[\G\](g: Generator[\G\]): Generator[\(T,G)\] = Nothing[\(T,G)\]

    mapReduce[\R\](_: T->R, _:(R,R)->R, z:R): R = z
    reduce(_:(T,T)->T, z:T):T = z
    reduce(r: Reduction[\T\]):T = r.empty()
    loop(f:T->()): () = ()
    opr IN(x:T, self): Boolean = false
    opr =(self,_:Nothing[\T\]) = true
end

(************************************************************
* Exception hierarchy
************************************************************)

trait Exception comprises { UncheckedException, CheckedException }
end

(* Exceptions which are not checked *)

trait UncheckedException extends Exception excludes CheckedException
end

object FailCalled(s:String) extends UncheckedException
  toString(): String = "FAIL: " s
end

object DivisionByZero extends UncheckedException
    toString(): String = "Division by zero"
end

object UnpastingError extends UncheckedException
    toString(): String = "Unpasting error"
end

object CallerViolation extends UncheckedException
    toString(): String = "Caller violation"
end

object CalleeViolation extends UncheckedException
    toString(): String = "Callee violation"
end

object TestFailure extends UncheckedException
    toString(): String = "Test failure"
end

object ContractHierarchyViolation extends UncheckedException
    toString(): String = "Contract hierarchy violation"
end

object NoEqualityOnFunctions extends UncheckedException
    toString(): String = "No equality on functions"
end

object InvalidRange extends UncheckedException
    toString(): String = "Invalid range"
end

object ForbiddenException(chain : Exception) extends UncheckedException
    toString(): String = "Forbidden exception"
end

(* Should this be called "IndexNotFound" instead? *)
object NotFound extends UncheckedException
    toString(): String = "Not found"
end

object IndexOutOfBounds(min:ZZ32,max:ZZ32,index:ZZ32) extends UncheckedException
    toString(): String = index " is outside the range " min " to " max
end

object NegativeLength extends UncheckedException
    toString(): String = "Negative length"
end

object IntegerOverflow extends UncheckedException
    toString(): String = "Integer Overflow"
end

object RationalComparisonError extends UncheckedException
    toString(): String = "Rational comparison error"
end

object FloatingComparisonError extends UncheckedException
    toString(): String = "Floating comparison error"
end

(* Checked Exceptions *)

trait CheckedException extends Exception excludes UncheckedException
end

object CastError extends CheckedException
    toString(): String = "Cast error"
end

object IOFailure extends CheckedException
    toString(): String = "I/O error"
end

object MatchFailure extends CheckedException
    toString(): String = "Match failure"
end

(* SetsNotDisjoint? *)
object DisjointUnionError extends CheckedException
    toString(): String = "Disjoint union error"
end

object APIMissing extends CheckedException
    toString(): String = "Api is missing"
end

object APINameCollision extends CheckedException
    toString(): String = "Api name collides with another"
end

object ExportedAPIMissing extends CheckedException
    toString(): String = "Exported api is missing"
end

object HiddenAPIMissing extends CheckedException
    toString(): String = "Hidden api is missing"
end

object TryAtomicFailure extends CheckedException
    toString(): String = "Try/atomic failure"
end

(* Should take a spawned thread as an argument *)
object AtomicSpawnSynchronization extends {UncheckedException}
    toString(): String = "Atomic spawn synchronization"
end

(************************************************************
* Array support
************************************************************)

oops(s:ZZ32, l:ZZ32, sz:ZZ32, got:ZZ32):() = do
   fail("Index of dimension " s " out of bounds; got " got
        " which is not in " l "#" sz);
end

trait HasRank extends Equality[\HasRank\] excludes { Number, AnyMaybe }
  (** \vspace{-4ex} NOT YET: %comprises Array[\T,E,I\]
                              where [\T,E,I\]{ T extends Array[\T,E,I\] }% *)
  rank():ZZ32
  opr =(self, other:HasRank): Boolean = false
end

(* Declared Rank-n-ness *)
trait Rank[\ nat n \] extends HasRank
  rank():ZZ32 = n
end

(* Potemkin exclusion traits.  Really we just want to say that
 * Rank[\n\] excludes Rank[\m\] where { m =/= n }, but we can't yet. *)

trait Rank1 extends { Rank[\1\]} excludes { Rank2, Rank3, Number, String }
end

trait Rank2 extends { Rank[\2\]} excludes { Rank3, Number, String }
end

trait Rank3 extends { Rank[\3\]} excludes { Number, String }
end

(* The trait Indexed_i[\n\] indicates that something has an i^th
 * dimension of size n.  In general anything which extends Indexed_i
 * must also extend Indexed_j for j < i. *)

trait Indexed1[\ nat n \] extends HasRank end

trait Indexed2[\ nat n \] extends HasRank end

trait Indexed3[\ nat n \] extends HasRank end

(** The indexed trait indicates that an object of type T can be
indexed using type I to obtain elements with type E.

An object i that's an instance of Indexed defines three basic things:
  The indexing operator opr [], which must be defined for every instance of
    the type.

  A suite of generators: i.indices generates the index space of the
    array.  i itself generates the values contained at those indices.
    i.indexValuePairs yields pairs of (index,value).  All of these
    share the same natural order.  It is necessary to define one of
    indices() and indexValuePairs(), in addition to generate() (but
    the latter requirement can be dispensed by instead extending
    DelegatedIndexed).

  A set of utility functions, assign, fill, and copy.  Only fill and
    copy need to be defined.
**)
trait Indexed[\E, I\] extends Generator[\E\]
    (** isEmpty indicates whether there are any valid indices.  It is
        defined as size=0 *)
    getter isEmpty(): Boolean = |self| = 0
    getter nonEmpty(): Boolean = NOT self.isEmpty
    (** %self.size% is equivalent to %|self|%. *)
    getter size(): ZZ32
    (** bounds() yields a range of indices that are valid for the
        indexed generator. *)
    getter bounds(): FullRange[\I,true\]
    (** indexValuePairs() generates the elements of the indexed object
        (exactly those elements that are generated by the object itself),
        but each element is paired with its index.  When we obtain
        (i,v) from indexValuePairs we know that:
           self[i] = v
           the i are distinct and i IN bounds()
           stripping away the i yields exactly the results of v <- self
        This generator attempts to follow the structure of the
        underlying object as closely as possible.  *)
    getter indexValuePairs(): Indexed[\(I,E),I\] =
        indices().map[\(I,E)\](fn (i:I): (I,E) => (i,self[i]))
    (** indices() yields the indices corresponding to the elements of
        the indexed object---it corresponds to the index component of
        indexValuePairs().  This may in general be a subset of all the
        valid indices represented by bounds().  This generator
        attempts to follow the structure of the underlying object as
        closely as possible. *)
    getter indices(): Indexed[\I,I\] =
        indexValuePairs().map[\I\](fn (i:I, e:E): I => i)

    (** %|self|% indicates the number of distinct valid indices that may
        be passed to indexing operations. *)
    opr |self| : ZZ32

    (** Indexing.  i IN bounds() must hold. *)
    opr[i:I] : E

    (** Locality of a particular index (NOT of its contents!). *)
    region(i:I): Region = Global

    (** Indexing by ranges.  The results are 0-based when the
        underlying index type has a notion of 0.  This ensures
        consistency of behavior between types such as vectors that
        *only* support 0 indexing and types such as arrays that permit
        other choices of lower bounds.  The easiest way to write the
        index by ranges operation for an instance of Indexed is to
        take advantage of indexing on the ranges themselves by writing
        (bounds())[r] in order to narrow and bounds check the range r
        and obtain a closed range of indices on the underlying
        data. **)
    opr[r:Range[\I\]] : Indexed[\E,I\]
    opr[_:OpenRange[\Any\]] : Indexed[\E,I\] = self[OpenRange[\I\]]

    (** Roughly speaking, ivmap(f) is equivalent to
        indexValuePairs.map(f).  However ivmap function isn't merely a
        convenient shortcut.  It's actually intended to create a copy
        of the underlying indexed structure when that is appropriate.

        The usual map function in Generator should do the same (and
        does for the instances in this library).  Copying can be bad
        for space, but is complexity-preserving if the mapped
        generator is used more than once. **)
    ivmap[\R\](f:(I,E)->R): Indexed[\R, I\] = indexValuePairs().map[\R\](f)
    map[\R\](f:E->R): Indexed[\R, I\] = SimpleMappedIndexed[\E,R,I\](self,f)

    (** indexOf(e) returns an index at which e can be found,
        or Nothing if no such index exists. **)
    indexOf(e:E): Maybe[\I\] =
        label found
            for (i,v) <- indexValuePairs(), v=e do
                exit found with Just[\I\](i)
            end
            Nothing[\I\]
        end found
end

trait ZeroIndexed[\E\] extends Indexed[\E,ZZ32\]
    bounds(): FullRange[\ZZ32,true\] = 0 # |self|
    zip[\F\](g:ZeroIndexed[\F\]):ZeroIndexed[\(E,F)\] =
        DefaultZip[\E,F\](self,g)
end

object DefaultZip[\E,F\](e:ZeroIndexed[\E\],f:ZeroIndexed[\F\])
        extends { ZeroIndexed[\(E,F)\], DelegatedIndexed[\(E,F),ZZ32\] }
    getter size(): ZZ32 = |self|
    getter indices(): Generator[\ZZ32\] =
        if |e| <= |f| then e.indices else f.indices end
    opr |self| : ZZ32 = |e| MIN |f|
    opr[i:ZZ32]:(E,F) = (e[i],f[i])
    opr[r:Range[\ZZ32\]] : ZeroIndexed[\(E,F)\] =
        DefaultZip[\E,F\](e[r],f[r])
end

trait LexicographicOrder[\T extends LexicographicOrder[\T,E\],E\]
        extends { StandardTotalOrder[\T\], ZeroIndexed[\E\] }
    opr CMP(self, other:T): TotalComparison =
        self.zip[\E\](other).generate[\TotalComparison\](
            LexicographicReduction,
            fn(a:E,b:E): TotalComparison => a CMP b) LEXICO
        (|self| CMP |other|)
    (** We give a specialized version of = because it can fail faster
        than CMP by checking sizes early. **)
    opr =(self,other:T): Boolean =
        (|self| = |other|) AND:
          (BIG AND [(a,b)<-self.zip[\E\](other)] a=b)
end

toArray[\E\](g:Indexed[\E,ZZ32\]): Array[\E,ZZ32\] = do
    bnds = g.bounds
    r = array[\E\](|bnds|).shift(bnds.lower())
    for (i,v) <- g.indexValuePairs do
        r.init(i,v)
    end
    r
  end

(** DelegatedIndexed is an Indexed generator that has recourse to
    another Indexed generator() internally.  By default this in turn
    is defined in terms of indexValuePairs().  Thus it's only
    necessary to define either indexValuePairs() or indices().

    This class is designed for convenience; it shouldn't be used as a
    type in runing code, but only as a supertype in lieu of Indexed.
**)
trait DelegatedIndexed[\E,I\] extends Indexed[\E,I\]
    getter generator(): Generator[\E\] =
        indexValuePairs().map[\E\](fn (i:I, e:E): E => e)
    getter size(): ZZ32 = |generator()|
    opr |self| : ZZ32 = |generator()|
    generate[\R\](r: Reduction[\R\], body: E->R): R =
        generator().generate[\R\](r,body)
    seq(self): SequentialGenerator[\E\] = seq(generator())
    cross[\G\](g: Generator[\G\]): Generator[\(E,G)\] =
        generator().cross[\G\](g)
    mapReduce[\R\](body: E->R, join:(R,R)->R, id:R): R =
        generator().mapReduce[\R\](body,join,id)
    reduce(j:(E,E)->E, z:E):E = generator().reduce(j,z)
    reduce(r: Reduction[\E\]):E = generator().reduce(r)
    loop(f:E->()): () = generator().loop(f)
end

(** The MutableIndexed trait is an indexed trait whose elements can be
    mutated using indexed assignment.  Right now we're using this type
    in a somewhat dangerous way, since eg Array1[\E,b0,s0\] extends
    both Indexed[\Array1[\E,b0,s0\],E,ZZ32\] and
    Indexed[\Array[\E,ZZ32\],E,ZZ32\].  We will need to find a
    solution to this at some point.
**)
trait MutableIndexed[\E, I\]
        extends { Indexed[\E,I\] }
    opr[i:I]:=(v:E) : ()

    (** For Ranged assignment, the extents of r and v.bounds must
        match, but the lower bounds need not. **)
    opr[r:Range[\I\]]:=(v:Indexed[\E,I\]) : ()
    opr[_:OpenRange[\Any\]]:=(v:Indexed[\E,I\]) : () =
        do self[bounds()] := v end
end

(* Array whose bounds are implicit rather than static, and which may
   be either mutable or immutable. *)
trait ReadableArray[\E,I\]
        extends { HasRank, Indexed[\E,I\], DelegatedIndexed[\E,I\] }
        comprises { Array[\E,I\], ImmutableArray[\E,I\] }
    (** CONCRETE GETTERS
        Default implementations of getters based on abstract methods
        below. **)
    getter indices(): Indexed[\I,I\] = bounds()
    getter indexValuePairs(): Indexed[\(I,E),I\] =
        zeroIndices().map[\(I,E)\](fn (i:I):(I,E) => (toIndex(i),get(i)))
    getter generator(): Indexed[\E,I\] =
        zeroIndices().map[\E\](fn (i:I):E => get(i))

    (** CONCRETE METHODS
        Default implementations of most array stuff based on the above.
        The things we can't provide are anything involving replica. **)

    opr[i:I]:E = get(offset(i))

    (** Initialize element at index i with value v.  This should occur
        once, before any other access or assignment occurs to element
        i.  An error will be signaled if an uninitialized element is
        read or an initialized element is re-initialized. **)
    init(i:I, v:E) = init0(offset(i),v)

    generate[\R\](r: Reduction[\R\], body: E->R): R =
        generator().generate[\R\](r,body)
    seq(self): SequentialGenerator[\E\] = seq(generator())

    (** 0-based non-bounds-checked indexing **)
    get(i:I): E
    init0(i:I, e:E): ()
    zeroIndices(): FullRange[\I,true\]
    (** Convert from base()-based indexing to 0-based indexing,
        performing bounds checking. **)
    offset(i:I): I
    (** Convert from 0-based indexing to base()-based indexing **)
    toIndex(i:I): I

    (** Indexed functionality with more specific type information **)
    opr[r:Range[\I\]] : ReadableArray[\E,I\]
    opr[_:OpenRange[\Any\]] : ReadableArray[\E,I\]
    ivmap[\R\](f:(I,E)->R): ReadableArray[\R, I\]
    map[\R\](f:E->R): ReadableArray[\R, I\]

    (** Shift the origin of an array.  This should yield a new view of
        the same array; ie initialization and/or update to either will
        be reflected in the other. **)
    shift(newOrigin:I):ReadableArray[\E,I\]

    (** Bulk initialization of an array using a given function or
        value.  These are defined with more specific self types in
        StandardImmutableArrayType. **)
    fill(f:I->E):ReadableArray[\E,I\]
    fill(v:E):ReadableArray[\E,I\]

    copy():ReadableArray[\E,I\]

    (** Create a fresh array structurally identical to the present
        one, but holding elements of type U. **)
    replica[\U\]():ReadableArray[\U,I\]

    opr =(self, other:HasRank): Boolean =
        typecase other of
            ReadableArray[\E,I\] => indices().generate[\Boolean\](AndReduction,
                                           fn (i:I):Boolean => (self[i]=other[i]))
            else => false
        end

end

trait ImmutableArray[\E,I\] extends { ReadableArray[\E,I\] }
        excludes { Array[\E,I\] }
    opr[r:Range[\I\]] : ImmutableArray[\E,I\]
    opr[_:OpenRange[\Any\]] : ImmutableArray[\E,I\]
    ivmap[\R\](f:(I,E)->R): ImmutableArray[\R, I\]
    map[\R\](f:E->R): ImmutableArray[\R, I\]
    shift(newOrigin:I):ImmutableArray[\E,I\]
    fill(f:I->E):ImmutableArray[\E,I\]
    fill(v:E):ImmutableArray[\E,I\]
    copy():ImmutableArray[\E,I\]
    replica[\U\]():ImmutableArray[\U,I\]

    (** Thaw array (return mutable copy) **)
    thaw():Array[\E,I\]
end

trait Array[\E,I\] extends { ReadableArray[\E,I\], MutableIndexed[\E,I\] }
    (** 0-based non-bounds-checked indexing **)
    put(i:I, e:E): ()
    opr[i:I]:=(v:E):() = put(offset(i),v)

    opr[r:Range[\I\]]:=(a:Indexed[\E,I\]):() = do
        a0 = a[a.bounds]        (* Make a have 0 origin *)
        s0 = self[r]            (* Make self have 0 origin *)
        if a0.bounds = s0.bounds then
            for (i,v) <- a0.indexValuePairs do
                s0.put(i,v)
            end
        else
            fail("Can't assign indexed with bounds " a.bounds //
                 " to differently-sized subarray " self.bounds)
        end
      end

    (** The following are repetition with better return type info. **)
    opr[r:Range[\I\]] : Array[\E,I\]
    opr[_:OpenRange[\Any\]] : Array[\E,I\]
    ivmap[\R\](f:(I,E)->R): Array[\R, I\]
    map[\R\](f:E->R): Array[\R, I\]
    shift(newOrigin:I):Array[\E,I\]
    fill(f:I->E):Array[\E,I\]
    fill(v:E):Array[\E,I\]
    assign(f:I->E):Array[\E,I\]
    copy():Array[\E,I\]
    replica[\U\]():Array[\U,I\]

    (** Freeze array (return mutable copy) **)
    freeze(): ImmutableArray[\E,I\]
end

(** Factory for arrays that returns an empty 0-indexed array of a given
    run-time-determined size. **)
array[\E\](x:ZZ32):Array[\E,ZZ32\] = __arr1(__Proxy[\E\],reflect(x))
array[\E\](x:ZZ32,y:ZZ32):Array[\E,(ZZ32,ZZ32)\] =
    __arr2(__Proxy[\E\],reflect(x), reflect(y))
array[\E\](x:ZZ32,y:ZZ32,z:ZZ32):Array[\E,(ZZ32,ZZ32,ZZ32)\] =
    __arr3(__Proxy[\E\],reflect(x), reflect(y), reflect(z))

(* This should be local to array, but we don't support local
   parametric methods in the language spec at all at the moment. *)
__arr1[\E, nat n\](w:__Proxy[\E\],x:N[\n\]):Array1[\E,0,n\] =
    array1[\E,n\]()
__arr2[\E, nat n, nat m\](w:__Proxy[\E\],x:N[\n\],y:N[\m\]):
    Array2[\E,0,n,0,m\] = array2[\E,n,m\]()
__arr3[\E, nat n, nat m, nat p\](w:__Proxy[\E\],x:N[\n\],y:N[\m\],z:N[\p\]):
    Array3[\E,0,n,0,m,0,p\] = array3[\E,n,m,p\]()

(** Factory for immutable arrays that returns an empty 0-indexed array
    of a given run-time-determined size. **)
immutableArray[\E\](x:ZZ32):ImmutableArray[\E,ZZ32\] =
    __imm1(__Proxy[\E\],reflect(x))
(*
immutableArray[\E\](x:ZZ32,y:ZZ32):ImmutableArray[\E,(ZZ32,ZZ32)\] =
    __imm2(__Proxy[\E\],reflect(x), reflect(y))
immutableArray[\E\](x:ZZ32,y:ZZ32,z:ZZ32):ImmutableArray[\E,(ZZ32,ZZ32,ZZ32)\] =
    __imm3(__Proxy[\E\],reflect(x), reflect(y), reflect(z))
*)

(* This should be local to immutableArray, but we don't support local
   parametric methods in the language spec at all at the moment. *)
__imm1[\E, nat n\](w:__Proxy[\E\],x:N[\n\]):ImmutableArray1[\E,0,n\] =
   immutableArray1[\E,n\]()
(*
__imm2[\E, nat n, nat m\](w:__Proxy[\E\],x:N[\n\],y:N[\m\]):
    ImmutableArray2[\E,0,n,0,m\] = immutableArray2[\E,n,m\]()
__imm3[\E, nat n, nat m, nat p\](w:__Proxy[\E\],x:N[\n\],y:N[\m\],z:N[\p\]):
    ImmutableArray3[\E,0,n,0,m,0,p\] = immutableArray3[\E,n,m,p\]()
*)

primitiveArray[\E\](x:ZZ32):Array[\E,ZZ32\] = __parr(__Proxy[\E\],reflect(x))
__parr[\E, nat n\](w:__Proxy[\E\],x:N[\n\]):PrimitiveArray[\E,n\] =
    PrimitiveArray[\E,n\]()

primitiveImmutableArray[\E\](x:ZZ32):ImmutableArray[\E,ZZ32\] =
    __piarr(__Proxy[\E\],reflect(x))
__piarr[\E, nat n\](w:__Proxy[\E\],x:N[\n\]):PrimImmutableArray[\E,n\] =
    PrimImmutableArray[\E,n\]()

(** NOTE: StandardImmutableArrayType is a parent of
    StandardMutableArrayType.  It therefore doesn't extend
    ImmutableArray as you might expect.  Other types that extend
    it should also extend ImmutableArray explicitly. **)
trait StandardImmutableArrayType[\T extends StandardImmutableArrayType[\T,E,I\],E,I\]
        extends { ReadableArray[\E,I\] }
    fill(f:I->E):T = do
        for i <- zeroIndices() do init0(i,f(toIndex(i))) end
        self
      end
    fill(v:E):T = do
        for i <- zeroIndices() do init0(i,v) end
        self
      end

    copy():T
end


trait StandardMutableArrayType[\T extends StandardMutableArrayType[\T,E,I\],E,I\]
        extends { StandardImmutableArrayType[\T,E,I\], Array[\E,I\] }
    assign(v:T):T = do
        for i <- zeroIndices() do put(i,v.get(i)) end
        self
      end
    assign(f:I->E):T = do
        for i <- zeroIndices() do put(i,f(toIndex(i))) end
        self
      end
end

(* Canonical partitioning of a positive number x into two pieces.  If
     (a,b) = partition(n)
   and n > 0 then 0 < a <= b,  n = a + b.
   As it turns out we choose a to be the largest power of 2 < n.
*)
partition(x:ZZ32):(ZZ32,ZZ32) = do
    m = partitionL(x)
    (m,x-m)
  end

(** A ReadableArray1[\T,b0,s0\] is an arbitrary 1-dimensional array
    whose s0 elements are of type T, and whose lowest index is b0.

    The natural order of all generators is from b0 to b0+s0-1. **)
trait ReadableArray1[\T, nat b0, nat s0\]
        extends { Indexed1[\s0\], Rank1, ReadableArray[\T,ZZ32\] }
        comprises { ImmutableArray1[\T,b0,s0\], Array1[\T,b0,s0\] }
    getter size():ZZ32 = s0
    getter bounds():ScalarRange[\ZZ32,true\] = ParRange[\ZZ32,true\](b0,s0)
    getter mutability():String
    getter toString() = do
        r = "[" b0 "#" s0 "]" mutability() " = ["
        if s0=0 then
            r "]"
        else
            f(i:ZZ32,t:T):String = " " t
            j(a:String,b:String):String = a b
            r indexValuePairs().mapReduce[\String\](f,j,"") " ]"
        end
      end

    opr |self| : ZZ32 = s0

    subarray[\nat b, nat s, nat o\]():ReadableArray1[\T, b, s\]

    (* Offset converts from b0 indexing to 0 indexing,
       bounds checking en route *)
    offset(i:ZZ32):ZZ32 = do
        r = i - b0
        if NOT (0 <= r < s0) then oops(1, b0, s0, i) end
        r
    end
    toIndex(i:ZZ32):ZZ32 = i + b0

    zeroIndices(): ScalarRange[\ZZ32,true\] = ParRange[\ZZ32,true\](0,s0)
end

trait ImmutableArray1[\T, nat b0, nat s0\]
    extends { StandardImmutableArrayType[\ImmutableArray1[\T,b0,s0\],T,ZZ32\],
              ImmutableArray[\T,ZZ32\], ReadableArray1[\T,b0,s0\] }
    getter mutability():String = "(immutable)"
    shift(o:ZZ32): ImmutableArray[\T,ZZ32\] =
        if o=0 then
            self
        else
            __subarrayI(self,reflect(o),N[\s0\],N[\b0\])
        end

    opr[r: Range[\ZZ32\]] : ImmutableArray[\T,ZZ32\] = do
        r' = (bounds())[r]
        z = N[\0\]
        s = reflect( |r'| )
        l = reflect(r'.lower())
        __subarrayI(self,z,s,l)
      end
    opr[_:OpenRange[\ZZ32\]] : ImmutableArray1[\T,0,s0\] =
        subarray[\0,s0,0\]()
    opr[_:OpenRange[\Any\]] : ImmutableArray1[\T,0,s0\] =
        subarray[\0,s0,0\]()

    (** subarray selects a subarray of this array based on static parameters.
        b#s are the new bounds of the array; o is
        the index of the subarray within the current array. **)
    subarray[\nat b, nat s, nat o\]():ImmutableArray1[\T, b, s\] = do
        boundsCheck = (bounds())[o#s]
        __ImmutableSubArray1[\T,b,s,b0,s0\](self,o-b0,1)
    end

    (* the replica method returns a replica of the array (similar layout
       etc.) but with a different element type. *)
    replica[\U\]():ImmutableArray1[\U,b0,s0\] = __immutableFactory1[\U,b0,s0\]()

    copy():ImmutableArray1[\T,b0,s0\] =
        replica[\T\]().fill(fn (i:ZZ32):T => get(i-b0))

    thaw():Array1[\T,b0,s0\] =
        __builtinFactory1[\T,b0,s0\]().fill(fn (i:ZZ32):T => get(i-b0))

    map[\R\](f:T->R): ImmutableArray1[\R,b0,s0\] =
        replica[\R\]().fill(fn (i:ZZ32):R => f(get(i-b0)))
    ivmap[\R\](f:(ZZ32,T)->R): ImmutableArray1[\R,b0,s0\] =
        replica[\R\]().fill(fn (i:ZZ32):R => f(i,get(i-b0)))
end

(** Array1[\T,b0,s0\] is a 1-dimension array whose s0 elements are of
    type T, and whose lowest index is b0. **)
trait Array1[\T, nat b0, nat s0\]
    extends { ReadableArray1[\T,b0,s0\],
              StandardMutableArrayType[\Array1[\T,b0,s0\],T,ZZ32\] }
    excludes {Number, String}

    getter mutability():String = ""

    shift(o:ZZ32): Array[\T,ZZ32\] =
        if o=0 then
            self
        else
            __subarray(self,reflect(o),N[\s0\],N[\b0\])
        end

    opr[r: Range[\ZZ32\]] : Array[\T,ZZ32\] = do
        r' = (bounds())[r]
        z = N[\0\]
        s = reflect( |r'| )
        l = reflect(r'.lower())
        __subarray(self,z,s,l)
      end
    opr[_:OpenRange[\ZZ32\]] : Array1[\T,0,s0\] =
        subarray[\0,s0,0\]()
    opr[_:OpenRange[\Any\]] : Array1[\T,0,s0\] =
        subarray[\0,s0,0\]()

    (** subarray selects a subarray of this array based on static parameters.
        b#s are the new bounds of the array; o is
        the index of the subarray within the current array. **)
    subarray[\nat b, nat s, nat o\]():Array1[\T, b, s\] = do
      boundsCheck = (bounds())[o#s]
      __SimpleSubArray1[\T,b,s,b0,s0\](self,o-b0,1)
    end

    (* the replica method returns a replica of the array (similar layout
       etc.) but with a different element type. *)
    replica[\U\]():Array1[\U,b0,s0\] = __builtinFactory1[\U,b0,s0\]()

    copy():Array1[\T,b0,s0\] =
        replica[\T\]().fill(fn (i:ZZ32):T => get(i-b0))

    freeze():ImmutableArray1[\T,b0,s0\] =
        __immutableFactory1[\T,b0,s0\]().fill(fn (i:ZZ32):T => get(i-b0))

    map[\R\](f:T->R): Array1[\R,b0,s0\] =
        replica[\R\]().fill(fn (i:ZZ32):R => f(get(i-b0)))
    ivmap[\R\](f:(ZZ32,T)->R): Array1[\R,b0,s0\] =
        replica[\R\]().fill(fn (i:ZZ32):R => f(i,get(i-b0)))
end

(* WORK AROUND absence of type inference in method calls. *)
__subarray[\T, nat b0, nat s0, nat b, nat s, nat o\]
          (it:Array1[\T,b0,s0\],
           _:N[\b\],_:N[\s\],_:N[\o\]):Array1[\T,b,s\] =
    it.subarray[\b,s,o\]()

__subarrayI[\T, nat b0, nat s0, nat b, nat s, nat o\]
          (it:ImmutableArray1[\T,b0,s0\],
           _:N[\b\],_:N[\s\],_:N[\o\]):ImmutableArray1[\T,b,s\] =
    it.subarray[\b,s,o\]()


(** A 1-D subarray of an Array1.
    b_a#s_a is the underlying array's bounds; b0#s0 is the
    result bounds.  o0 is the starting offset in the underlying array,
    defined in 0-indexed space.  m0 the stride within that array.
    Invariant: (s0-1) m0 + o0 - b_a) < s_a *)
object __SimpleSubArray1[\T, nat b0, nat s0, nat b_a, nat s_a\]
       (arr:Array1[\T,b_a,s_a\], o0:ZZ32, m0:ZZ32)
       extends Array1[\T, b0, s0\]
    index(i:ZZ32):T = i m0 + o0
    get(i:ZZ32):T = arr.get(index(i))
    put(i:ZZ32,v:T): () = arr.put(index(i),v)
    init0(i:ZZ32,v:T): () = arr.init0(index(i),v)
    subarray[\nat b, nat s, nat o\]():Array1[\T, b, s\] = do
        boundsCheck = (bounds())[o#s]
        o_n : ZZ32 = o - b0
        __SimpleSubArray1[\T, b, s, b_a, s_a\](arr, index(o_n), m0)
      end
end

(** The same as above, but immutable. **)
object __ImmutableSubArray1[\T, nat b0, nat s0, nat b_a, nat s_a\]
       (arr:ImmutableArray1[\T,b_a,s_a\], o0:ZZ32, m0:ZZ32)
       extends ImmutableArray1[\T, b0, s0\]
    index(i:ZZ32):T = i m0 + o0
    get(i:ZZ32):T = arr.get(index(i))
    put(i:ZZ32,v:T): () = arr.put(index(i),v)
    init0(i:ZZ32,v:T): () = arr.init0(index(i),v)
    subarray[\nat b, nat s, nat o\]():ImmutableArray1[\T, b, s\] = do
        boundsCheck = (bounds())[o#s]
        o_n : ZZ32 = o - b0
        __ImmutableSubArray1[\T, b, s, b_a, s_a\](arr, index(o_n), m0)
      end
end

trait AnyVector end

trait Vector[\T extends Number, nat s0\]
        extends { AnyVector, Array1[\T,0,s0\], AdditiveGroup[\Vector[\T,s0\]\] }
        excludes { AnyMultiplicativeRing }
    opr +(self, v:Vector[\T,s0\]): Vector[\T,s0\] =
        ivmap[\T\](fn (i:ZZ32, e: T):T => e + v.get(i))
    opr -(self, v:Vector[\T,s0\]): Vector[\T,s0\] =
        ivmap[\T\](fn (i:ZZ32, e: T):T => e - v.get(i))
    opr -(self): Vector[\T,s0\] = map[\T\](fn (e: T):T => - e)
    scale(t: T): Vector[\T,s0\] = map[\T\](fn (v) => t v)
    pmul(v: Vector[\T,s0\]): Vector[\T,s0\] =
        ivmap[\T\](fn (i:ZZ32, e: T):T => e v.get(i))
    dot(v: Vector[\T,s0\]): T =
        SUM [(i,me_i)<-indexValuePairs()] me_i v.get(i)
end

object __DefaultVector[\T, nat s0\]() extends Vector[\T,s0\]
  mem: PrimitiveArray[\T,s0\] = PrimitiveArray[\T,s0\]()
  get(i:ZZ32):T = mem.get(i)
  put(i:ZZ32, v:T) = mem.put(i,v)
  init0(i:ZZ32, v:T) = mem.init0(i,v)
  replica[\U\]() = array1[\U,s0\]()
end

(* builtinFactory1 must be a non-overloaded 0-parameter factory for
   1-D arrays.  The type parameters are enshrined in LHSEvaluator.java
   and NonPrimitive.java; the factory name is enshrined in
   WellKnownNames.java.  There must be some factory, named in this
   file, with this type signature.  A similar thing is true for
   K-dimensional array types. *)
__builtinFactory1[\T, nat b0, nat s0\]():Array1[\T,b0,s0\] = do
    r = array1[\T,s0\]()
    typecase _ = N[\b0\] of
        N[\0\] => r
        else => r.subarray[\b0,s0,0\]()
    end
  end

(* immutableFactory1 is a non-overloaded 0-parameter factory for
   0-indexed 1-D arrays.  It is also mentioned in WellKnownNames as it
   is used to allocate storage for varargs. *)
__immutableFactory1[\T, nat b0, nat s0\]():ReadableArray1[\T,b0,s0\] = do
    r = immutableArray1[\T,s0\]()
    typecase _ = N[\b0\] of
        N[\0\] => r
        else => r.subarray[\b0,s0,0\]()
    end
  end

(* TODO: fix when Number is covariant. *)
array1[\T, nat s0\]():Array1[\T,0,s0\] =
    typecase _ = __Proxy[\T\] of
        __Proxy[\ZZ32\] => vector[\T,s0\]()
        __Proxy[\ZZ64\] => vector[\T,s0\]()
        __Proxy[\RR64\] => vector[\T,s0\]()
        __Proxy[\Number\] => vector[\T,s0\]()
        else => PrimitiveArray[\T,s0\]()
    end
array1[\T, nat s0\](v:T):Array1[\T,0,s0\] = array1[\T,s0\]().fill(v)
array1[\T, nat s0\](f:ZZ32->T):Array1[\T,0,s0\] = array1[\T,s0\]().fill(f)

immutableArray1[\T, nat s0\](): ImmutableArray1[\T,0,s0\] =
    PrimImmutableArray[\T,s0\]()

(* vector is the same as array1, but specialized to numeric type arguments *)
vector[\T extends Number, nat s0\]():Vector[\T,s0\] = __DefaultVector[\T,s0\]()
vector[\T extends Number, nat s0\](v:T):Vector[\T,s0\] =
    vector[\T,s0\]().fill(v)
vector[\T extends Number, nat s0\](f:ZZ32->T):Vector[\T,s0\] =
    vector[\T,s0\]().fill(f)

pmul[\ T extends Number, nat k \]
    (a : Vector[\T,k\], b : Vector[\T,k\]):Vector[\T,k\] = a.pmul(b)

opr DOT[\ T extends Number, nat n \]
       (me : Vector[\T,n\], other : Vector[\T,n\]):T = me.dot(other)

opr juxtaposition[\ T extends Number, nat n \]
     (me : Vector[\T,n\], other : Vector[\T,n\]):T = me.dot(other)

opr DOT[\ T extends Number, nat n \]
       (me : Vector[\T,n\], other : T) : Vector[\T,n\] = me.scale(other)

opr juxtaposition[\ T extends Number, nat n \]
     (me : Vector[\T,n\], other : T) : Vector[\T,n\] = me.scale(other)

opr DOT[\ T extends Number, nat n \]
        (other : T, me : Vector[\T,n\]) : Vector[\T,n\] = me.scale(other)

opr juxtaposition[\ T extends Number, nat n \]
     (other : T, me : Vector[\T,n\]) : Vector[\T,n\] = me.scale(other)

squaredNorm[\T extends Number, nat s0\](a:Vector[\T,s0\]):T = a.dot(a)

opr ||[\ T extends Number, nat k \]me : Vector[\T,k\]|| : RR64 = SQRT squaredNorm(me)

(** Array2[\T,b0,s0,b1,s1\] is the type of 2-dimensional arrays of
    element type T, with size s0 in the first dimension and s1 in the
    second dimension and lowest index (b0,b1).  Natural order for all
    generators in each dimension is from b to b+s-1; the overall order
    of elements need only be consistent with the cross product of
    these orderings (see Generator.cross()). **)
trait Array2[\T, nat b0, nat s0, nat b1, nat s1\]
    extends { Indexed1[\s0\], Indexed2[\s1\], Rank2,
              StandardMutableArrayType[\Array2[\T,b0,s0,b1,s1\],T,(ZZ32,ZZ32)\] }
    excludes { Number, String }
  getter size():ZZ32 = s0 s1
  getter bounds():Tuple2Range[\ZZ32,ZZ32,true\] = Tuple2Range[\ZZ32,ZZ32,true\](b0,b1,s0,s1)
  getter toString() = do
    r : String := "[" b0 "#" s0 "," b1 "#" s1 "] ="
    row(i) =
      for j <- seq(0#s1) do
         r := r " " get(i,j)
      end
    if s0 = 0 then
      r " []"
    else
      r := r "\n["
      row(0)
      for i <- seq(1#(s0-1)) do
        r := r "\n "
        row(i)
      end
      r " ]"
    end
  end
  opr |self| : ZZ32 = s0 s1
  (* Translate from b0,b1-indexing to 0-indexing, checking bounds. *)
  offset(t_1:(ZZ32,ZZ32)):(ZZ32,ZZ32) = do
    (a0,a1) = t_1
    c0 = a0 - b0
    c1 = a1 - b1
    if NOT (0 <= c0 < s0) then
        oops(1,b0,s0,a0)
    elif NOT (0 <= c1 < s1) then
        oops(2,b1,s1,a1)
    else
        (c0,c1)
    end
  end
  toIndex(t_1:(ZZ32,ZZ32)):(ZZ32,ZZ32) = do (a0,a1)=t_1; (a0+b0,a1+b1) end
  opr[x:ZZ32,y:ZZ32]:=(v:T):() = do self[ (x,y) ] := v end
  opr[r:Range[\(ZZ32,ZZ32)\]]: Array[\T,(ZZ32,ZZ32)\] = do
      r' = (bounds())[r]
      (z0,z1) = r'.extent()
      (l0,l1) = r'.lower()
      __subarray(self, N[\0\],reflect(z0), N[\0\], reflect(z1),
                 reflect(l0),reflect(l1))
    end
  opr[_:OpenRange[\ZZ32\]] : Array2[\T,0,s0,0,s1\] =
      subarray[\0,s0,0,s1,0,0\]()
  opr[_:OpenRange[\Any\]] : Array2[\T,0,s0,0,s1\] =
      subarray[\0,s0,0,s1,0,0\]()
  shift(t_1:(ZZ32,ZZ32)): Array[\T,(ZZ32,ZZ32)\] = do
      (o0,o1)=t_1
      if o0=0 AND o1=0 then
          self
      else
          __subarray(self,reflect(o0),N[\s0\], reflect(o1),N[\s1\], N[\b0\],N[\b1\])
      end
    end

  (** 2-D subarray given static subarray parameters.
      (bo1,bo2)#(so1,so2) are output bounds.
      The result is the subarray starting at (o0,o1) in the original array.
   **)
  subarray[\nat bo0, nat so0, nat bo1, nat so1, nat o0, nat o1\]
          (): Array2[\T,bo0,so0,bo1,so1\] = do
      boundsCheck = (bounds())[ (o0,o1)#(so0,so1) ]
      SubArray2[\T,bo0,so0,bo1,so1,b0,s0,b1,s1\](self,1,o0,1,o1)
    end

  zeroIndices():Tuple2Range[\ZZ32,ZZ32,true\] = Tuple2Range[\ZZ32,ZZ32,true\](0,0,s0,s1)

  replica[\U\]():Array2[\U,b0,s0,b1,s1\] =
      __builtinFactory2[\U,b0,s0,b1,s1\]()
  copy():Array2[\T,b0,s0,b1,s1\] =
      self.replica[\T\]().fill(fn (i:ZZ32,j:ZZ32):T => get(i-b0,j-b1))
  put(t1:(ZZ32, ZZ32), v:T) : ()
  get(t1:(ZZ32, ZZ32)):T
  t():Array2[\T,b1,s1,b0,s0\] = TransposedArray2[\T,b1,s1,b0,s0\](self)
  (* Copied here for better return type information. *)
  map[\R\](f:T->R): Array2[\R,b0,s0,b1,s1\] =
      replica[\R\]().fill(fn (i:ZZ32,j:ZZ32):R => f(get(i-b0,j-b1)))
  ivmap[\R\](f:((ZZ32,ZZ32),T)->R): Array2[\R,b0,s0,b1,s1\] =
      replica[\R\]().fill(fn (i:ZZ32,j:ZZ32):R => f((i,j),get(i-b0,j-b1)))

  freeze():ImmutableArray[\T,(ZZ32,ZZ32)\] = fail("Freeze not defined yet!")
end

(** Hoisted out here to work around the lack of method type inference. **)
__subarray[\T, nat b0, nat s0, nat b1, nat s1,
           nat bo0, nat so0, nat bo1, nat so1, nat o0, nat o1\]
        (this:Array2[\T,b0,s0,b1,s1\],
         _:N[\bo0\],_:N[\so0\],_:N[\bo1\],_:N[\so1\],
         _:N[\o0\],_:N[\o1\]): Array2[\T,bo0,so0,bo1,so1\] = do
(*
    println("__subarray[\\" b0 "," s0 "," b1 "," s1 ", "
            bo0 "," so0 "," bo1 "," so1 ", " o0 "," o1 "\\]()")
*)
    this.subarray[\bo0,so0,bo1,so1,o0,o1\]()
  end

(** Default array is column-major, but we could switch it. **)
object __DefaultArray2[\T, nat b0, nat s0, nat b1, nat s1\]()
        extends Array2[\T, b0, s0, b1, s1\]
    mem:PrimitiveArray[\T, (s0 s1) \] = PrimitiveArray[\T, (s0 s1) \]()
    init0(t_1:(ZZ32,ZZ32), v:T) : () = do (i,j)=t_1; mem.init0(i s1 + j, v) end
    put(t_1:(ZZ32,ZZ32), v:T) : () = do (i,j)=t_1; mem.put(i s1 + j, v) end
    get(t_1:(ZZ32,ZZ32)):T = do (i,j)=t_1; mem.get(i s1 + j) end
end

(** Transposes the index of the underlying array mem, and transposes
    its natural order as well. **)
object TransposedArray2[\T, nat b0, nat s0, nat b1, nat s1\]
                       (mem:Array2[\T,b1,s1,b0,s0\])
    extends Array2[\T, b0, s0, b1, s1\]
  replica[\U\]():Array2[\U,b0,s0,b1,s1\] = mem.replica[\U\]().t()

  init0(t_1:(ZZ32,ZZ32), v:T) : () = do (i,j)=t_1; mem.init0((j,i),v) end
  put(t_1:(ZZ32,ZZ32), v:T) : () = do (i,j)=t_1; mem.put((j,i),v) end
  get(t_1:(ZZ32,ZZ32)):T = do (i,j)=t_1; mem.get(j,i) end
  t() = mem
end

(** Simple 2-D subarray, used when we have no other knowledge of the
    underlying structure.  If we're on top of a PrimitiveArray we can
    really do much better than this.
    (b0,b1)#(s0,s1) are the bounds of the subarray.
    (bu0,bu1)#(su0,su1) are the bounds of the underlying array.
    m_i, o_i are multiplier and offset (0-based) of dimension i.
**)
object SubArray2[\T, nat b0, nat s0, nat b1, nat s1,
                     nat bu0, nat su0, nat bu1, nat su1\]
                (mem: Array2[\T,bu0,su0,bu1,su1\],
                 m0:ZZ32, o0:ZZ32, m1:ZZ32, o1:ZZ32)
        extends Array2[\T,b0,s0,b1,s1\]
    index(a0:ZZ32, a1:ZZ32): (ZZ32,ZZ32) = (m0 a0 + o0, m1 a1 + o1)
    init0(t_1:(ZZ32,ZZ32), v:T): () = mem.init0(index(t_1),v)
    put(t_1:(ZZ32,ZZ32), v:T): () = mem.put(index(t_1),v)
    get(t_1:(ZZ32,ZZ32)): T = mem.get(index(t_1))
end

trait AnyMatrix end

trait Matrix[\T extends Number, nat s0, nat s1\]
        extends { AnyMatrix, Array2[\T, 0, s0, 0, s1\], AdditiveGroup[\Matrix[\T,s0,s1\]\] }
        excludes { AnyMultiplicativeRing }
    opr +(self, v:Matrix[\T,s0,s1\]): Matrix[\T,s0,s1\] =
        ivmap[\T\](fn (i:(ZZ32,ZZ32),e:T):T => e + v.get(i))
    opr -(self, v:Matrix[\T,s0,s1\]): Matrix[\T,s0,s1\] =
        ivmap[\T\](fn (i:(ZZ32,ZZ32),e:T):T => e - v.get(i))
    opr -(self): Matrix[\T,s0,s1\] = map[\T\](fn (e:T):T => - e)
    scale(t_1: T): Matrix[\T,s0,s1\] = map[\T\](fn (e:T):T => t_1 e)
    mul[\ nat s2 \](other: Matrix[\T,s1,s2\]): Matrix[\T,s0,s2\] = do
        res = matrix[\T,s0,s2\]()
        mma(a:ZZ32,i:ZZ32,b:ZZ32,j:ZZ32,c:ZZ32,k:ZZ32):() =
            if k>=i AND k>=j then
              if k=1 then
                pr : T = get(a,b) other.get(b,c)
                (* If this were atomic, we could parallelize j-partition. *)
                res.put((a,c), res.get(a,c) + pr)
              else
                (k0,k1) = partition(k)
                (mma(a,i,b,j,c,k0),mma(a,i,b,j,c+k0,k1))
              end
            elif j>=i then
                (j0,j1) = partition(j)
                mma(a,i,b,j0,c,k)
                mma(a,i,b+j0,j1,c,k)
            else
                (i0,i1) = partition(i)
                (mma(a,i0,b,j,c,k),mma(a+i0,i1,b,j,c,k))
            end
        mm(a:ZZ32,i:ZZ32,b:ZZ32,j:ZZ32,c:ZZ32,k:ZZ32):() =
            if k>=i AND k>=j then
              if k=1 then
                res.put((a,c), get(a,b) other.get(b,c))
              else
                (k0,k1) = partition(k)
                (mm(a,i,b,j,c,k0),mm(a,i,b,j,c+k0,k1))
              end
            elif j>=i then
                (j0,j1) = partition(j)
                mm(a,i,b,j0,c,k)
                mma(a,i,b+j0,j1,c,k)
            else
                (i0,i1) = partition(i)
                (mm(a,i0,b,j,c,k),mm(a+i0,i1,b,j,c,k))
            end
        if s0=0 OR s1=0 OR s2=0 then
          res
        else
          mm(0,s0,0,s1,0,s2)
          res
        end
      end
    rmul(v: Vector[\T,s1\]): Vector[\T,s0\] = do
        row(i:ZZ32):T =
            SUM[(j,v_j)<-v.indexValuePairs] get(i,j) v_j
        vector[\T,s0\]().fill(row)
      end
    lmul(v: Vector[\T,s0\]): Vector[\T,s1\] = do
        col(i:ZZ32):T =
            SUM[(j,v_j)<-v.indexValuePairs] v_j get(j,i)
        vector[\T,s1\]().fill(col)
      end
    t(): Matrix[\T,s1,s0\] = TransposedMatrix[\T,s1,s0\](self)
end

(* Default matrix should match default array in column vs row-major. *)
object __DefaultMatrix[\T, nat s0, nat s1\]()
    extends Matrix[\T, s0, s1\]
  mem:PrimitiveArray[\T, (s0 s1) \] = PrimitiveArray[\T, (s0 s1) \]()
  init0(t_1:(ZZ32,ZZ32), v:T) : () = do (i,j)=t_1; mem.init0(i s1 + j, v) end
  put(t_1:(ZZ32,ZZ32), v:T) : () = do (i,j)=t_1; mem.put(i s1 + j, v) end
  get(t_1:(ZZ32,ZZ32)):T = do (i,j)=t_1; mem.get(i s1 + j) end
end

object TransposedMatrix[\T, nat s0, nat s1\](mem:Matrix[\T,s1,s0\])
        extends Matrix[\T, s0, s1\]
    replica[\U\]():Array2[\U,0,s0,0,s1\] = mem.replica[\U\]().t()

    init0(t_1:(ZZ32,ZZ32), v:T) : () = do (i,j)=t_1; mem.init0((j,i),v) end
    put(t_1:(ZZ32,ZZ32), v:T) : () = do (i,j)=t_1; mem.put((j,i),v) end
    get(t_1:(ZZ32,ZZ32)):T = do (i,j)=t_1; mem.get(j,i) end
    t(): Matrix[\T,s1,s0\] = mem
    add(v:TransposedMatrix[\T,s0,s1\]): Matrix[\T,s0,s1\] = mem.add(v.t()).t()
    subtract(v:Matrix[\T,s0,s1\]): Matrix[\T,s0,s1\] = mem.subtract(v.t()).t()
    negate(): Matrix[\T,s0,s1\] = mem.negate().t()
    scale(f: T): Matrix[\T,s0,s1\] = mem.scale(f).t()
    rmul(v: Vector[\T,s1\]): Vector[\T,s0\] = mem.lmul(v)
    lmul(v: Vector[\T,s0\]): Vector[\T,s1\] = mem.rmul(v)
(*  Can't overload generic methods yet, but this is preferable.
    mul[\nat s2\](v:TransposedMatrix[\T,s1,s2\]): Matrix[\T,s0,s2\] =
        v.t().mul(mem).t()
*)
end

__builtinFactory2[\T,nat b0,nat s0,nat b1,nat s1\]():Array2[\T,b0,s0,b1,s1\] =
    if b0=0 AND b1=0 then
        array2[\T,s0,s1\]()
    else
        __DefaultArray2[\T,b0,s0,b1,s1\]()
    end

(* array2 is a factory for 0-based 2-D arrays. *)
(* TODO: fix when Number is covariant. *)
array2[\T, nat s0, nat s1\]():Array2[\T,0,s0,0,s1\] =
    typecase _ = __Proxy[\T\] of
        __Proxy[\ZZ32\] => matrix[\T,s0,s1\]()
        __Proxy[\ZZ64\] => matrix[\T,s0,s1\]()
        __Proxy[\RR64\] => matrix[\T,s0,s1\]()
        __Proxy[\Number\] => matrix[\T,s0,s1\]()
        else => __DefaultArray2[\T,0,s0,0,s1\]()
    end
array2[\T, nat s0, nat s1\](v:T):Array2[\T,0,s0,0,s1\] =
  array2[\T,s0,s1\]().fill(v)
array2[\T, nat s0, nat s1\](f:(ZZ32,ZZ32)->T):Array2[\T,0,s0,0,s1\] =
  array2[\T,s0,s1\]().fill(f)

(* matrix is the same as array1, but specialized to numeric type
   arguments, except that the default value (if given) is used to
   construct a multiple of the identity matrix. *)
matrix[\T extends Number, nat s0, nat s1\]():Matrix[\T,s0,s1\] =
  __DefaultMatrix[\T,s0,s1\]()
matrix[\T extends Number, nat s0, nat s1\](v:T):Matrix[\T,s0,s1\] =
  array2[\T,s0,s1\]().fill(fn (x:ZZ32,y:ZZ32):T => if x=y then v else 0 end)

(* Matrix multiplication; used to use a cache-oblivious algorithm, but
   we ran into trouble due to lack of support for atomic increment of
   matrix elements. *)
opr DOT[\ T extends Number, nat n, nat m, nat p\]
       (me:Matrix[\T,n,m\], other:Matrix[\T,m,p\]): Matrix[\T,n,p\] =
        me.mul[\p\](other)

opr juxtaposition[\ T extends Number, nat n, nat m, nat p\]
     (me:Matrix[\T,n,m\], other:Matrix[\T,m,p\]): Matrix[\T,n,p\] =
        me.mul[\p\](other)


(* matrix-vector multiplication *)
opr DOT[\ T extends Number, nat n, nat m \]
       (me:Matrix[\T,n,m\], v:Vector[\T,m\]):Vector[\T,n\] = me.rmul(v)

opr juxtaposition[\ T extends Number, nat n, nat m \]
     (me:Matrix[\T,n,m\], v:Vector[\T,m\]):Vector[\T,n\] = me.rmul(v)

(* vector-matrix multiplication *)
opr DOT[\ T extends Number, nat n, nat m \]
       (v:Vector[\T,n\], me:Matrix[\T,n,m\]):Vector[\T,m\] = me.lmul(v)

opr juxtaposition[\ T extends Number, nat n, nat m \]
     (v:Vector[\T,n\], me:Matrix[\T,n,m\]):Vector[\T,m\] = me.lmul(v)

opr DOT[\ T extends Number, nat n, nat m \]
       (me : Matrix[\T,n,m\], other : T) : Matrix[\T,n,m\] = me.scale(other)

opr juxtaposition[\ T extends Number, nat n, nat m \]
     (me : Matrix[\T,n,m\], other : T) : Matrix[\T,n,m\] = me.scale(other)

opr DOT[\ T extends Number, nat n, nat m \]
       (other : T, me : Matrix[\T,n,m\]) : Matrix[\T,n,m\] = me.scale(other)

opr juxtaposition[\ T extends Number, nat n, nat m \]
     (other : T, me : Matrix[\T,n,m\]) : Matrix[\T,n,m\] = me.scale(other)

(** Array3[\T,b0,s0,b1,s1,b2,s2\] is the type of 3-dimensional arrays
    of element type T, with size s_i in the i^th dimension and lowest
    index (b0,b1,b2).  Natural order for all generators in each
    dimension is from b to b+s-1; the overall order of elements need
    only be consistent with the cross product of these orderings (see
    Generator.cross()). **)
trait Array3[\T, nat b0, nat s0, nat b1, nat s1, nat b2, nat s2\]
    extends { Indexed1[\s0\], Indexed2[\s1\], Indexed3[\s2\], Rank3,
              StandardMutableArrayType[\Array3[\T,b0,s0,b1,s1,b2,s2\],T,
                                        (ZZ32,ZZ32,ZZ32)\] }
    excludes { Number, String }

    getter size():ZZ32 = s0 s1 s2
    getter bounds():Tuple3Range[\ZZ32,ZZ32,ZZ32,true\] =
        Tuple3Range[\ZZ32,ZZ32,ZZ32,true\](b0,b1,b2,s0,s1,s2)

    getter toString():String = do
      r : String := "[" b0 "#" s0 "," b1 "#" s1 "," b2 "#" s2 "] ="
      row(i,k) =
          for j <- seq(0#s1) do
            r := r " " self[b0+i,b1+j,b2+k]
          end
      plane(k) =
        if s1 > 0 then
          row(0,k)
          for i <- seq(1#(s0-1)) do
            r := r "\n "
            row(i,k)
          end
        end
      if s0=0 then
        r " []"
      else
        r := r "\n["
        plane(0)
        for k <- seq(1#(s2-1)) do
          r := r " ;;\n "
          plane(k)
        end
        r " ]"
      end
    end

    opr |self| : ZZ32 = s0 s1 s2

    (* Again, offset performs bounds checking and shifts to 0 indexing. *)
    offset(t:(ZZ32,ZZ32,ZZ32)):(ZZ32,ZZ32,ZZ32) = do
      (a0,a1,a2)=t
      c0 = a0 - b0; c1 = a1 - b1; c2 = a2 - b2
      if NOT (0 <= c0 < s0) then oops(1, b0, s0, a0); end
      if NOT (0 <= c1 < s1) then oops(2, b1, s1, a1); end
      if NOT (0 <= c2 < s2) then oops(3, b2, s2, a2); end
      (c0,c1,c2)
    end
    toIndex(t:(ZZ32,ZZ32,ZZ32)):(ZZ32,ZZ32,ZZ32) =
        do (a0,a1,a2) = t; (a0+b0,a1+b1,a2+b2) end

    (* And get and put are 0-indexed without bounds checks. *)
    put(t:(ZZ32,ZZ32,ZZ32), v:T) : ()
    get(t:(ZZ32,ZZ32,ZZ32)):T

    opr[i:ZZ32, j:ZZ32, k:ZZ32] := (v:T) = do self[ (i,j,k) ] := v end
    opr[r:Range[\(ZZ32,ZZ32,ZZ32)\]]: Array[\T,(ZZ32,ZZ32,ZZ32)\] = do
        r' = (bounds())[r]
        (z0,z1,z2) = r'.extent()
        (l0,l1,l2) = r'.lower()
        __subarray(self, N[\0\],reflect(z0), N[\0\],reflect(z1),
                   N[\0\],reflect(z2), reflect(l0),reflect(l1), reflect(l2))
      end
    opr[_:OpenRange[\ZZ32\]] : Array3[\T,0,s0,0,s1,0,s2\] =
        subarray[\0,s0,0,s1,0,s2,0,0,0\]()
    opr[_:OpenRange[\Any\]] : Array3[\T,0,s0,0,s1,0,s2\] =
        subarray[\0,s0,0,s1,0,s2,0,0,0\]()
    shift(t:(ZZ32,ZZ32,ZZ32)): Array[\T,(ZZ32,ZZ32)\] = do
        (o0,o1,o2)=t
        if o0=0 AND o1=0 AND o2=0 then
            self
        else
            __subarray(self,
                       reflect(o0),N[\s0\], reflect(o1),N[\s1\],
                       reflect(o2),N[\s2\], N[\b0\],N[\b1\],N[\b2\])
        end
      end

    (** 3-D subarray given static subarray parameters.
        (bo1,bo2)#(so1,so2) are output bounds.
        The result is the subarray starting at (o0,o1) in the original array.
     **)
    subarray[\nat bo0, nat so0, nat bo1, nat so1, nat bo2, nat so2,
              nat o0, nat o1, nat o2\]
            (): Array3[\T,bo0,so0,bo1,so1,bo2,so2\] = do
        boundsCheck = (bounds())[ (o0,o1,o2)#(so0,so1,so2) ]
        SubArray3[\T,bo0,so0,bo1,so1,bo2,so2,b0,s0,b1,s1,b2,s2\](self,1,o0,1,o1,1,o2)
      end

    zeroIndices():Tuple3Range[\ZZ32,ZZ32,ZZ32,true\] =
        Tuple3Range[\ZZ32,ZZ32,ZZ32,true\](0,0,0,s0,s1,s2)

    replica[\U\]():Array3[\U,b0,s0,b1,s1,b2,s2\] =
        __DefaultArray3[\U,b0,s0,b1,s1,b2,s2\]()
    copy():Array3[\T,b0,s0,b1,s1,b2,s2\] =
        self.replica[\T\]().fill(fn (i:ZZ32,j:ZZ32,k:ZZ32):T =>
                                     get(i-b0,j-b1,k-b2))
    map[\R\](f:T->R): Array3[\R,b0,s0,b1,s1,b2,s2\] =
        replica[\R\]().fill(fn (i:ZZ32,j:ZZ32,k:ZZ32):R =>
                                f(get(i-b0,j-b1,k-b2)))
    ivmap[\R\](f:((ZZ32,ZZ32,ZZ32),T)->R): Array3[\R,b0,s0,b1,s1,b2,s2\] =
        replica[\R\]().fill(fn (i:ZZ32,j:ZZ32,k:ZZ32):R =>
                                f((i,j,k),get(i-b0,j-b1,k-b2)))

    freeze():ImmutableArray[\T,(ZZ32,ZZ32,ZZ32)\] =
        fail("Freeze not defined yet!")
end

(** hoisted out here to work around absence of type inference for methods. **)
__subarray[\T, nat b0, nat s0, nat b1, nat s1, nat b2, nat s2,
            nat bo0, nat so0, nat bo1, nat so1, nat bo2, nat so2,
            nat o0, nat o1, nat o2\]
      (this:Array3[\T,b0,s0,b1,s1,b2,s2\],
       _:N[\bo0\],_:N[\so0\],_:N[\bo1\],_:N[\so1\],_:N[\bo2\],_:N[\so2\],
       _:N[\o0\],_:N[\o1\],_:N[\o2\]): Array3[\T,bo0,so0,bo1,so1,bo2,so2\] =
    this.subarray[\bo0,so0,bo1,so1,bo2,so2,o0,o1,o2\]()

object __DefaultArray3[\T, nat b0, nat s0, nat b1, nat s1, nat b2, nat s2\]() extends
                                        Array3[\T, b0, s0, b1, s1, b2, s2\]
  mem:PrimitiveArray[\T,(s0 (s1 s2))\] = PrimitiveArray[\T,(s0 (s1 s2))\]()

  ofs(i:ZZ32,j:ZZ32,k:ZZ32):ZZ32 = (i s1 + j) s2 + k

  init0(t:(ZZ32,ZZ32,ZZ32), v:T) : () = mem.init0(ofs(t),v)
  put(t:(ZZ32,ZZ32,ZZ32), v:T) : () = mem.put(ofs(t),v)
  get(t:(ZZ32,ZZ32,ZZ32)) : T = mem.get(ofs(t))
end

(** Simple 3-D subarray, used when we have no other knowledge of the
    underlying structure.  If we're on top of a PrimitiveArray we can
    really do much better than this.
    (b0,b1,b2)#(s0,s1,s2) are the bounds of the subarray.
    (bu0,bu1,bu2)#(su0,su1,su2) are the bounds of the underlying array.
    m_i, o_i are multiplier and offset (0-based) of dimension i.
**)
object SubArray3[\T, nat b0, nat s0, nat b1, nat s1, nat b2, nat s2,
                     nat bu0, nat su0, nat bu1, nat su1, nat bu2, nat su2\]
                (mem: Array3[\T,bu0,su0,bu1,su1,bu2,su2\],
                 m0:ZZ32,o0:ZZ32, m1:ZZ32,o1:ZZ32, m2:ZZ32, o2:ZZ32)
        extends Array3[\T,b0,s0,b1,s1,b2,s2\]
    index(a0:ZZ32, a1:ZZ32, a2:ZZ32): (ZZ32,ZZ32,ZZ32) =
        (m0 a0 + o0, m1 a1 + o1, m2 a2 + o2)
    init0(t:(ZZ32,ZZ32,ZZ32), v:T): () = mem.init0(index(t),v)
    put(t:(ZZ32,ZZ32,ZZ32), v:T): () = mem.put(index(t),v)
    get(t:(ZZ32,ZZ32,ZZ32)): T = mem.get(index(t))
end

__builtinFactory3[\T, nat b0, nat s0, nat b1, nat s1, nat b2, nat s2\]():
        Array3[\T,b0,s0,b1,s1,b2,s2\] =
  __DefaultArray3[\T,b0,s0,b1,s1,b2,s2\]()

array3[\T,nat s0, nat s1, nat s2\]():Array3[\T,0,s0,0,s1,0,s2\] =
  __DefaultArray3[\T,0,s0,0,s1,0,s2\]()

(*************

trait Monoid[\ T, opr OPLUS \]
  where { T extends Monoid[\ T, OPLUS \] }
    zero() : T
    opr OPLUS(self, other:T):T
end

*************)

(************************************************************
* Reductions
************************************************************)

trait Reduction[\L\]
    empty(): L
    join(a: L, b: L): L
end

(** Invariants:
    join must be associative with identity empty
    unlift(lift(x)) = x
 **)
trait ActualReduction[\R,L\] extends Reduction[\L\]
    getter toString():String
    lift(r: Any): L
    unlift(l:L): R
    (** If this reduction left-distributes over r, return a pair of
        reductions with the same lift and unlift **)
    leftDistribute(r: Reduction[\R\]): SomeReductionPair[\R\] =
        distribute(r)
    (** If this reduction right-distributes over r, return a pair of
        reductions with the same lift and unlift **)
    rightDistribute(r: Reduction[\R\]): SomeReductionPair[\R\] =
        distribute(r)
    (** If this reduction distributes over r, return a pair of
        reductions with the same lift and unlift **)
    distribute[\K\](r: ActualReduction[\R,K\]): SomeReductionPair[\R\] =
        NoReductionPair[\R\]
end

trait PossibleReductionPair[\R\] extends Condition[\SomeReductionPair[\R\]\]
    comprises { NoReductionPair[\R\], SomeReductionPair[\R\] }
end

object NoReductionPair[\R\] extends PossibleReductionPair[\R\]
    getter holds(): Boolean = false
    getter cond[\G\](t:PossibleReductionPair[\R\]->G, e:()->G): G = e()
end

trait SomeReductionPair[\R\] extends PossibleReductionPair[\R\]
    getter holds(): Boolean = true
    getter cond[\G\](t:PossibleReductionPair[\R\]->G, e:()->G): G = t(self)
    getter outer(): Reduction[\R\]
    getter inner(): Reduction[\R\]
end

trait ReductionPair[\R,L\] extends SomeReductionPair[\R\]
    getter outer(): ActualReduction[\R,L\]
    getter inner(): ActualReduction[\R,L\]
end

(** The usual lifting to Maybe for identity-less operators **)
trait AssociativeReduction[\R\] extends ActualReduction[\R,AnyMaybe\]
    empty(): Nothing[\R\] = Nothing[\R\]
    join(a: AnyMaybe, b: AnyMaybe): AnyMaybe =
        if av <- a then
            if bv <- b then
                Just(simpleJoin(av,bv))
            else
                a
            end
        else
            b
        end
    simpleJoin(a:Any, b:Any): Any
    lift(r:Any): AnyMaybe = Just(r)
    unlift(r:AnyMaybe): R = r.get
end

trait CommutativeReduction[\R\] extends AssociativeReduction[\R\] end

(** Monoids don't require a special lift and unlift operation. **)
trait MonoidReduction[\R\] extends ActualReduction[\R,R\]
    lift(r:Any): R = r
    unlift(r:R): R = r
end

trait CommutativeMonoidReduction[\R\] extends MonoidReduction[\R\] end

trait ReductionWithZeroes[\R,L\] extends ActualReduction[\R,L\]
    isLeftZero(l:L): Boolean = isZero(l)
    isRightZero(l:L): Boolean = isZero(l)
    isZero(l:L): Boolean = false
end

trait BigOperator[\I,O,R,L\]
    getter reduction(): ActualReduction[\R,L\]
    getter body(): I->R
    getter unwrap(): R->O
end

object BigReduction[\R,L\](r:ActualReduction[\R,L\]) extends BigOperator[\R,R,R,L\]
    getter reduction(): ActualReduction[\R,L\] = r
    getter body(): R->R = fn x => x
    getter unwrap(): R->R = fn x => x
end

object Comprehension[\I,O,R,L\](u: R->O, r: ActualReduction[\R,L\], singleton:I->R)
        extends BigOperator[\I,O,R,L\]
    getter reduction(): ActualReduction[\R,L\] = r
    getter body(): I->R = singleton
    getter unwrap(): R->O = u
end

(** VoidReduction is usually done for effect, so we pretend that
    the completion performs the effects.  This rules out things
    distributing over void (that would change the number of effects in
    our program) but not void distributing over other things. **)
object VoidReduction extends { CommutativeMonoidReduction[\()\] }
    getter toString() = "VoidReduction"
    empty(): () = ()
    join(a: (), b: ()): () = ()
end

(* Hack to permit any Number to work non-parametrically. *)
object SumReduction extends CommutativeMonoidReduction[\Number\]
    getter toString() = "SumReduction"
    empty(): Number = 0
    join(a: Number, b: Number): Number = a+b
end

opr SUM[\T extends Number\](): Comprehension[\T,Number,Number,Number\] =
    Comprehension[\T,Number,Number,Number\](fn x => x, SumReduction, cast[\Number\])

object ProdReduction extends CommutativeMonoidReduction[\Number\]
    getter toString() = "ProdReduction"
    empty(): Number = 1
    join(a:Number, b:Number): Number = a b
end

opr PROD[\T extends Number\](): Comprehension[\T,Number,Number,Number\] =
    Comprehension[\T,Number,Number,Number\](fn x => x, ProdReduction, cast[\Number\])

object MinReduction[\T extends StandardMin[\T\]\] extends CommutativeReduction[\T\]
    getter toString() = "MinReduction"
    simpleJoin(a, b) = a MIN b
end

opr BIG MIN[\T extends StandardMin[\T\]\](): BigReduction[\T,AnyMaybe\] =
    BigReduction[\T,AnyMaybe\](MinReduction[\T\])

object MaxReduction[\T extends StandardMax[\T\]\] extends CommutativeReduction[\T\]
    getter toString() = "MaxReduction"
    simpleJoin(a, b) = a MAX b
end

opr BIG MAX[\T extends StandardMax[\T\]\](): BigReduction[\T,AnyMaybe\] =
    BigReduction[\T,AnyMaybe\](MaxReduction[\T\])

opr BIG MINNUM(): BigReduction[\RR64,RR64\] =
    BigReduction[\RR64,RR64\](
        MapReduceReduction[\RR64\](fn (a:RR64,b:RR64):RR64 => a MINNUM b, 0.0/0.0))

opr BIG MAXNUM(): BigReduction[\RR64,RR64\] =
    BigReduction[\RR64,RR64\](
        MapReduceReduction[\RR64\](fn (a:RR64,b:RR64):RR64 => a MAXNUM b, 0.0/0.0))

(** AndReduction and OrReduction take advantage of natural zeroes for early exit. **)
object AndReduction
        extends { CommutativeMonoidReduction[\Boolean\],
                  ReductionWithZeroes[\Boolean,Boolean\] }
    getter toString() = "AndReduction"
    empty(): Boolean = true
    join(a: Boolean, b: Boolean): Boolean = a AND b
    isZero(a:Boolean): Boolean = NOT a
end

opr BIG AND[\T\](): BigReduction[\Boolean,Boolean\] =
    BigReduction[\Boolean,Boolean\](AndReduction)

object OrReduction
        extends { CommutativeMonoidReduction[\Boolean\],
                  ReductionWithZeroes[\Boolean,Boolean\] }
    getter toString() = "OrReduction"
    empty(): Boolean = false
    join(a: Boolean, b: Boolean): Boolean = a OR b
    isZero(a:Boolean): Boolean = a
end

opr BIG OR[\T\]():BigReduction[\Boolean, Boolean\] =
    BigReduction[\Boolean,Boolean\](OrReduction)

(** A reduction performing String concatenation **)
object StringReduction extends MonoidReduction[\String\]
    getter toString() = "StringReduction"
    empty(): String = ""
    join(a:String, b:String): String = a || b
end

(** A reduction performing String concatenation with a space **)
object SpaceReduction extends MonoidReduction[\String\]
    getter toString() = "StringReduction"
    empty(): String = ""
    join(a:String, b:String): String = a ||| b
end

(** A reduction performing String concatenation with newline separation **)
object NewlineReduction extends AssociativeReduction[\String\]
    getter toString() = "NewlineReduction"
    simpleJoin(a, b) = a // b
end

(** This operator performs string concatenation, first converting
    its inputs (of type Any) to String if necessary. **)
opr BIG ||(): Comprehension[\Any,String,String,String\] =
    Comprehension[\Any,String,String,String\](
        identity[\String\],StringReduction,fn (x:Any)=> "" x)

(** This operator performs string concatenation, first converting
    its inputs (of type Any) to String if necessary, and separating
    non-empty components by a space. **)
opr BIG |||(): Comprehension[\Any,String,String,String\] =
    Comprehension[\Any,String,String,String\](
        identity[\String\],SpaceReduction,fn (x:Any)=> "" x)

(** This operator performs string concatenation with newline
    separation, first converting its inputs (of type Any) to String if
    necessary. **)
opr BIG //(): Comprehension[\Any,String,String,AnyMaybe\] =
    Comprehension[\Any,String,String,AnyMaybe\](
        fn x => x,NewlineReduction,fn (x:Any)=> "" x)

(** A %MapReduceReduction% takes an associative binary function %j% on
    arguments of type %R%, and the identity of that function %z%, and
    returns the corresponding reduction. **)
object MapReduceReduction[\R\](j:(R,R)->R, z:R) extends MonoidReduction[\R\]
    getter toString()="mapReduceReduction"
    empty(): R = z
    join(a:R, b:R): R = (j)(a,b)
end

(** A %MIMapReduceReduction% takes an associative binary function %j% on
    arguments of type %R%, and the identity of that function %z%, and
    returns the corresponding reduction. **)
object MIMapReduceReduction[\R\](j:(Any,Any)->R, z:Any) extends MonoidReduction[\Any\]
    getter toString():String="MIMapReduceReduction from embiggen"
    empty(): R = z
    join(a:Any, b:Any): R = (j)(a,b)
end

(** %embiggen% takes a type %T% and an operation (either an operator
    %OP% or binary function %f% on type %T%), along with the identity
    %z% of the operation, and returns a function suitable as the right-hand
    side of the definition of the corresponding BIG operator. **)
(*
embiggen[\T,opr OP\](z:Any, g:(Reduction[\Any\],T->Any) -> Any) : T =
    g(MIOprReduction[\T,OP\](z), fn (x) => x)
*)
embiggen[\T\](j:(Any,Any)->T, z:T) : Comprehension[\T,T,Any,Any\] =
    Comprehension[\T,T,Any,Any\](fn (x) => x, MIMapReduceReduction[\T\](j,z), fn (x) => x)

(** Helpers for maps and cross products of generators.  These can be
    quite a bit more sophisticated (for example, we can hoist maps
    outwards if we think that'd be useful), but let's get this much
    working first. *)
trait SomeMappedGenerator[\F\] extends Generator[\F\] end

trait MappedGenerator[\E,F\] extends SomeMappedGenerator[\F\]
    getter g(): Generator[\E\]
    getter f(): E -> F
    getter toString(): String = g() ".map(f)"
    generate[\R\](r: Reduction[\R\], m: F->R): R =
        g().generate[\R\](r, m COMPOSE f())
    reduce(r: Reduction[\F\]): F =
        g().generate[\F\](r, f())
    map[\G\](f': F->G): SimpleMappedGenerator[\E,G\] =
        SimpleMappedGenerator[\E,G\](g(), f' COMPOSE f())
    seq(self) = SimpleMappedSeqGenerator[\E,F\](seq(g()),f())
end

object SimpleMappedGenerator[\E,F\](g0: Generator[\E\], f0: E->F)
        extends MappedGenerator[\E,F\]
    getter g() = g0
    getter f() = f0
end

object SimpleMappedIndexed[\E,F,I\](g0: Indexed[\E,I\], f0: E->F)
        extends { MappedGenerator[\E,F\], Indexed[\F,I\] }
    getter size(): ZZ32 = |g()|
    getter g() = g0
    getter f() = f0
    getter bounds(): Range[\I\] = g().bounds
    getter indexValuePairs(): Indexed[\(I,F),I\] =
        g().indexValuePairs.map[\(I,F)\](fn (i:I, e:E): F => (i,(f())(e)))
    getter indices(): Indexed[\I,I\] = g().indices
    opr |self| : ZZ32 = |g()|
    opr[i:I] : F = (f())(g()[i])
    opr[r:Range[\I\]] : Indexed[\F,I\] =
        SimpleMappedIndexed[\E,F,I\](g()[r],f())

    ivmap[\R\](ff:(I,F)->R): Indexed[\R,I\] =
        g0().ivmap[\R\](fn (i:I, e:E) => ff(i, (f())(e)))
    map[\G\](f': F->G): SimpleMappedIndexed[\E,G,I\] =
        SimpleMappedIndexed[\E,G,I\](g(), f' COMPOSE f())
end

object SimpleMappedSeqGenerator[\E,F\](g0: SequentialGenerator[\E\], f0: E->F)
        extends { MappedGenerator[\E,F\], SequentialGenerator[\F\] }
    getter g() = g0
    getter f() = f0
    getter toString() = "seq(" g0 ".map(f))"
    seq(self): SimpleMappedSeqGenerator[\E,F\] = self
end

trait FilterGenerator[\E\] extends Generator[\E\]
    getter g(): Generator[\E\]
    getter p(): E -> Condition[\()\]
    getter toString(): String = g() ".filter(p)"
    generate[\R\](r:Reduction[\R\], m: E->R): R =
        g().generate[\R\](r, fn(e:E):R => if (p())(e).holds then m(e) else r.empty() end)
    reduce(r: Reduction[\E\]): E =
        g().generate[\E\](r, fn(e:E):E => if (p())(e).holds then e else r.empty() end)
    filter(p': E -> Condition[\()\]): SimpleFilterGenerator[\E\] =
        SimpleFilterGenerator[\E\](g(), p() ANDCOND p')
    seq(self) = SimpleSeqFilterGenerator[\E\](seq(g()),p())
end

object SimpleFilterGenerator[\E\](g0:Generator[\E\], p0: E->Condition[\()\])
        extends FilterGenerator[\E\]
    getter g(): Generator[\E\]       = g0
    getter p(): E -> Condition[\()\]  = p0
end

object SimpleSeqFilterGenerator[\E\](g0: SequentialGenerator[\E\], p0: E->Condition[\()\])
        extends { FilterGenerator[\E\], SequentialGenerator[\E\] }
    getter g(): Generator[\E\]       = g0
    getter p(): E -> Condition[\()\]  = p0
    getter toString() = "seq(" g0 ".filter(p))"
    seq(self): SimpleSeqFilterGenerator[\E\] = self
end

trait NestedGenerator[\E,F\] extends Generator[\F\]
    getter g(): Generator[\E\]
    getter f(): E -> Generator[\F\]
    getter toString(): String
    generate[\R\](r: Reduction[\R\], m: F->R): R =
        g().generate[\R\](r,fn (e:E):R => (f())(e).generate[\R\](r,m))
    mapReduce[\R\](body: F->R, join:(R,R)->R, id:R): R =
        g().mapReduce[\R\](
            fn (e:E): R => (f())(e).mapReduce[\R\](body,join,id),
            join, id)
    reduce(r: Reduction[\F\]): F =
        g().generate[\F\](r,fn (e:E):F => (f())(e).reduce(r))
    reduce(j:(F,F)->F, z:F):F =
        g().mapReduce[\F\](fn (e:E): F => (f())(e).reduce(j,z), j, z)
    loop(body:F->()): () =
        g().loop(fn (e:E) => (f())(e).loop(body))
    map[\G\](h:F->G): Generator[\G\] =
        g().nest[\G\](fn (e:E): Generator[\G\] => (f())(e).map[\G\](h))
    nest[\G\](h:F->Generator[\G\]): Generator[\G\] =
        g().nest[\G\](fn (e:E):Generator[\G\] => (f())(e).nest[\G\](h))
end

object SimpleNestedGenerator[\E,F\](g0: Generator[\E\], f0: E->Generator[\F\])
        extends { NestedGenerator[\E,F\] }
    getter g() = g0
    getter f() = f0
    getter toString() = g() ".nest(f)"
    seq(self) =
        typecase ff = f() of
            (E -> SequentialGenerator[\F\]) =>
                SimpleNestedSeqGenerator[\E,F\](seq(g()),ff)
            else =>
                SimpleNestedSeqGenerator[\E,F\](seq(g()),
                     fn (e:E): SequentialGenerator[\F\] => seq(self.f()(e)))
        end
end

object SimpleNestedSeqGenerator[\E,F\]
        (g0: SequentialGenerator[\E\], f0: E->SequentialGenerator[\F\])
        extends { NestedGenerator[\E,F\], SequentialGenerator[\F\] }
    getter g() = g0
    getter f() = f0
    getter toString() = "seq(" g() ".nest(f))"
end

trait PairGenerator[\E,F\] extends Generator[\(E,F)\]
  comprises { SimplePairGenerator[\E,F\], SimplePairSeqGenerator[\E,F\] }
    getter e(): Generator[\E\]
    getter f(): Generator[\F\]
    getter toString(): String
    getter size():ZZ32 = |self|
    opr |self| : ZZ32 = |e()| |f()|
    generate[\R\](r: Reduction[\R\], m:(E,F)->R): R =
        e().generate[\R\](r, fn (a: E): R =>
                              f().generate[\R\](r, fn (b: F): R => m(a,b)))
end

object SimplePairGenerator[\E,F\](e0: Generator[\E\], f0: Generator[\F\])
        extends PairGenerator[\E,F\]
    getter e() = e0
    getter f() = f0
    getter toString()=(e0 ".cross(" f0 ")" )
    seq(self) = SimplePairSeqGenerator[\E,F\](seq(e0),seq(f0))
end

object SimplePairSeqGenerator[\E,F\](e0: Generator[\E\], f0: Generator[\F\])
        extends { PairGenerator[\E,F\], SequentialGenerator[\(E,F)\] }
    getter e() = e0
    getter f() = f0
    getter toString()=("seq(" e ".cross(" f "))" )
end

(** Helper for serializing generators naively.  This code should make
 *  obvious that naive seq is *VERY INEFFICIENT*.  It
 *  constructs a function closure whose size is proportional to
 *  %|g|%, and then executes that closure.  This trick is old hat
 *  to lambda-calculus wonks, but pretty unfamiliar to the common man.
 *
 *  Basically each element takes in the accumulated value a from the element
 *  to its left.  This is joined with the result of mapping on the
 *  current element value, and that result is returned.  Join simply
 *  reverse-composes the functions for its subtrees, so the output of the left
 *  subtree is fed to the right subtree.
 *
 *  Note that a similar trick can be used to reverse and sequentialize
 *  a generator (use forward function composition rather than reverse
 *  composition, and flipping around the join at the leaves). *)
object NaiveSeqGenerator[\E\](g: Generator[\E\])
        extends SequentialGenerator[\E\]
    getter size() = |g|
    getter toString():String = "naive seq(" g.toString ")"
    opr |self| : ZZ32 = |g|
    generate[\R\](r: Reduction[\R\], m:E->R): R = do
      rcompose(f:R->R, f':R->R): R->R = fn (x:R):R => f'(f(x))
      id(x:R):R = x
      mp(x:E):R->R = fn (a:R):R => r.join(a,m(x))
      f : R -> R = g.mapReduce[\R->R\](mp,rcompose,id)
      mt : R = r.empty()
      f mt
    end
end

(************************************************************
* Ranges
************************************************************)

(** Ranges in general represent uses of the # and : operators.
    It's mostly subtypes of Range that are interesting.

    The partial order on ranges describes containment:
      a < b iff all points in a are strictly contained in b.
 **)
trait Range[\T\]
    extends StandardPartialOrder[\Range[\T\]\]
    comprises { RangeWithLower[\T\], RangeWithUpper[\T\],
                RangeWithExtent[\T\], PartialRange[\T\] }
    excludes { Number }
    
    opr[r:Range[\T\]]: Range[\T\] = fail("Unrecognized Range " r)
    opr[_:OpenRange[\Any\]] : Range[\T\] = self
    opr[other:LowerRange[\T\]] : Range[\T\]
    opr[other:UpperRange[\T\]] : Range[\T\]
    opr[other:ExtentRange[\T\]] : Range[\T\]
    opr[other:CompleteRange[\T\]] : CompleteRange[\T\]
    opr =(self,_:Range[\T\]): Boolean = false
end

trait PartialRange[\T\] extends Range[\T\]
    comprises { OpenRange[\T\],
                LowerRange[\T\], UpperRange[\T\], ExtentRange[\T\] }
    excludes { CompleteRange[\T\] }
end

object OpenRange[\T\] extends { Range[\T\], PartialRange[\T\] }
    toString():String = "#"
    opr[_:OpenRange[\T\]] : OpenRange[\T\] = self

    opr[other:LowerRange[\T\]] : LowerRange[\T\] = other
    opr[other:UpperRange[\T\]] : UpperRange[\T\] = other
    opr[other:ExtentRange[\T\]] : ExtentRange[\T\] = other
    opr[other:CompleteRange[\T\]] : CompleteRange[\T\] = other
    opr =(self,_:OpenRange[\T\]): Boolean = true
    opr CMP(self,other:Range[\T\]): TotalComparison =
        typecase other of
            OpenRange[\T\] => EqualTo
            else => GreaterThan
        end
end

opr PARTIAL_LEXICO(a:Comparison, b:Comparison) =
    typecase a of
        Unordered => Unordered
        EqualTo => b
        else =>
            typecase b of
                Unordered => Unordered
                else => a
            end
    end

opr PARTIAL_LEXICO(a:Comparison, b:()->Comparison) =
    typecase a of
        Unordered => Unordered
        EqualTo => b()
        else =>
            typecase _ = b() of
                Unordered => Unordered
                else => a
            end
    end

(** Non-traditional partial ordering on tuples, CMP ordering elsewhere. **)
opr PCMP[\A,B\](t1:(A,B), t2:(A,B)): Comparison = do
    (a1,b1)=t1
    (a2,b2)=t2
    (a1 CMP a2) PARTIAL_LEXICO: (b1 CMP b2)
  end

opr PCMP[\A,B,C\](t1:(A,B,C), t2:(A,B,C)): Comparison = do
    (a1,b1,c1)=t1
    (a2,b2,c2)=t2
    (a1 CMP a2) PARTIAL_LEXICO: (b1 CMP b2) PARTIAL_LEXICO (c1 CMP c2)
  end

opr PCMP(a:Number, b:Number): TotalComparison = a CMP b

(** Non-traditional subtraction of bounds, elementwise on tuples **)
opr SUBTR[\A,B\](t1:(A,B), t2:(A,B)): (A,B) = do
    (a1,b1) = t1
    (a2,b2) = t2
    (a1 SUBTR a2, b1 SUBTR b2)
  end

opr SUBTR[\A,B,C\](t1:(A,B,C), t2:(A,B,C)): (A,B,C) = do
    (a1,b1,c1) = t1
    (a2,b2,c2) = t2
    (a1 SUBTR a2, b1 SUBTR b2, c1 SUBTR c2)
  end

opr SUBTR[\T extends MultiplicativeRing[\T\]\](a:T, b:T):T = a.one + a - b

(** join upper and lower bounds comparisons (total or partial) to
    yield partial containment. **)

trait RangeWithLower[\T\] extends Range[\T\]
        comprises { LowerRange[\T\], CompleteRange[\T\] }
    getter lower():T
end

rangeBoundError[\T\](r1:Range[\T\], r2: Range[\T\]): Range[\T\] =
    fail(r1 "[" r2 "] does not contain bound of indexing range")

object LowerRange[\T\](lo:T) extends { RangeWithLower[\T\], PartialRange[\T\] }
    getter lower():T = lo
    toString():String = lower() "#"
    opr[_:OpenRange[\T\]] : LowerRange[\T\] = self
    opr[other:LowerRange[\T\]] : LowerRange[\T\] =
        if other.lower >= self.lower then
            other
        else
            rangeBoundError[\T\](self,other)
        end
    opr[other:ExtentRange[\T\]] : FullRange[\T,true\] = self.lower # other.extent()
    opr[other:UpperRange[\T\]] : FullRange[\T,true\] =
        if other.upper >= self.lower then
            self.lower : other.upper
        else
            rangeBoundError[\T\](self,other)
        end
    opr[other:CompleteRange[\T\]] : CompleteRange[\T\] =
        if other.lower >= self.lower then
            other
        else
            rangeBoundError[\T\](self,other)
        end
    opr =(self, x:LowerRange[\T\]): Boolean = self.lower = x.lower
    opr CMP(self, other:Range[\T\]): Comparison =
        typecase other of
            OpenRange[\T\] => LessThan
            LowerRange[\T\] => other.lower PCMP self.lower
            CompleteRange[\T\] =>
                typecase c = other.lower PCMP self.lower of
                    EqualTo => GreaterThan
                    LessThan => Unordered
                    else => c
                end
            else => Unordered
        end
end

trait RangeWithUpper[\T\] extends Range[\T\]
        comprises { UpperRange[\T\], CompleteRange[\T\] }
    getter upper():T
end

object UpperRange[\T\](up:T) extends { RangeWithUpper[\T\], PartialRange[\T\] }
    getter upper():T = up
    toString():String = ":" upper()
    opr[_:OpenRange[\T\]] : UpperRange[\T\] = self
    opr[other:LowerRange[\T\]] : FullRange[\T,true\] = other.lower:self:upper
    opr[other:UpperRange[\T\]] : UpperRange[\T\] =
        if other.upper <= self.upper then
            other
        else
            rangeBoundError[\T\](self,other)
        end
    opr[other:ExtentRange[\T\]] : FullRange[\T,true\] =
        (self.upper SUBTR other.extent) # other.extent
    opr[other:CompleteRange[\T\]] : CompleteRange[\T\] =
        if other.upper <= self.upper then
            other
        else
            rangeBoundError[\T\](self,other)
        end
    opr =(self,x:UpperRange[\T\]): Boolean = upper() = x.upper()
    opr CMP(self, other:Range[\T\]): Comparison =
        typecase other of
            OpenRange[\T\]  => LessThan
            UpperRange[\T\] => self.upper PCMP other.upper
            CompleteRange[\T\]  =>
                typecase c = self.upper PCMP other.upper of
                    EqualTo => GreaterThan
                    LessThan => Unordered
                    else => c
                end
            else => Unordered
        end
end

trait RangeWithExtent[\T\] extends Range[\T\]
        comprises { ExtentRange[\T\], CompleteRange[\T\] }
    getter extent():T
    toString():String = "#" self.extent
end

object ExtentRange[\T\](ex:T) extends { RangeWithExtent[\T\], PartialRange[\T\] }
    getter extent():T = ex
    opr[_:OpenRange[\T\]] : ExtentRange[\T\] = self
    opr[other:LowerRange[\T\]] : FullRange[\T,true\] = other.lower # self.extent
    opr[other:UpperRange[\T\]] : FullRange[\T,true\] =
        (other.upper SUBTR self.extent) # self.extent
    opr[other:ExtentRange[\T\]] : ExtentRange[\T\] =
        if other.extent <= self.extent then
            other
        else
            rangeBoundError[\T\](self,other)
        end
    opr[other:CompleteRange[\T\]] : CompleteRange[\T\] =
        if other.extent <= self.extent then
            other
        else
            rangeBoundError[\T\](self,other)
        end
    opr =(self,x:ExtentRange[\T\]): Boolean = self.extent = x.extent
    opr CMP(self, other:Range[\T\]): Comparison =
        typecase other of
            OpenRange[\T\]   => LessThan
            ExtentRange[\T\] => self.extent PCMP other.extent
            else => Unordered
        end
end

(** The type %CompleteRange[\T\]% is really just %FullRange[\T,flag\]% for
    either setting of %flag%.  We need it for stuff like %typecase%
    and %extends% clauses since we can't bind universally quantified
    variables in those settings.  Since essentially all the methods are
    independent of the setting of flag we include them here. **)
trait CompleteRange[\T\]
        extends { RangeWithLower[\T\], RangeWithUpper[\T\],
                  RangeWithExtent[\T\], Indexed[\T,T\] }
        comprises { FullRange[\T,true\], FullRange[\T,false\] }
    getter indices(): FullRange[\T,true\] = bounds()
    getter isBounded(): Boolean
    getter isSized(): Boolean

    (** The following two methods are used in bounds checking code to
        permit ranges that are just outside the edge of the given
        range.   Because they're intended for internal use we ignore
        the setting of flag and return a CompleteRange. **)
    widenUpper(): CompleteRange[\T\]
    widenLower(): CompleteRange[\T\]

    opr[r:Range[\T\]]: CompleteRange[\T\] = fail("Unrecognized Range " r)
    opr[_:OpenRange[\T\]]: CompleteRange[\T\] = self
    (** Square-bracket indexing on a FullRange restricts that range to
         the range provided.  Restriction should behave as follows:
            - Restriction to an OpenRange is the identity.
            - An UpperRange or ExtentRange restrict the upper bound and
              extent of the range.
            - A LowerRange restricts the lower bound and extent of the range.
        Note that this makes it compatible with the square-bracket
        indexing of the Indexed trait. **)
    opr[r:LowerRange[\T\]]: CompleteRange[\T\] =
        if r.lower = self.lower then
            (* Important, also handles empty range case. *)
            self
        elif r.lower IN widenUpper() then
            r.lower:self.upper
        else
            rangeBoundError[\T\](self,r)
        end
    opr[r:UpperRange[\T\]]: CompleteRange[\T\] =
        if r.upper = self.upper then
            (* Important, also handles empty range case. *)
            self
        elif r.upper IN widenLower() then
            self.lower:r.upper
        else
            rangeBoundError[\T\](self,r)
        end
    opr[r:ExtentRange[\T\]]: CompleteRange[\T\] =
        if r.extent <= self.extent then
            self.lower#r.extent
        else
            fail(self "[" r "] is too small in some dimension")
        end
    opr[r:CompleteRange[\T\]]: CompleteRange[\T\] =
        if r.lower < self.lower then
            rangeBoundError[\T\](self,r)
        elif r.upper > self.upper then
            rangeBoundError[\T\](self,r)
        else
            r
        end

    toString():String =
        if self.isBounded then
            self.lower ":" self.upper
        else
            self.lower "#" self.extent
        end

    opr =(self, other:CompleteRange[\T\]): Boolean =
    (*    isBounded()=other.isBounded() AND 
         APB 2008.08.20: this seems bogus: isBounded tells us only whether we were created 
         using an extent or an upper bound, and should not effect the comparison *)
         lower()=other.lower() AND extent() = other.extent()    
         
    opr CMP(self, other:Range[\T\]): Comparison =
        typecase other of
            CompleteRange[\T\] =>
                a = self.lower PCMP other.lower
                b = self.upper PCMP other.upper
                typecase _ = (a,b) of
                    (Unordered,Comparison) => Unordered
                    (Comparison,Unordered) => Unordered
                    (LessThan,LessThan) => Unordered
                    (GreaterThan,GreaterThan) => Unordered
    (*  This case seems unncessary; the (EqualTo,Comparison) case will cover us    --- APB  2008.08.20
                    (EqualTo,EqualTo) => EqualTo
                        if self.isBounded = other.isBounded then
                            EqualTo
                        else
                            Unordered
                        end
   *)
                    (EqualTo,Comparison) => b
                    (GreaterThan,Comparison) => LessThan
                    else (* LessThan,Comparison *) => GreaterThan
                end
            else => INVERSE (other CMP self)
        end
end

(** A %FullRange% is either bounded (generated from a lower and upper
    bound, as in %l:u%) or sized (generated from a lower bound and a
    size, as in %l#s%).  These two kinds of range behave identically
    in most settings; the only time they behave differently is when we
    construct a strided range: both %l:u:i% and %l#s:i% generate
    elements %l, l+i, l+2i, ...%, but the former is bounded and
    generates no element larger than %i% and the latter is sized and
    generates exactly %s% elements.

    Data structure bounds are always bounded. **)
trait FullRange[\T, bool bounded\] extends CompleteRange[\T\]
        comprises { ScalarRange[\T,bounded\], TupleRange[\T,bounded\] }
    getter isBounded(): Boolean = bounded
    getter isSized(): Boolean = NOT bounded
end

rangeIndexFailure[\T\](r:CompleteRange[\T\],i:T): () =
    fail(r "[" i "] out of range.")

(** ParRanges and SeqRanges, the simplest and most basic of all
    generators.  They generate |self| numbers, incrementing by 1,
    starting with lower() (that's the canonical ordering). *)
trait ScalarRange[\N extends Integral[\N\], bool bounded\]
        extends { FullRange[\N, bounded\] }
        comprises { ParRange[\N, bounded\], SeqRange[\N, bounded\] }
    getter bounds():FullRange[\N,true\] = ParRange[\N,true\](0,extent())
    getter upper():N = lower() + extent() - 1
    getter size():ZZ32 = |self|
    widenUpper(): ScalarRange[\N,false\] = lower()#(extent()+1)
    widenLower(): ScalarRange[\N,false\] = (lower()-1)#(extent()+1)
    opr |self| : ZZ32 = narrow(0 MAX extent())
    opr[i:N]: N =
        if 0 <= i < extent() then
            i + lower()
        else
            rangeIndexFailure[\N\](self,i)
        end
    opr IN(a:N, self): Boolean = do
        t = a - lower()
        0 <= t < extent()
      end
end

value object ParRange[\N extends Integral[\N\], bool bounded\](lo:N, ex:N)
        extends ScalarRange[\N,bounded\]
    getter extent() = ex
    getter lower() = lo
    getter indexValuePairs() =
        ParRange[\N,bounded\](0,ex).map[\(N,N)\](fn (n:N):(N,N) => (n,n+lo))

    seq(self) = SeqRange[\N,bounded\](lo,ex)
    generate[\R\](r: Reduction[\R\], m: N->R): R =
        if ex < 1 then
            r.empty()
        else
            gen(l,s) = if s>=2 then
                           (s1,s2) = partition(s)
                           r.join(gen(l,s1),gen(l+s1,s2))
                       else
                           m(l)
                       end
            gen(lo,ex)
        end
    mapReduce[\R\](m: N->R, j:(R,R)->R, z:R): R =
        if ex < 1 then
            z
        else
            gen(l,s) = if s>=2 then
                           (s1,s2) = partition(s)
                           j(gen(l,s1),gen(l+s1,s2))
                       else
                           m(l)
                       end
            gen(lo,ex)
        end
    reduce(j:(N,N)->N, z:N):N =
        if ex < 1 then
            z
        else
            gen(l,s) = if s>=2 then
                           (s1,s2) = partition(s)
                           j(gen(l,s1),gen(l+s1,s2))
                       else
                           l
                       end
            gen(lo,ex)
        end
    reduce(r: Reduction[\N\]):N =
        if ex < 1 then
            r.empty()
        else
            gen(l,s) = if s>=2 then
                           (s1,s2) = partition(s)
                           r.join(gen(l,s1),gen(l+s1,s2))
                       else
                           l
                       end
            gen(lo,ex)
        end
    loop(f:N->()): () =
        if ex < 1 then
            ()
        else
            gen(l,s) = if s>=2 then
                           (s1,s2) = partition(s)
                           (gen(l,s1),gen(l+s1,s2))
                           ()
                       else
                           f(l)
                       end
            gen(lo,ex)
        end
end
value object SeqRange[\N extends Integral[\N\], bool bounded\](lo:N, ex:N)
        extends { ScalarRange[\N,bounded\], SequentialGenerator[\N\] }
    getter extent() = ex
    getter lower() = lo
    getter toString() =
        ("seq(" lo (if bounded then ":" upper() else "#" ex ")"))
    getter indexValuePairs() =
        SeqRange[\N,bounded\](0,ex).map[\(N,N)\](fn (n:N):(N,N) => (n,n+lo))
    generate[\R\](r: Reduction[\R\], m: N->R): R = do
        h = lo + ex
        l : N := lo
        acc : R = r.empty()
        (* DON'T DEFINE THIS USING RECURSION!!! *)
        while l < h do
            b = m(l)
            acc := r.join(acc, b)
            l += 1
        end
        acc
      end
    loop(f: N->()): () = do
        h = lo + ex
        l : N := lo
        (* DON'T DEFINE THIS USING RECURSION!!! *)
        while l < h do
            f(l)
            l += 1
        end
      end
end

trait TupleRange[\T, bool bounded\] extends FullRange[\T,bounded\]
end

object Tuple2Range[\N extends Integral[\N\], M extends Integral[\M\], bool bounded\]
                  (l1:N,l2:M,x1:N,x2:M)
        extends { TupleRange[\(N,M),bounded\], DelegatedIndexed[\(N,M),(N,M)\] }
    getter bounds():Tuple2Range[\N,M,true\] = Tuple2Range[\N,M,true\](0,0,x1,x2)
    getter extent():(N,M) = (x1,x2)
    getter lower():(N,M) = (l1,l2)
    getter upper():(N,M) = (l1+x1-1, l2+x2-1)
    getter size():(N,M) = |self|
    getter generator():Generator[\(N,M)\] =
        ParRange[\N,bounded\](l1,x1).cross[\M\](ParRange[\M,bounded\](l2,x2))
    getter toString(): String =
        ("(" l1 "," l2
         (if bounded then
             (u1,u2) = upper()
             "):(" u1 "," u2 ")"
         else
             ")#(" x1 "," x2 ")"
         end))
    getter indexValuePairs() =
        bounds().map[\((N,M),(N,M))\](
            fn (n:N,m:M):((N,M),(N,M)) => ((n,m),(n+l1,m+l2)))
    widenUpper(): Tuple2Range[\N,M,true\] = Tuple2Range[\N,M,true\](l1,l2,x1+1,x2+1)
    widenLower(): Tuple2Range[\N,M,true\] = Tuple2Range[\N,M,true\](l1-1,l2-1,x1+1,x2+1)
    opr |self| : ZZ32 = narrow((0 MAX x1) (0 MAX x2))
    opr[i:(N,M)]: (N,M) = do
        (i1,i2) = i
        if 0 <= i1 < x1 AND 0 <= i2 < x2 then
            (i1+l1,i2+l2)
        else
            rangeIndexFailure[\(N,M)\](self,i)
        end
      end
    opr IN(t:(N,M), self): Boolean = do
        (v1,v2) = t
        0 <= v1-l1 < x1 AND 0 <= v2-l2 < x2
      end
    opr =(self,r:Tuple2Range[\N,M,bounded\]): Boolean =
        l1=r.l1 AND l2=r.l2 AND x1=r.x1 AND x2=r.x2
end

object Tuple3Range[\N1 extends Integral[\N1\], N2 extends Integral[\N2\], N3 extends Integral[\N3\],
                   bool bounded\]
                  (l1:N1,l2:N2,l3:N3,x1:N1,x2:N2,x3:N3)
        extends { TupleRange[\(N1,N2,N3),bounded\],
                  DelegatedIndexed[\(N1,N2,N3),(N1,N2,N3)\] }
    getter bounds():Tuple3Range[\N1,N2,N3,true\] =
        Tuple3Range[\N1,N2,N3,true\]((0,0,0),(x1,x2,x3))
    getter extent():(N1,N2,N3) = (x1,x2,x3)
    getter lower():(N1,N2,N3) = (l1,l2,l3)
    getter upper():(N1,N2,N3) = (l1+x1-1,l2+x2-1,l3+x3-1)
    getter size():(N1,N2,N3) = |self|
    getter generator():Generator[\(N1,N2,N3)\] =
        ParRange[\N1,bounded\](l1,x1).cross[\N2\](ParRange[\N2,bounded\](l2,x2)).cross[\N3\](
            ParRange[\N3,bounded\](l3,x3)).map[\(N1,N2,N3)\](
                fn (t:(N1,N2),n3:N3):(N1,N2,N3) => do (n1,n2)=t; (n1,n2,n3) end)
    getter toString(): String =
        ("(" l1 "," l2 "," l3
         (if bounded then
             (u1,u2,u3) = upper()
             "):(" u1 "," u2 "," u3 ")"
         else
             ")#(" x1 "," x2 "," x3 ")"
         end))
    getter indexValuePairs() =
        (0#x1).cross(0#x2).cross(0#x3).map[\((N1,N2,N3),(N1,N2,N3))\](
            fn (t:(N1,N2),n3:N3):((N1,N2,N3),(N1,N2,N3)) =>
                do (n1,n2)=t; ((n1,n2,n3),(n1+l1,n2+l2,n3+l3)) end)
    widenUpper(): Tuple3Range[\N1,N2,N3,true\] =
        Tuple3Range[\N1,N2,N3,true\](l1,l2,l3,x1+1,x2+1,x3+1)
    widenLower(): Tuple3Range[\N1,N2,N3,true\] =
        Tuple3Range[\N1,N2,N3,true\](l1-1,l2-1,l3-1,x1+1,x2+1,x3+1)
    opr |self| : (N1,N2,N3) = narrow((0 MAX x1) (0 MAX x2) (0 MAX x3))
    opr[i:(N1,N2,N3)]: (N1,N2,N3) = do
        (i1,i2,i3) = i
        if 0 <= i1 < x1 AND 0 <= i2 < x2 AND 0 <= i3 < x3 then
            (i1+l1,i2+l2,i3+l3)
        else
            rangeIndexFailure[\(N1,N2,N3)\](self,i)
        end
      end
    opr IN(t:(N1,N2,N3), self): Boolean = do
        (v1,v2,v3) = t
        0 <= v1-l1 < x1 AND 0 <= v2-l2 < x2 AND 0 <= v3-l3 < x3
      end
    opr =(self,r:Tuple3Range[\N1,N2,N3,bounded\]): Boolean =
        l1=r.l1 AND l2=r.l2 AND x1=r.x1 AND x2=r.x2 AND l3=r.l3 AND x3=r.x3
end

(** Helpers for # to get the type instantiation "right".  We pass in bogus ZZ32's to
    ensure that the result type is at least ZZ32. **)
sizedRange[\I extends AnyIntegral\](_:I,lo:I,ex:I): ParRange[\I,false\] = ParRange[\I,false\](lo,ex)
sized2Range[\I extends AnyIntegral, J extends AnyIntegral\]
           (_:I,_:J,l1:I,l2:J,ex1:I,ex2:J): Tuple2Range[\I,J,false\] =
    Tuple2Range[\I,J,false\](l1,l2,ex1,ex2)
sized3Range[\I extends AnyIntegral, J extends AnyIntegral, K extends AnyIntegral\]
     (_:I,_:J,_:K,l1:I,l2:J,l3:K,ex1:I,ex2:J,ex3:K): Tuple3Range[\I,J,K,false\] =
    Tuple3Range[\I,J,K,false\](l1,l2,l3,ex1,ex2,ex3)

(** The # and : operators serve as factories for parallel ranges. **)
opr #[\I extends AnyIntegral\](lo:I, ex:I): ParRange[\I,false\] = sizedRange(0 asif ZZ32,lo,ex)
(*
opr #(lo:IntLiteral, ex:IntLiteral): ParRange[\ZZ32\] = ParRange[\ZZ32,false\](lo,ex)
*)
opr #[\I extends AnyIntegral, J extends AnyIntegral\]
     (lo:(I,J), ex:(I,J)): Tuple2Range[\I,J,false\] =
    do (l1,l2)=lo; (x1,x2)=ex; sized2Range(0 asif ZZ32,0 asif ZZ32,l1,l2,x1,x2) end
opr #[\I extends AnyIntegral, J extends AnyIntegral, K extends AnyIntegral\]
     (lo:(I,J,K), ex:(I,J,K)): Tuple3Range[\I,J,K,false\] =
    do (l1,l2,l3)=lo; (x1,x2,x3)=ex
       sized3Range(0 asif ZZ32,0 asif ZZ32,0 asif ZZ32,l1,l2,l3,x1,x2,x3)
    end

(** Helpers for : to get the type instantiation "right".  We pass in bogus ZZ32's to
    ensure that the result type is at least ZZ32.  **)
boundedRange[\I extends AnyIntegral\](_:I,lo:I,ex:I): ParRange[\I,true\] = ParRange[\I,true\](lo,ex)
bounded2Range[\I extends AnyIntegral, J extends AnyIntegral\]
           (_:I,_:J,l1:I,l2:J,ex1:I,ex2:J): Tuple2Range[\I,J,true\] =
    Tuple2Range[\I,J,true\](l1,l2,ex1,ex2)
bounded3Range[\I extends AnyIntegral, J extends AnyIntegral, K extends AnyIntegral\]
     (_:I,_:J,_:K,l1:I,l2:J,l3:K,ex1:I,ex2:J,ex3:K): Tuple3Range[\I,J,K,true\] =
    Tuple3Range[\I,J,K,true\](l1,l2,l3,ex1,ex2,ex3)

opr :[\I extends AnyIntegral\](lo:I, hi:I): ParRange[\I,true\] =
        boundedRange(0 asif ZZ32,lo,hi-lo+1)
(*
opr :(lo:IntLiteral, hi:IntLiteral): ParRange[\ZZ32,true\] =
        ParRange[\ZZ32,true\](lo,hi-lo+1)
*)
opr :[\I extends AnyIntegral, J extends AnyIntegral\]
     (lo:(I,J), hi:(I,J)): Tuple2Range[\I,J,true\] =
    do (l1,l2)=lo; (h1,h2)=hi; bounded2Range(0 asif ZZ32,0 asif ZZ32,l1,l2,h1-l1+1,h2-l2+1) end
opr :[\I extends AnyIntegral, J extends AnyIntegral, K extends AnyIntegral\]
     (lo:(I,J,K), hi:(I,J,K)): Tuple3Range[\I,J,K,true\] = do
    (l1,l2,l3)=lo; (h1,h2,h3)=hi
    bounded3Range(0 asif ZZ32,0 asif ZZ32,0 asif ZZ32,l1,l2,l3,h1-l1+1,h2-l2+1,h3-l3+1)
  end

lowerRange[\I extends AnyIntegral\](_:I,x:I):LowerRange[\I\] = LowerRange[\I\](x)
lower2Range[\I extends AnyIntegral, J extends AnyIntegral\](_:I,_:J,x:I,y:J):
    LowerRange[\(I,J)\] = LowerRange[\(I,J)\](x,y)
lower3Range[\I extends AnyIntegral, J extends AnyIntegral, K extends AnyIntegral\]
           (_:I,_:J,_:K,x:I,y:J,z:K): LowerRange[\(I,J,K)\] =
    LowerRange[\(I,J,K)\](x,y,z)

extentRange[\I extends AnyIntegral\](_:I,x:I):ExtentRange[\I\] = ExtentRange[\I\](x)
extent2Range[\I extends AnyIntegral, J extends AnyIntegral\](_:I,_:J,x:I,y:J):
    ExtentRange[\(I,J)\] = ExtentRange[\(I,J)\](x,y)
extent3Range[\I extends AnyIntegral, J extends AnyIntegral, K extends AnyIntegral\]
            (_:I,_:J,_:K,x:I,y:J,z:K): ExtentRange[\(I,J,K)\] =
    ExtentRange[\(I,J,K)\](x,y,z)

upperRange[\I extends AnyIntegral\](_:I,x:I):UpperRange[\I\] = UpperRange[\I\](x)
upper2Range[\I extends AnyIntegral, J extends AnyIntegral\](_:I,_:J,x:I,y:J):
    UpperRange[\(I,J)\] = UpperRange[\(I,J)\](x,y)
upper3Range[\I extends AnyIntegral, J extends AnyIntegral, K extends AnyIntegral\]
           (_:I,_:J,_:K,x:I,y:J,z:K): UpperRange[\(I,J,K)\] =
    UpperRange[\(I,J,K)\](x,y,z)

(** Factories for incomplete ranges **)
opr (x:I)#[\I extends AnyIntegral\] : LowerRange[\I\] = lowerRange(0 asif ZZ32, x)
opr (x:I,y:J)#[\I extends AnyIntegral, J extends AnyIntegral\] : LowerRange[\(I,J)\] =
    lower2Range(0 asif ZZ32,0 asif ZZ32, x,y)
opr (x:I,y:J,z:K)#[\I extends AnyIntegral, J extends AnyIntegral, K extends AnyIntegral\] :
         LowerRange[\(I,J,K)\] = lower3Range(0 asif ZZ32,0 asif ZZ32,0 asif ZZ32, x,y,z)
opr (x:I):[\I\] : LowerRange[\I\] = x#
opr #[\I extends AnyIntegral\](x:I) : ExtentRange[\I\] = extentRange(0 asif ZZ32, x)
opr #[\I extends AnyIntegral, J extends AnyIntegral\](xy:(I,J)) : ExtentRange[\(I,J)\] =
    do (x,y) = xy; extent2Range(0 asif ZZ32,0 asif ZZ32, x,y) end
opr #[\I extends AnyIntegral, J extends AnyIntegral, K extends AnyIntegral\](xyz:(I,J,K)) :
         ExtentRange[\(I,J,K)\] =
    do (x,y,z) = xyz; extent3Range(0 asif ZZ32,0 asif ZZ32,0 asif ZZ32, x,y,z) end
opr :[\I extends AnyIntegral\](x:I) : UpperRange[\I\] = upperRange(0 asif ZZ32, x)
opr :[\I extends AnyIntegral, J extends AnyIntegral\](xy:(I,J)) : UpperRange[\(I,J)\] =
    do (x,y) = xy; upper2Range(0 asif ZZ32,0 asif ZZ32, x,y) end
opr :[\I extends AnyIntegral, J extends AnyIntegral, K extends AnyIntegral\](xyz:(I,J,K)) :
         UpperRange[\(I,J,K)\] =
    do (x,y,z) = xyz; upper3Range(0 asif ZZ32,0 asif ZZ32,0 asif ZZ32, x,y,z) end
(*
opr (x:I)#[\I\] : LowerRange[\I\] = LowerRange[\I\](x)
opr (x:I):[\I\] : LowerRange[\I\] = LowerRange[\I\](x)
opr #[\I\](x:I) : ExtentRange[\I\] = ExtentRange[\I\](x)
opr :[\I\](x:I) : UpperRange[\I\] = UpperRange[\I\](x)
*)
opr #() : OpenRange[\Any\] = OpenRange[\Any\]
opr :() : OpenRange[\Any\] = OpenRange[\Any\]

trait String extends { StandardTotalOrder[\String\], ZeroIndexed[\Char\], DelegatedIndexed[\Char,ZZ32\] }
             excludes {Number, Char, AnyMultiplicativeRing, AnyList}
    getter size() : ZZ32
    getter toString() : String = self
    getter indices() : FullRange[\ZZ32,true\] = self.bounds
    (* This is the "dumb" generator.  It always works, but can be "Big O" inefficient, e.g., for
       CatString trees.  Sub-traits should override it with equivalent but more efficeint code. *)
    getter generator() : Generator[\Char\] =
        indices().map[\Char\](fn (i:ZZ32):Char => self[i])       
    getter depth() = 0
    
    showStructure() = showStructure(0)
     
    opr |self| : ZZ32 = self.size
       
    opr CASE_INSENSITIVE_CMP(self, other:String): TotalComparison =
      ((BIG LEXICO[i <- 0#(|self| MIN |other|) ] self.get(i) CASE_INSENSITIVE_CMP other.get(i)) LEXICO (|self| CMP |other|))
    
    opr CMP(self, other:String): TotalComparison =
      ((BIG LEXICO[i <- 0#(|self| MIN |other|) ] self.get(i) CMP other.get(i)) LEXICO (|self| CMP |other|))
      
    opr [i:ZZ32]: Char =
        (** As a convenience, we permit LowerRange indexing to go 1 past the bounds
        of the String, returning the empty String, in order to permit some convenient
        String-trimming idioms. **)
        if 0 <= i < |self|
        then self.get(i)
        else 
            (* println "Index " i " not in bounds " self.bounds *)
            throw IndexOutOfBounds(0,(|self|)-1,i)
        end

    opr[r0:Range[\ZZ32\]] : String =  do
        r1 = self.bounds[r0]     (* This will complete r0 if it is incomplete, and throw a bounds exception when necessary *)
        self.uncheckedSubstring(r1)
      end 

    (** The operator %||% with at least one String argument converts to String and
        appends **)
    opr ||(self, b:String):String
    opr ||(self, b:Number):String = self || b.toString
    opr ||(self, c:Char):String
    opr ||(self, b:()):String = self || "()"
    opr ||(self, b:(Any,Any)):String = do (i,j) = b; self || "(" || i || "," j || ")" end
    opr ||(self, b:(Any,Any,Any)):String =
        do (i,j,k) = b; self || "(" || i || "," j || "," k || ")" end
    opr ||(self, b:(Any,Any,Any,Any)):String =
        do (i,j,k,l) = b; self || "(" || i || "," j || "," k || "," l || ")" end
    opr ||(self, b:(Any,Any,Any,Any,Any)):String =
        do (i,j,k,l,m) = b; self || "(" || i || "," j || "," k || "," l || "," m || ")" end
    opr ||(self, b:(Any,Any,Any,Any,Any,Any)):String =
        do (i,j,k,l,m,n) = b; self || "(" || i || "," j || "," k || "," l || "," m || "," n || ")" end
    opr ||(a:Any, self):String = a.toString || self
    opr ||(self, b:Any):String = self || b.toString

    (** The operator %|||% with at least one String argument converts to string,
        then appends with a whitespace separator unless one of the two arguments is
        empty.  If there is an empty argument, the other argument is returned. **)
    opr |||(self, b:String): String =
        if |self| = 0 then b
        elif |b| = 0 then self
        else self || " " || b
        end
    opr |||(self, b:Any): String = self ||| (""||b)
    opr |||(a:Any, self): String = (""||a) ||| self

    (** Right now for backward compatibility juxtaposition works like %||% **)
    opr juxtaposition(a:Any, self):String = (""||a) || self
    opr juxtaposition(self, b:String):String = self || b
    opr juxtaposition(self, b:Any):String = self || b

    (** opr %//% appends with a single newline separator.  Note that as
        with %||% and %|||% at least one argument must be a String.
        There is a postfix version in FortressLibrary that performs
        string conversion, but the prefix version must occur before a
        String. **)
    opr //(self) : String = "\n"||self
    opr //(self, a:String): String = self||"\n"||a
    opr //(self, a:Any): String    = self||"\n"||a
    opr //(a:Any, self): String    = a||"\n"||self

    (** opr %///% appends with a double newline separator.  Note that, as
        with %||% and %|||%, at least one argument must be a String.
        There is a postfix version in FortressLibrary that performs
        string conversion, but the prefix version must occur before a
        String. **)
    opr ///(self) : String = "\n\n"||self
    opr ///(self, a:String): String = self||"\n\n"||a
    opr ///(self, a:Any): String    = self||"\n\n"||a
    opr ///(a:Any, self): String    = a||"\n\n"||self

  end

(***********************************************************
* Numeric primitives
************************************************************)

nanoTime():ZZ64 = builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Long$NanoTime")
printTaskTrace():() = builtinPrimitive("com.sun.fortress.interpreter.glue.prim.StringPrim$PrintTaskTrace")

__globalTimeInformation: ZZ64 := 0
recordTime(dummy: Any): () = do __globalTimeInformation := nanoTime() end
printTime(dummy: Any): () = do
    r = nanoTime()
    e = r - __globalTimeInformation
    __globalTimeInformation := r
    secs: ZZ64 = (e+500000) DIV 1000000
    println("Operation took " secs "ms")
  end

fromRawBits(a:ZZ64):RR64 = builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$FromRawBits")

random(a:Number):RR64 = builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Float$Random")

char(a:ZZ32):Char = builtinPrimitive("com.sun.fortress.interpreter.glue.prim.Char$Chr")


print(a:Any):() = print(a.toString)
println(a:Any):() = println(a.toString)
(* 0-argument versions handle passing of () to single-argument versions. *)
print():() = print("")
println():() = println("")
print(a: JavaString): () = 
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.StringPrim$Print")
println(a: JavaString): () = 
        builtinPrimitive("com.sun.fortress.interpreter.glue.prim.StringPrim$Println")

forDigit(x:ZZ32, radix:ZZ32): Maybe[\Char\] =
    if (radix<2 OR radix>16 OR x<0 OR x>=radix) then
        Nothing[\Char\]
    elif (radix <= 12) then
        forDigit(x, "0123456789xe")
    else
        forDigit(x, "0123456789abcdef")
    end

forDigit(x:ZZ32, radixString:String): Maybe[\Char\] =
    if( 0 <= x < |radixString| ) then
        Just[\Char\](radixString[x])
    else
        Nothing[\Char\]
    end


(** Convert string to arbitrary integral type.  String may contain
    leading - or +.  Valid radixes are 2,8,10,12, and 16.  No spaces
    may occur in the string, and it must have at least one digit.

    Right now we're oblivious to overflow!  This ought to be fixed.

    This ought to support arbitrary integer types, but right now
    there's no clean way to convert all the arithmetic to use a
    provided type without having an instance of that type in hand.
    Unsigned types in particular are a bit sticky.
 **)
strToInt(s: String, radix: ZZ32): ZZ32 = do
    (i0, neg) =
        if s[0] = '-' then (1, true)
        elif s[0] = '+' then (1, false)
        else (0,false) end
    n = |s| - 1
    r = SUM[ i <- i0:n] if d <- digit(s[i], radix) then
                          d radix^(n-i)
                      else
                          fail("strToInt(" s.toString "," radix "): non-digit " s[i].toString)
                      end
    if neg then -r else r end
  end

(** Radix-10 digit conversion.  All caveats of the flexible-radix conversion above apply. *)
strToInt(s: String): ZZ32 = strToInt(s,10)

match(regex:String,some:String):Boolean = builtinPrimitive("com.sun.fortress.interpreter.glue.prim.StringPrim$Match")

(** postfix %//% appends a single newline separator.  The postfix version
    performs implicit conversion to String, which is why it's a
    standalone top-level function.  Infix and prefix versions are defined by String. **)
opr (x:Any)// : String = x||"\n"

(** postfix %///% appends a double newline separator.  The postfix version
    performs implicit conversion to String, which is why it's a
    standalone top-level function.  Infix and prefix versions are defined by String. **)
opr (x:Any)/// : String = x||"\n\n"

(* A way to get environment information from inside of fortress *)
getEnvironment(name:String, defaultValue:String):String = builtinPrimitive("com.sun.fortress.interpreter.glue.prim.StringPrim$GetEnvironment")

(* These are useful temporary hacks for debugging multi-threaded programs *)
printThreadInfo(a:String):() = builtinPrimitive("com.sun.fortress.interpreter.glue.prim.StringPrim$PrintThreadInfo")
printThreadInfo(a:Number):() = builtinPrimitive("com.sun.fortress.interpreter.glue.prim.StringPrim$PrintThreadInfo")
throwError(a:String):() = builtinPrimitive("com.sun.fortress.interpreter.glue.prim.StringPrim$ThrowError")

opr SEQV(a:Any, b:Any):Boolean =
    builtinPrimitive("com.sun.fortress.interpreter.glue.prim.AnyPrim$SEquiv")

opr  XOR(a:Boolean, b:Boolean):Boolean = if a then NOT b else b end
opr  OR(a:Boolean, b:Boolean):Boolean = if a then true else b end
opr AND(a:Boolean, b:Boolean):Boolean = if a then b else false end
opr  OR(a:Boolean, b:()->Boolean):Boolean = if a then true else b() end
opr AND(a:Boolean, b:()->Boolean):Boolean = if a then b() else false end
opr NOT(a:Boolean):Boolean = if a then false else true end
opr ->(a: Boolean, b:Boolean):Boolean = if a then b else true end
opr ->(a: Boolean, b:()->Boolean):Boolean = if a then b() else true end
opr <->(a: Boolean, b:Boolean):Boolean = a=b

true : Boolean = () SEQV ()
false : Boolean = true SEQV ()

opr +[\T extends Number\](x:T):T = x

(* Shouldn't these operators have to extend something? A,B,C? *)

opr =[\A,B\](t1:(A,B), t2:(A,B)): Boolean = do
    (a1,b1) = t1
    (a2,b2) = t2
    (a1=a2) AND: (b1=b2)
  end

opr <[\A,B\](t1:(A,B), t2:(A,B)): Boolean = do
    (a1,b1) = t1
    (a2,b2) = t2
    typecase _ = a1 CMP a2 of
        LessThan => true
        EqualTo => b1 < b2
        GreaterThan => false
    end
  end

opr <=[\A,B\](t1:(A,B), t2:(A,B)): Boolean = do
    (a1,b1) = t1
    (a2,b2) = t2
    typecase _ = a1 CMP a2 of
        LessThan => true
        EqualTo => b1 <= b2
        GreaterThan => false
    end
  end

opr >[\A,B\](t1:(A,B), t2:(A,B)): Boolean = do
    (a1,b1) = t1
    (a2,b2) = t2
    typecase _ = a1 CMP a2 of
        LessThan => false
        EqualTo => b1 > b2
        GreaterThan => true
    end
  end

opr >=[\A,B\](t1:(A,B), t2:(A,B)): Boolean = do
    (a1,b1) = t1
    (a2,b2) = t2
    typecase _ = a1 CMP a2 of
        LessThan => false
        EqualTo => b1 >= b2
        GreaterThan => true
    end
  end

opr CMP[\A,B\](t1:(A,B), t2:(A,B)): Comparison = do
    (a1,b1) = t1
    (a2,b2) = t2
    (a1 CMP a2) LEXICO: (b1 CMP b2)
  end

opr =[\A,B,C\](t1:(A,B,C), t2:(A,B,C)): Boolean = do
    (a1,b1,c1) = t1
    (a2,b2,c2) = t2
    (a1=a2) AND: (b1=b2) AND: (c1=c2)
  end

opr <[\A,B,C\](t1:(A,B,C), t2:(A,B,C)): Boolean = do
    (a1,b1,c1) = t1
    (a2,b2,c2) = t2
    typecase _ = (a1 CMP a2) LEXICO: (b1 CMP b2) of
        LessThan => true
        EqualTo => c1 < c2
        GreaterThan => false
    end
  end

opr <=[\A,B,C\](t1:(A,B,C), t2:(A,B,C)): Boolean = do
    (a1,b1,c1) = t1
    (a2,b2,c2) = t2
    typecase _ = (a1 CMP a2) LEXICO: (b1 CMP b2) of
        LessThan => true
        EqualTo => c1 <= c2
        GreaterThan => false
    end
  end

opr >[\A,B,C\](t1:(A,B,C), t2:(A,B,C)): Boolean = do
    (a1,b1,c1) = t1
    (a2,b2,c2) = t2
    typecase _ = (a1 CMP a2) LEXICO: (b1 CMP b2) of
        LessThan => false
        EqualTo => c1 > c2
        GreaterThan => true
    end
  end

opr >=[\A,B,C\](t1:(A,B,C), t2:(A,B,C)): Boolean = do
    (a1,b1,c1) = t1
    (a2,b2,c2) = t2
    typecase _ = (a1 CMP a2) LEXICO: (b1 CMP b2) of
        LessThan => false
        EqualTo => c1 >= c2
        GreaterThan => true
    end
  end

opr CMP[\A,B,C\](t1:(A,B,C), t2:(A,B,C)): Comparison = do
    (a1,b1,c1) = t1
    (a2,b2,c2) = t2
    (a1 CMP a2) LEXICO: (b1 CMP b2) LEXICO: (c1 CMP c2)
  end

object Box[\T\](var val : T)
end

end