asymptote 1.86

git-svn-id: svn://tug.org/texlive/trunk@15218 c570f23f-e606-0410-a88d-b1316a301751

1 files changed, 2351 insertions, 0 deletions

diff --git a/Build/source/utils/asymptote/runarray.cc b/Build/source/utils/asymptote/runarray.cc
new file mode 100644
index 00000000000..c67e69db6f5
--- /dev/null
+++ b/Build/source/utils/asymptote/runarray.cc
@@ -0,0 +1,2351 @@
+/***** Autogenerated from runarray.in; changes will be overwritten *****/
+
+#line 1 "runtimebase.in"
+/*****
+ * runtimebase.in
+ * Andy Hammerlindl  2009/07/28
+ *
+ * Common declarations needed for all code-generating .in files.
+ *
+ *****/
+
+
+#line 1 "runarray.in"
+/*****
+ * runarray.in
+ *
+ * Runtime functions for array operations.
+ *
+ *****/
+
+#line 1 "runtimebase.in"
+#include "stack.h"
+#include "types.h"
+#include "builtin.h"
+#include "entry.h"
+#include "errormsg.h"
+#include "array.h"
+#include "triple.h"
+#include "callable.h"
+
+using vm::stack;
+using vm::error;
+using vm::array;
+using vm::callable;
+using types::formal;
+using types::function;
+using camp::triple;
+
+#define PRIMITIVE(name,Name,asyName) using types::prim##Name;
+#include <primitives.h>
+#undef PRIMITIVE
+
+typedef double real;
+
+void unused(void *);
+
+namespace run {
+array *copyArray(array *a);
+array *copyArray2(array *a);
+array *copyArray3(array *a);
+
+double *copyArrayC(const array *a, size_t dim=0);
+double *copyArray2C(const array *a, bool square=true, size_t dim2=0);
+
+triple *copyTripleArrayC(const array *a, size_t dim=0);
+triple *copyTripleArray2C(const array *a, bool square=true, size_t dim2=0);
+double *copyTripleArray2Components(array *a, bool square=true, size_t dim2=0);
+}
+
+function *realRealFunction();
+
+// Return the component of vector v perpendicular to a unit vector u.
+inline triple perp(triple v, triple u)
+{
+  return v-dot(v,u)*u;
+}
+
+#define CURRENTPEN processData().currentpen
+
+#line 20 "runarray.in"
+#include "array.h"
+#include "arrayop.h"
+#include "triple.h"
+#include "path3.h"
+#include "Delaunay.h"
+
+#ifdef HAVE_LIBFFTW3
+#include "fftw++.h"
+#endif
+
+using namespace camp;
+using namespace vm;
+
+typedef array boolarray;
+typedef array Intarray;
+typedef array Intarray2;
+typedef array realarray;
+typedef array realarray2;
+typedef array pairarray;
+typedef array triplearray2;
+
+using types::booleanArray;
+using types::IntArray;
+using types::IntArray2;
+using types::realArray;
+using types::realArray2;
+using types::pairArray;
+using types::tripleArray2;
+
+typedef callable callableReal;
+
+void outOfBounds(const char *op, size_t len, Int n)
+{
+  ostringstream buf;
+  buf << op << " array of length " << len << " with out-of-bounds index " << n;
+  error(buf);
+}
+
+inline item& arrayRead(array *a, Int n)  
+{
+  size_t len=checkArray(a);
+  bool cyclic=a->cyclic();
+  if(cyclic && len > 0) n=imod(n,len);
+  else if(n < 0 || n >= (Int) len) outOfBounds("reading",len,n);
+  return (*a)[(unsigned) n];
+}
+
+// Helper function to create deep arrays.
+static array* deepArray(Int depth, Int *dims)
+{
+  assert(depth > 0);
+  
+  if (depth == 1) {
+    return new array(dims[0]);
+  } else {
+    Int length = dims[0];
+    depth--; dims++;
+
+    array *a = new array(length);
+
+    for (Int index = 0; index < length; index++) {
+      (*a)[index] = deepArray(depth, dims);
+    }
+    return a;
+  }
+}
+
+namespace run {
+array *Identity(Int n)
+{
+  size_t N=(size_t) n;
+  array *c=new array(N);
+  for(size_t i=0; i < N; ++i) {
+    array *ci=new array(N);
+    (*c)[i]=ci;
+    for(size_t j=0; j < N; ++j)
+      (*ci)[j]=0.0;
+    (*ci)[i]=1.0;
+  }
+  return c;
+}
+}
+
+static const char *incommensurate="Incommensurate matrices";
+static const char *singular="Singular matrix";
+static size_t *pivot,*Row,*Col;
+
+namespace run {
+
+array *copyArray(array *a)
+{
+  size_t size=checkArray(a);
+  array *c=new array(size);
+  for(size_t i=0; i < size; i++) 
+    (*c)[i]=(*a)[i];
+  return c;
+}
+
+inline size_t checkdimension(const array *a, size_t dim)
+{
+  size_t size=checkArray(a);
+  if(dim && size != dim) {
+    ostringstream buf;
+    buf << "array of length " << dim << " expected";
+    error(buf);
+  }
+  return size;
+}
+
+double *copyArrayC(const array *a, size_t dim)
+{
+  size_t size=checkdimension(a,dim);
+  double *c=new double[size];
+  for(size_t i=0; i < size; i++) 
+    c[i]=read<double>(a,i);
+  return c;
+}
+
+triple *copyTripleArrayC(const array *a, size_t dim)
+{
+  size_t size=checkdimension(a,dim);
+  triple *c=new triple[size];
+  for(size_t i=0; i < size; i++) 
+    c[i]=read<triple>(a,i);
+  return c;
+}
+
+array *copyArray2(array *a)
+{
+  size_t size=checkArray(a);
+  array *c=new array(size);
+  for(size_t i=0; i < size; i++) {
+    array *ai=read<array*>(a,i);
+    size_t aisize=checkArray(ai);
+    array *ci=new array(aisize);
+    (*c)[i]=ci;
+    for(size_t j=0; j < aisize; j++) 
+      (*ci)[j]=(*ai)[j];
+  }
+  return c;
+}
+
+array *copyArray3(array *a)
+{
+  size_t size=checkArray(a);
+  array *c=new array(size);
+  for(size_t i=0; i < size; i++) {
+    array *ai=read<array*>(a,i);
+    size_t aisize=checkArray(ai);
+    array *ci=new array(aisize);
+    (*c)[i]=ci;
+    for(size_t j=0; j < aisize; j++) {
+      array *aij=read<array*>(ai,j);
+      size_t aijsize=checkArray(aij);
+      array *cij=new array(aijsize);
+      (*ci)[j]=cij;
+      for(size_t k=0; k < aijsize; k++) 
+        (*cij)[k]=(*aij)[k];
+    }
+  }
+  return c;
+}
+
+double *copyArray2C(const array *a, bool square, size_t dim2)
+{
+  size_t n=checkArray(a);
+  size_t m=(square || n == 0) ? n : checkArray(read<array*>(a,0));
+  if(n > 0 && dim2 && m != dim2) {
+    ostringstream buf;
+    buf << "second matrix dimension must be " << dim2;
+    error(buf);
+  }
+  
+  double *c=new double[n*m];
+  for(size_t i=0; i < n; i++) {
+    array *ai=read<array*>(a,i);
+    size_t aisize=checkArray(ai);
+    if(aisize == m) {
+      double *ci=c+i*m;
+      for(size_t j=0; j < m; j++) 
+        ci[j]=read<double>(ai,j);
+    } else
+      error(square ? "matrix must be square" : "matrix must be rectangular");
+  }
+  return c;
+}
+
+triple *copyTripleArray2C(const array *a, bool square, size_t dim2)
+{
+  size_t n=checkArray(a);
+  size_t m=(square || n == 0) ? n : checkArray(read<array*>(a,0));
+  if(n > 0 && dim2 && m != dim2) {
+    ostringstream buf;
+    buf << "second matrix dimension must be " << dim2;
+    error(buf);
+  }
+  
+  triple *c=new triple[n*m];
+  for(size_t i=0; i < n; i++) {
+    array *ai=read<array*>(a,i);
+    size_t aisize=checkArray(ai);
+    if(aisize == m) {
+      triple *ci=c+i*m;
+      for(size_t j=0; j < m; j++) 
+        ci[j]=read<triple>(ai,j);
+    } else
+      error(square ? "matrix must be square" : "matrix must be rectangular");
+  }
+  return c;
+}
+
+double *copyTripleArray2Components(array *a, bool square, size_t dim2)
+{
+  size_t n=checkArray(a);
+  size_t m=(square || n == 0) ? n : checkArray(read<array*>(a,0));
+  if(n > 0 && dim2 && m != dim2) {
+    ostringstream buf;
+    buf << "second matrix dimension must be " << dim2;
+    error(buf);
+  }
+  
+  size_t nm=n*m;
+  double *cx=new double[3*nm];
+  double *cy=cx+nm;
+  double *cz=cx+2*nm;
+  for(size_t i=0; i < n; i++) {
+    array *ai=read<array*>(a,i);
+    size_t aisize=checkArray(ai);
+    if(aisize == m) {
+      double *xi=cx+i*m;
+      double *yi=cy+i*m;
+      double *zi=cz+i*m;
+      for(size_t j=0; j < m; j++) {
+        triple v=read<triple>(ai,j);
+        xi[j]=v.getx();
+        yi[j]=v.gety();
+        zi[j]=v.getz();
+      }
+    } else
+      error(square ? "matrix must be square" : "matrix must be rectangular");
+  }
+  return cx;
+}
+
+triple operator *(const array& t, const triple& v)
+{
+  size_t n=checkArray(&t);
+  if(n != 4) error(incommensurate);
+  array *t0=read<array*>(t,0);
+  array *t1=read<array*>(t,1);
+  array *t2=read<array*>(t,2);
+  array *t3=read<array*>(t,3);
+  
+  if(checkArray(t0) != 4 || checkArray(t1) != 4 || 
+     checkArray(t2) != 4 || checkArray(t3) != 4)
+    error(incommensurate);
+
+  double x=v.getx();
+  double y=v.gety();
+  double z=v.getz();
+  
+  double f=read<real>(t3,0)*x+read<real>(t3,1)*y+read<real>(t3,2)*z+
+    read<real>(t3,3);
+  if(f == 0.0) run::dividebyzero();
+  f=1.0/f;
+  
+  return triple((read<real>(t0,0)*x+read<real>(t0,1)*y+read<real>(t0,2)*z+
+                 read<real>(t0,3))*f,
+                (read<real>(t1,0)*x+read<real>(t1,1)*y+read<real>(t1,2)*z+
+                 read<real>(t1,3))*f,
+                (read<real>(t2,0)*x+read<real>(t2,1)*y+read<real>(t2,2)*z+
+                 read<real>(t2,3))*f);
+}
+
+triple multshiftless(const array& t, const triple& v)
+{
+  size_t n=checkArray(&t);
+  if(n != 4) error(incommensurate);
+  array *t0=read<array*>(t,0);
+  array *t1=read<array*>(t,1);
+  array *t2=read<array*>(t,2);
+  array *t3=read<array*>(t,3);
+  
+  if(checkArray(t0) != 4 || checkArray(t1) != 4 || 
+     checkArray(t2) != 4 || checkArray(t3) != 4)
+    error(incommensurate);
+
+  double x=v.getx();
+  double y=v.gety();
+  double z=v.getz();
+  
+  double f=read<real>(t3,0)*x+read<real>(t3,1)*y+read<real>(t3,2)*z+
+    read<real>(t3,3);
+  if(f == 0.0) run::dividebyzero();
+  f=1.0/f;
+  
+  return triple((read<real>(t0,0)*x+read<real>(t0,1)*y+read<real>(t0,2)*z)*f,
+                (read<real>(t1,0)*x+read<real>(t1,1)*y+read<real>(t1,2)*z)*f,
+                (read<real>(t2,0)*x+read<real>(t2,1)*y+read<real>(t2,2)*z)*f);
+}
+
+double norm(double *a, size_t n) 
+{
+  if(n == 0) return 0.0;
+  double M=fabs(a[0]);
+  for(size_t i=1; i < n; ++i)
+    M=::max(M,fabs(a[i]));
+  return M;
+}
+
+double norm(triple *a, size_t n) 
+{
+  if(n == 0) return 0.0;
+  double M=a[0].abs2();
+  for(size_t i=1; i < n; ++i)
+    M=::max(M,a[i].abs2());
+  return sqrt(M);
+}
+
+}
+
+static inline void inverseAllocate(size_t n)
+{
+  pivot=new size_t[n];
+  Row=new size_t[n];
+  Col=new size_t[n];
+}
+
+static inline void inverseDeallocate()
+{
+  delete[] pivot;
+  delete[] Row;
+  delete[] Col;
+}
+
+callable *Func;
+stack *FuncStack;
+double wrapFunction(double x)
+{
+  FuncStack->push(x);
+  Func->call(FuncStack);
+  return pop<double>(FuncStack);
+}
+
+callable *compareFunc;
+bool compareFunction(const vm::item& i, const vm::item& j)
+{
+  FuncStack->push(i);
+  FuncStack->push(j);
+  compareFunc->call(FuncStack);
+  return pop<bool>(FuncStack);
+}
+
+void checkSquare(array *a) 
+{
+  size_t n=checkArray(a);
+  for(size_t i=0; i < n; i++)
+    if(checkArray(read<array*>(a,i)) != n)
+      error("matrix a must be square");
+}
+
+// Crout's algorithm for computing the LU decomposition of a square matrix.
+// cf. routine ludcmp (Press et al.,  Numerical Recipes, 1991).
+Int LUdecompose(double *a, size_t n, size_t* index, bool warn=true)
+{
+  double *vv=new double[n];
+  Int swap=1;
+  for(size_t i=0; i < n; ++i) {
+    double big=0.0;
+    double *ai=a+i*n;
+    for(size_t j=0; j < n; ++j) {
+      double temp=fabs(ai[j]);
+      if(temp > big) big=temp;
+    }
+    if(big == 0.0) {
+      delete[] vv;
+      if(warn) error(singular);
+      else return 0;
+    }
+    vv[i]=1.0/big;
+  }
+  for(size_t j=0; j < n; ++j) {
+    for(size_t i=0; i < j; ++i) {
+      double *ai=a+i*n;
+      double sum=ai[j];
+      for(size_t k=0; k < i; ++k) {
+        sum -= ai[k]*a[k*n+j];
+      }
+      ai[j]=sum;
+    }
+    double big=0.0;
+    size_t imax=j;
+    for(size_t i=j; i < n; ++i) {
+      double *ai=a+i*n;
+      double sum=ai[j];
+      for(size_t k=0; k < j; ++k)
+        sum -= ai[k]*a[k*n+j];
+      ai[j]=sum;
+      double temp=vv[i]*fabs(sum);
+      if(temp >= big) {
+        big=temp;
+        imax=i;
+      }
+    }
+    double *aj=a+j*n;
+    double *aimax=a+imax*n;
+    if(j != imax) {
+      for(size_t k=0; k < n; ++k) {
+        double temp=aimax[k];
+        aimax[k]=aj[k];
+        aj[k]=temp;
+      }
+      swap *= -1;
+      vv[imax]=vv[j];
+    }
+    if(index) 
+      index[j]=imax;
+    if(j != n) {
+      double denom=aj[j];
+      if(denom == 0.0) {
+        delete[] vv;
+        if(warn) error(singular);
+        else return 0;
+      }
+      for(size_t i=j+1; i < n; ++i)
+        a[i*n+j] /= denom;
+    }
+  }
+  delete[] vv;
+  return swap;
+}
+
+namespace run {
+void dividebyzero(size_t i)
+{
+  ostringstream buf;
+  if(i > 0) buf << "array element " << i << ": ";
+  buf << "Divide by zero";
+  error(buf);
+}
+  
+void integeroverflow(size_t i)
+{
+  ostringstream buf;
+  if(i > 0) buf << "array element " << i << ": ";
+  buf << "Integer overflow";
+  error(buf);
+}
+}
+
+// Autogenerated routines:
+
+
+
+namespace run {
+// Create an empty array.
+#line 474 "runarray.in"
+void emptyArray(stack *Stack)
+{
+#line 475 "runarray.in"
+  {Stack->push<array*>(new array(0)); return;}
+}
+
+// Create a new array (technically a vector).
+// This array will be multidimensional.  First the number of dimensions
+// is popped off the stack, followed by each dimension in reverse order.
+// The array itself is technically a one dimensional array of one
+// dimension arrays and so on.
+#line 484 "runarray.in"
+void newDeepArray(stack *Stack)
+{
+  Int depth=vm::pop<Int>(Stack);
+#line 485 "runarray.in"
+  assert(depth > 0);
+
+  Int *dims = new Int[depth];
+
+  for (Int index = depth-1; index >= 0; index--) {
+    Int i=pop<Int>(Stack);
+    if(i < 0) error("cannot create a negative length array");
+    dims[index]=i;
+  }
+
+  array *a=deepArray(depth, dims);
+  delete[] dims;
+  {Stack->push<array*>(a); return;}
+}
+
+// Creates an array with elements already specified.  First, the number
+// of elements is popped off the stack, followed by each element in
+// reverse order.
+#line 504 "runarray.in"
+void newInitializedArray(stack *Stack)
+{
+  Int n=vm::pop<Int>(Stack);
+#line 505 "runarray.in"
+  assert(n >= 0);
+
+  array *a = new array(n);
+
+  for (Int index = n-1; index >= 0; index--)
+    (*a)[index] = pop(Stack);
+
+  {Stack->push<array*>(a); return;}
+}
+
+// Similar to newInitializedArray, but after the n elements, append another
+// array to it.
+#line 518 "runarray.in"
+void newAppendedArray(stack *Stack)
+{
+  Int n=vm::pop<Int>(Stack);
+  array* tail=vm::pop<array*>(Stack);
+#line 519 "runarray.in"
+  assert(n >= 0);
+
+  array *a = new array(n);
+
+  for (Int index = n-1; index >= 0; index--)
+    (*a)[index] = pop(Stack);
+  
+  copy(tail->begin(), tail->end(), back_inserter(*a));
+
+  {Stack->push<array*>(a); return;}
+}
+
+// The function T[] array(int n, T value, int depth=0) produces a array of n
+// copies of x, where each copy is copied up to depth.
+#line 534 "runarray.in"
+void newDuplicateArray(stack *Stack)
+{
+  Int depth=vm::pop<Int>(Stack,Int_MAX);
+  item value=vm::pop(Stack);
+  Int n=vm::pop<Int>(Stack);
+#line 535 "runarray.in"
+  if(n < 0) error("cannot create a negative length array");
+  if(depth < 0) error("cannot copy to a negative depth");
+
+  {Stack->push<array*>(new array(n, value, depth)); return;}
+}
+
+// Read an element from an array. Checks for initialization & bounds.
+#line 543 "runarray.in"
+void arrayRead(stack *Stack)
+{
+  Int n=vm::pop<Int>(Stack);
+  array * a=vm::pop<array *>(Stack);
+#line 544 "runarray.in"
+  item& i=arrayRead(a,n);
+  if (i.empty()) {
+    ostringstream buf;
+    buf << "read uninitialized value from array at index " << n;
+    error(buf);
+  }
+  {Stack->push(i); return;}
+}
+
+// Slice a substring from an array.
+#line 555 "runarray.in"
+void arraySliceRead(stack *Stack)
+{
+  Int right=vm::pop<Int>(Stack);
+  Int left=vm::pop<Int>(Stack);
+  array * a=vm::pop<array *>(Stack);
+#line 556 "runarray.in"
+  checkArray(a);
+  {Stack->push(a->slice(left, right)); return;}
+}
+
+// Slice a substring from an array.  This implements the cases a[i:] and a[:]
+// where the endpoint is not given, and assumed to be the length of the array.
+#line 563 "runarray.in"
+void arraySliceReadToEnd(stack *Stack)
+{
+  Int left=vm::pop<Int>(Stack);
+  array * a=vm::pop<array *>(Stack);
+#line 564 "runarray.in"
+  size_t len=checkArray(a);
+  {Stack->push(a->slice(left, (Int)len)); return;}
+}
+
+// Read an element from an array of arrays. Check bounds and initialize
+// as necessary.
+#line 571 "runarray.in"
+void arrayArrayRead(stack *Stack)
+{
+  Int n=vm::pop<Int>(Stack);
+  array * a=vm::pop<array *>(Stack);
+#line 572 "runarray.in"
+  item& i=arrayRead(a,n);
+  if (i.empty()) i=new array(0);
+  {Stack->push(i); return;}
+}
+
+// Write an element to an array.  Increase size if necessary.
+#line 579 "runarray.in"
+void arrayWrite(stack *Stack)
+{
+  Int n=vm::pop<Int>(Stack);
+  array * a=vm::pop<array *>(Stack);
+  item value=vm::pop(Stack);
+#line 580 "runarray.in"
+  size_t len=checkArray(a);
+  bool cyclic=a->cyclic();
+  if(cyclic && len > 0) n=imod(n,len);
+  else {
+    if(cyclic) outOfBounds("writing cyclic",len,n);
+    if(n < 0) outOfBounds("writing",len,n);
+    if(len <= (size_t) n)
+      a->resize(n+1);
+  }
+  (*a)[n] = value;
+  {Stack->push(value); return;}
+}
+
+#line 594 "runarray.in"
+void arraySliceWrite(stack *Stack)
+{
+  Int right=vm::pop<Int>(Stack);
+  Int left=vm::pop<Int>(Stack);
+  array * dest=vm::pop<array *>(Stack);
+  array * src=vm::pop<array *>(Stack);
+#line 595 "runarray.in"
+  checkArray(src);
+  checkArray(dest);
+  dest->setSlice(left, right, src);
+  {Stack->push<array*>(src); return;}
+}
+
+#line 602 "runarray.in"
+void arraySliceWriteToEnd(stack *Stack)
+{
+  Int left=vm::pop<Int>(Stack);
+  array * dest=vm::pop<array *>(Stack);
+  array * src=vm::pop<array *>(Stack);
+#line 603 "runarray.in"
+  checkArray(src);
+  size_t len=checkArray(dest);
+  dest->setSlice(left, (Int) len, src);
+  {Stack->push<array*>(src); return;}
+}
+
+// Returns the length of an array.
+#line 611 "runarray.in"
+void arrayLength(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+#line 612 "runarray.in"
+  {Stack->push<Int>((Int) checkArray(a)); return;}
+}
+
+// Returns an array of integers representing the keys of the array.
+#line 617 "runarray.in"
+void arrayKeys(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+#line 618 "runarray.in"
+  size_t size=checkArray(a);
+
+  array *keys=new array();
+  for (size_t i=0; i<size; ++i) {
+    item& cell = (*a)[i];
+    if (!cell.empty())
+      keys->push((Int)i);
+  }
+
+  {Stack->push<array*>(keys); return;}
+}
+
+// Return the cyclic flag for an array.
+#line 632 "runarray.in"
+void arrayCyclicFlag(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+#line 633 "runarray.in"
+  checkArray(a);
+  {Stack->push<bool>(a->cyclic()); return;}
+}
+
+#line 638 "runarray.in"
+void arraySetCyclicFlag(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+  bool b=vm::pop<bool>(Stack);
+#line 639 "runarray.in"
+  checkArray(a);
+  a->cyclic(b);
+  {Stack->push<bool>(b); return;}
+}
+
+// Check to see if an array element is initialized.
+#line 646 "runarray.in"
+void arrayInitializedHelper(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+  Int n=vm::pop<Int>(Stack);
+#line 647 "runarray.in"
+  size_t len=checkArray(a);
+  bool cyclic=a->cyclic();
+  if(cyclic && len > 0) n=imod(n,len);
+  else if(n < 0 || n >= (Int) len) {Stack->push<bool>(false); return;}
+  item&i=(*a)[(unsigned) n];
+  {Stack->push<bool>(!i.empty()); return;}
+}
+
+// Returns the initialize method for an array.
+#line 657 "runarray.in"
+void arrayInitialized(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+#line 658 "runarray.in"
+  {Stack->push<callable*>(new thunk(new bfunc(arrayInitializedHelper),a)); return;}
+}
+
+// The helper function for the cyclic method that sets the cyclic flag.
+#line 663 "runarray.in"
+void arrayCyclicHelper(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+  bool b=vm::pop<bool>(Stack);
+#line 664 "runarray.in"
+  checkArray(a);
+  a->cyclic(b);
+}
+
+// Set the cyclic flag for an array.
+#line 670 "runarray.in"
+void arrayCyclic(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+#line 671 "runarray.in"
+  {Stack->push<callable*>(new thunk(new bfunc(arrayCyclicHelper),a)); return;}
+}
+
+// The helper function for the push method that does the actual operation.
+#line 676 "runarray.in"
+void arrayPushHelper(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+  item x=vm::pop(Stack);
+#line 677 "runarray.in"
+  checkArray(a);
+  a->push(x);
+  {Stack->push(x); return;}
+}
+
+// Returns the push method for an array.
+#line 684 "runarray.in"
+void arrayPush(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+#line 685 "runarray.in"
+  {Stack->push<callable*>(new thunk(new bfunc(arrayPushHelper),a)); return;}
+}
+
+// The helper function for the append method that appends b to a.
+#line 690 "runarray.in"
+void arrayAppendHelper(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+  array * b=vm::pop<array *>(Stack);
+#line 691 "runarray.in"
+  checkArray(a);
+  size_t size=checkArray(b);
+  for(size_t i=0; i < size; i++)
+    a->push((*b)[i]);
+}
+
+// Returns the append method for an array.
+#line 699 "runarray.in"
+void arrayAppend(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+#line 700 "runarray.in"
+  {Stack->push<callable*>(new thunk(new bfunc(arrayAppendHelper),a)); return;}
+}
+
+// The helper function for the pop method.
+#line 705 "runarray.in"
+void arrayPopHelper(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+#line 706 "runarray.in"
+  size_t asize=checkArray(a);
+  if(asize == 0) 
+    error("cannot pop element from empty array");
+  {Stack->push(a->pop()); return;}
+}
+
+// Returns the pop method for an array.
+#line 714 "runarray.in"
+void arrayPop(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+#line 715 "runarray.in"
+  {Stack->push<callable*>(new thunk(new bfunc(arrayPopHelper),a)); return;}
+}
+
+// The helper function for the insert method.
+#line 720 "runarray.in"
+void arrayInsertHelper(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+  array * x=vm::pop<array *>(Stack);
+  Int i=vm::pop<Int>(Stack);
+#line 721 "runarray.in"
+  size_t asize=checkArray(a);
+  checkArray(x);
+  if(a->cyclic() && asize > 0) i=imod(i,asize);
+  if(i < 0 || i > (Int) asize) 
+    outOfBounds("inserting",asize,i);
+  (*a).insert((*a).begin()+i,(*x).begin(),(*x).end());
+}
+
+// Returns the insert method for an array.
+#line 731 "runarray.in"
+void arrayInsert(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+#line 732 "runarray.in"
+  {Stack->push<callable*>(new thunk(new bfunc(arrayInsertHelper),a)); return;}
+}
+
+// Returns the delete method for an array.
+#line 737 "runarray.in"
+void arrayDelete(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+#line 738 "runarray.in"
+  {Stack->push<callable*>(new thunk(new bfunc(arrayDeleteHelper),a)); return;}
+}
+
+#line 742 "runarray.in"
+void arrayAlias(stack *Stack)
+{
+  array * b=vm::pop<array *>(Stack);
+  array * a=vm::pop<array *>(Stack);
+#line 743 "runarray.in"
+  {Stack->push<bool>(a==b); return;}
+}
+
+// Return array formed by indexing array a with elements of integer array b
+#line 748 "runarray.in"
+void arrayIntArray(stack *Stack)
+{
+  array * b=vm::pop<array *>(Stack);
+  array * a=vm::pop<array *>(Stack);
+#line 749 "runarray.in"
+  size_t asize=checkArray(a);
+  size_t bsize=checkArray(b);
+  array *r=new array(bsize);
+  bool cyclic=a->cyclic();
+  for(size_t i=0; i < bsize; i++) {
+    Int index=read<Int>(b,i);
+    if(cyclic && asize > 0) index=imod(index,asize);
+    else
+      if(index < 0 || index >= (Int) asize)
+        outOfBounds("reading",asize,index);
+    (*r)[i]=(*a)[index];
+  }
+  {Stack->push<array*>(r); return;}
+}
+
+// returns the complement of the integer array a in {0,2,...,n-1},
+// so that b[complement(a,b.length)] yields the complement of b[a].
+#line 767 "runarray.in"
+// Intarray* complement(Intarray *a, Int n);
+void gen_runarray31(stack *Stack)
+{
+  Int n=vm::pop<Int>(Stack);
+  Intarray * a=vm::pop<Intarray *>(Stack);
+#line 768 "runarray.in"
+  size_t asize=checkArray(a);
+  array *r=new array(0);
+  bool *keep=new bool[n];
+  for(Int i=0; i < n; ++i) keep[i]=true;
+  for(size_t i=0; i < asize; ++i) {
+    Int j=read<Int>(a,i);
+    if(j >= 0 && j < n) keep[j]=false;
+  }
+  for(Int i=0; i < n; i++)
+    if(keep[i]) r->push(i);
+  
+  delete[] keep;
+  {Stack->push<Intarray*>(r); return;}
+}
+
+// Generate the sequence {f(i) : i=0,1,...n-1} given a function f and integer n
+#line 785 "runarray.in"
+void arraySequence(stack *Stack)
+{
+  Int n=vm::pop<Int>(Stack);
+  callable * f=vm::pop<callable *>(Stack);
+#line 786 "runarray.in"
+  if(n < 0) n=0;
+  array *a=new array(n);
+  for(Int i=0; i < n; ++i) {
+    Stack->push(i);
+    f->call(Stack);
+    (*a)[i]=pop(Stack);
+  }
+  {Stack->push<Intarray*>(a); return;}
+}
+
+// Return the array {0,1,...n-1}
+#line 798 "runarray.in"
+// Intarray* sequence(Int n);
+void gen_runarray33(stack *Stack)
+{
+  Int n=vm::pop<Int>(Stack);
+#line 799 "runarray.in"
+  if(n < 0) n=0;
+  array *a=new array(n);
+  for(Int i=0; i < n; ++i) {
+    (*a)[i]=i;
+  }
+  {Stack->push<Intarray*>(a); return;}
+}
+
+// Apply a function to each element of an array
+#line 809 "runarray.in"
+void arrayFunction(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+  callable * f=vm::pop<callable *>(Stack);
+#line 810 "runarray.in"
+  size_t size=checkArray(a);
+  array *b=new array(size);
+  for(size_t i=0; i < size; ++i) {
+    Stack->push((*a)[i]);
+    f->call(Stack);
+    (*b)[i]=pop(Stack);
+  }
+  {Stack->push<array*>(b); return;}
+}
+
+#line 821 "runarray.in"
+void arraySort(stack *Stack)
+{
+  callable * f=vm::pop<callable *>(Stack);
+  array * a=vm::pop<array *>(Stack);
+#line 822 "runarray.in"
+  array *c=copyArray(a);
+  compareFunc=f;
+  FuncStack=Stack;
+  stable_sort(c->begin(),c->end(),compareFunction);
+  {Stack->push<array*>(c); return;}
+}
+
+#line 830 "runarray.in"
+// bool all(boolarray *a);
+void gen_runarray36(stack *Stack)
+{
+  boolarray * a=vm::pop<boolarray *>(Stack);
+#line 831 "runarray.in"
+  size_t size=checkArray(a);
+  bool c=true;
+  for(size_t i=0; i < size; i++)
+    if(!get<bool>((*a)[i])) {c=false; break;}
+  {Stack->push<bool>(c); return;}
+}
+
+#line 839 "runarray.in"
+// boolarray* !(boolarray* a);
+void gen_runarray37(stack *Stack)
+{
+  boolarray* a=vm::pop<boolarray*>(Stack);
+#line 840 "runarray.in"
+  size_t size=checkArray(a);
+  array *c=new array(size);
+  for(size_t i=0; i < size; i++)
+    (*c)[i]=!read<bool>(a,i);
+  {Stack->push<boolarray*>(c); return;}
+}
+
+#line 848 "runarray.in"
+// Int sum(boolarray *a);
+void gen_runarray38(stack *Stack)
+{
+  boolarray * a=vm::pop<boolarray *>(Stack);
+#line 849 "runarray.in"
+  size_t size=checkArray(a);
+  Int sum=0;
+  for(size_t i=0; i < size; i++)
+    sum += read<bool>(a,i) ? 1 : 0;
+  {Stack->push<Int>(sum); return;}
+}
+
+#line 857 "runarray.in"
+void arrayCopy(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+#line 858 "runarray.in"
+  {Stack->push<array*>(copyArray(a)); return;}
+}
+
+#line 862 "runarray.in"
+void arrayConcat(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+#line 863 "runarray.in"
+  // a is an array of arrays to be concatenated together.
+  // The signature is
+  //   T[] concat(... T[][] a);
+
+  size_t numArgs=checkArray(a);
+  size_t resultSize=0;
+  for (size_t i=0; i < numArgs; ++i) {
+    resultSize += checkArray(a->read<array *>(i));
+  }
+
+  array *result=new array(resultSize);
+
+  size_t ri=0;
+  for (size_t i=0; i < numArgs; ++i) {
+    array *arg=a->read<array *>(i);
+    size_t size=checkArray(arg);
+
+    for (size_t j=0; j < size; ++j) {
+      (*result)[ri]=(*arg)[j];
+      ++ri;
+    }
+  }
+
+  {Stack->push<array*>(result); return;}
+}
+
+#line 890 "runarray.in"
+void array2Copy(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+#line 891 "runarray.in"
+  {Stack->push<array*>(copyArray2(a)); return;}
+}
+
+#line 895 "runarray.in"
+void array3Copy(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+#line 896 "runarray.in"
+  {Stack->push<array*>(copyArray3(a)); return;}
+}
+
+#line 900 "runarray.in"
+void array2Transpose(stack *Stack)
+{
+  array * a=vm::pop<array *>(Stack);
+#line 901 "runarray.in"
+  size_t asize=checkArray(a);
+  array *c=new array(0);
+  for(size_t i=0; i < asize; i++) {
+    size_t ip=i+1;
+    array *ai=read<array*>(a,i);
+    size_t aisize=checkArray(ai);
+    size_t csize=checkArray(c);
+    if(csize < aisize) {
+      c->resize(aisize);
+      for(size_t j=csize; j < aisize; j++) {
+        (*c)[j]=new array(ip);
+      }
+    }
+    for(size_t j=0; j < aisize; j++) {
+      array *cj=read<array*>(c,j);
+      if(checkArray(cj) < ip) cj->resize(ip);
+      (*cj)[i]=(*ai)[j];
+    }
+  }
+  {Stack->push<array*>(c); return;}
+}
+
+// a is a rectangular 3D array; perm is an Int array indicating the type of
+// permutation  (021 or 120, etc; original is 012).
+// Transpose by sending respective members to the permutated locations:
+// return the array obtained by putting a[i][j][k] into position perm{ijk}. 
+#line 928 "runarray.in"
+void array3Transpose(stack *Stack)
+{
+  array * perm=vm::pop<array *>(Stack);
+  array * a=vm::pop<array *>(Stack);
+#line 929 "runarray.in"
+  const size_t DIM=3;
+
+  if(checkArray(perm) != DIM) {
+    ostringstream buf;
+    buf << "permutation array must have length " << DIM;
+    error(buf);
+  }
+  
+  size_t* size=new size_t[DIM];
+  for(size_t i=0; i < DIM; ++i) size[i]=DIM;
+  
+  for(size_t i=0; i < DIM; ++i) {
+    Int p=read<Int>(perm,i);
+    size_t P=(size_t) p;
+    if(p < 0 || P >= DIM) {
+      ostringstream buf;
+      buf << "permutation index out of range: " << p;
+      error(buf);
+    }
+    size[P]=P;
+  }
+  
+  for(size_t i=0; i < DIM; ++i)
+    if(size[i] == DIM) error("permutation indices must be distinct");
+  
+  static const char *rectangular=
+    "3D transpose implemented for rectangular matrices only";
+  
+  size_t isize=size[0]=checkArray(a);
+  array *a0=read<array*>(a,0);
+  size[1]=checkArray(a0);
+  array *a00=read<array*>(a0,0);
+  size[2]=checkArray(a00);
+  for(size_t i=0; i < isize; i++) {
+    array *ai=read<array*>(a,i);
+    size_t jsize=checkArray(ai);
+    if(jsize != size[1]) error(rectangular);
+    for(size_t j=0; j < jsize; j++) {
+      array *aij=read<array*>(ai,j);
+      if(checkArray(aij) != size[2]) error(rectangular);
+    }
+  }
+  
+  size_t perm0=(size_t) read<Int>(perm,0);
+  size_t perm1=(size_t) read<Int>(perm,1);
+  size_t perm2=(size_t) read<Int>(perm,2);
+  
+  size_t sizep0=size[perm0];
+  size_t sizep1=size[perm1];
+  size_t sizep2=size[perm2];
+  
+  array *c=new array(sizep0);
+  for(size_t i=0; i < sizep0; ++i) {
+    array *ci=new array(sizep1);
+    (*c)[i]=ci;
+    for(size_t j=0; j < sizep1; ++j) {
+      array *cij=new array(sizep2);
+      (*ci)[j]=cij;
+    }
+  }
+  
+  size_t* i=new size_t[DIM];
+  
+  for(i[0]=0; i[0] < size[0]; ++i[0]) {
+    array *a0=read<array*>(a,i[0]);
+    for(i[1]=0; i[1] < size[1]; ++i[1]) {
+      array *a1=read<array*>(a0,i[1]);
+      for(i[2]=0; i[2] < size[2]; ++i[2]) {
+        array *c0=read<array*>(c,i[perm0]);
+        array *c1=read<array*>(c0,i[perm1]);
+        (*c1)[i[perm2]]=read<real>(a1,i[2]);
+      }
+    }
+  }
+  
+  delete [] i;  
+  delete [] size;  
+
+  {Stack->push<array*>(c); return;}
+}
+
+// In a boolean array, find the index of the nth true value or -1 if not found
+// If n is negative, search backwards.
+#line 1013 "runarray.in"
+// Int find(boolarray *a, Int n=1);
+void gen_runarray45(stack *Stack)
+{
+  Int n=vm::pop<Int>(Stack,1);
+  boolarray * a=vm::pop<boolarray *>(Stack);
+#line 1014 "runarray.in"
+  size_t size=checkArray(a);
+  Int j=-1;
+  if(n > 0)
+    for(size_t i=0; i < size; i++)
+      if(read<bool>(a,i)) {
+        n--; if(n == 0) {j=(Int) i; break;}
+      }
+  if(n < 0)
+    for(size_t i=size; i > 0;)
+      if(read<bool>(a,--i)) {
+        n++; if(n == 0) {j=(Int) i; break;}
+      }
+  {Stack->push<Int>(j); return;}
+}
+
+// construct vector obtained by replacing those elements of b for which the
+// corresponding elements of a are false by the corresponding element of c.
+#line 1032 "runarray.in"
+void arrayConditional(stack *Stack)
+{
+  array * c=vm::pop<array *>(Stack);
+  array * b=vm::pop<array *>(Stack);
+  array * a=vm::pop<array *>(Stack);
+#line 1033 "runarray.in"
+  size_t size=checkArray(a);
+  array *r=new array(size);
+  if(b && c) {
+    checkArrays(a,b);
+    checkArrays(b,c);
+    for(size_t i=0; i < size; i++)
+      (*r)[i]=read<bool>(a,i) ? (*b)[i] : (*c)[i];
+  } else {
+    r->clear();
+    if(b) {
+      checkArrays(a,b);
+      for(size_t i=0; i < size; i++)
+        if(read<bool>(a,i)) r->push((*b)[i]);
+    } else if(c) {
+      checkArrays(a,c);
+      for(size_t i=0; i < size; i++)
+        if(!read<bool>(a,i)) r->push((*c)[i]);
+    }
+  }
+  {Stack->push<array*>(r); return;}
+}
+
+// Return an n x n identity matrix.
+#line 1057 "runarray.in"
+// realarray2* identity(Int n);
+void gen_runarray47(stack *Stack)
+{
+  Int n=vm::pop<Int>(Stack);
+#line 1058 "runarray.in"
+  {Stack->push<realarray2*>(Identity(n)); return;}
+}
+
+// Return the diagonal matrix with diagonal entries given by a.
+#line 1063 "runarray.in"
+void diagonal(stack *Stack)
+{
+  realarray * a=vm::pop<realarray *>(Stack);
+#line 1064 "runarray.in"
+  size_t n=checkArray(a);
+  array *c=new array(n);
+  for(size_t i=0; i < n; ++i) {
+    array *ci=new array(n);
+    (*c)[i]=ci;
+    for(size_t j=0; j < i; ++j)
+      (*ci)[j]=0.0;
+    (*ci)[i]=read<real>(a,i);
+    for(size_t j=i+1; j < n; ++j)
+      (*ci)[j]=0.0;
+  }
+  {Stack->push<realarray2*>(c); return;}
+}
+
+// Return the inverse of an n x n matrix a using Gauss-Jordan elimination.
+#line 1080 "runarray.in"
+// realarray2* inverse(realarray2 *a);
+void gen_runarray49(stack *Stack)
+{
+  realarray2 * a=vm::pop<realarray2 *>(Stack);
+#line 1081 "runarray.in"
+  a=copyArray2(a);
+  size_t n=checkArray(a);
+  checkSquare(a);
+  
+  inverseAllocate(n);
+  
+  for(size_t i=0; i < n; i++)
+    pivot[i]=0;
+ 
+  size_t col=0, row=0;
+  // This is the main loop over the columns to be reduced.
+  for(size_t i=0; i < n; i++) {
+    real big=0.0;
+    // This is the outer loop of the search for a pivot element.
+    for(size_t j=0; j < n; j++) {
+      array *aj=read<array*>(a,j);
+      if(pivot[j] != 1) {
+        for(size_t k=0; k < n; k++) {
+          if(pivot[k] == 0) {
+            real temp=fabs(read<real>(aj,k));
+            if(temp >= big) {
+              big=temp;
+              row=j;
+              col=k;
+            }
+          } else if(pivot[k] > 1) {
+            inverseDeallocate();
+            error(singular);
+          }
+        }
+      }
+    }
+    ++(pivot[col]);
+    
+    // Interchange rows, if needed, to put the pivot element on the diagonal.
+    array *acol=read<array*>(a,col);
+    if(row != col) {
+      array *arow=read<array*>(a,row);
+      for(size_t l=0; l < n; l++) {
+        real temp=read<real>(arow,l);
+        (*arow)[l]=read<real>(acol,l);
+        (*acol)[l]=temp;
+      }
+    }
+    
+    Row[i]=row; 
+    Col[i]=col;
+
+    // Divide the pivot row by the pivot element.
+    real denom=read<real>(acol,col);
+    if(denom == 0.0) {
+      inverseDeallocate();
+      error(singular);
+    }
+    real pivinv=1.0/denom;
+    (*acol)[col]=1.0;
+    for(size_t l=0; l < n; l++) 
+      (*acol)[l]=read<real>(acol,l)*pivinv;
+    
+    // Reduce all rows except for the pivoted one.
+    for(size_t k=0; k < n; k++) {
+      if(k != col) { 
+        array *ak=read<array*>(a,k);
+        real akcol=read<real>(ak,col);
+        (*ak)[col]=0.0;
+        for(size_t l=0; l < n; l++)
+          (*ak)[l]=read<real>(ak,l)-read<real>(acol,l)*akcol;
+      }
+    }
+  }
+  
+  // Unscramble the inverse matrix in view of the column interchanges.
+  for(size_t l=n; l > 0;) {
+    l--;
+    size_t r=Row[l];
+    size_t c=Col[l];
+    if(r != c) {
+      for(size_t k=0; k < n; k++) {
+        array *ak=read<array*>(a,k);
+        real temp=read<real>(ak,r);
+        (*ak)[r]=read<real>(ak,c);
+        (*ak)[c]=temp;
+      }
+    }
+  }
+  inverseDeallocate();
+  {Stack->push<realarray2*>(a); return;}
+}
+
+// Solve the linear equation ax=b by LU decomposition, returning the
+// solution x, where a is an n x n matrix and b is an array of length n.
+// If no solution exists, return an empty array.
+#line 1174 "runarray.in"
+// realarray* solve(realarray2 *a, realarray *b, bool warn=true);
+void gen_runarray50(stack *Stack)
+{
+  bool warn=vm::pop<bool>(Stack,true);
+  realarray * b=vm::pop<realarray *>(Stack);
+  realarray2 * a=vm::pop<realarray2 *>(Stack);
+#line 1175 "runarray.in"
+  size_t n=checkArray(a);
+  
+  if(n == 0) {Stack->push<realarray*>(new array(0)); return;}
+  
+  size_t m=checkArray(b);
+  if(m != n) error(incommensurate);
+  
+  real *A=copyArray2C(a);
+  size_t *index=new size_t[n];
+  
+  if(LUdecompose(A,n,index,warn) == 0)
+    {Stack->push<realarray*>(new array(0)); return;}
+
+  array *x=new array(n);
+  
+  real *B=copyArrayC(b);
+  
+  for(size_t i=0; i < n; ++i) {
+    size_t ip=index[i];
+    real sum=B[ip];
+    B[ip]=B[i];
+    real *Ai=A+i*n;
+    for(size_t j=0; j < i; ++j)
+      sum -= Ai[j]*B[j];
+    B[i]=sum;
+  }
+  
+  for(size_t i=n; i > 0;) {
+    --i;
+    real sum=B[i];
+    real *Ai=A+i*n;
+    for(size_t j=i+1; j < n; ++j)
+      sum -= Ai[j]*B[j];
+    B[i]=sum/Ai[i];
+  }
+  
+  for(size_t i=0; i < n; ++i)
+    (*x)[i]=B[i];
+
+  delete[] index;
+  delete[] B;
+  delete[] A;
+  
+  {Stack->push<realarray*>(x); return;}
+}
+
+// Solve the linear equation ax=b by LU decomposition, returning the
+// solution x, where a is an n x n matrix and b is an n x m matrix.
+// If no solution exists, return an empty array.
+#line 1225 "runarray.in"
+// realarray2* solve(realarray2 *a, realarray2 *b, bool warn=true);
+void gen_runarray51(stack *Stack)
+{
+  bool warn=vm::pop<bool>(Stack,true);
+  realarray2 * b=vm::pop<realarray2 *>(Stack);
+  realarray2 * a=vm::pop<realarray2 *>(Stack);
+#line 1226 "runarray.in"
+  size_t n=checkArray(a);
+  
+  if(n == 0) {Stack->push<realarray2*>(new array(0)); return;}
+  
+  if(checkArray(b) != n) error(incommensurate);
+  size_t m=checkArray(read<array*>(b,0));
+  
+  real *A=copyArray2C(a);
+  real *B=copyArray2C(b,false);
+  
+  size_t *index=new size_t[n];
+  
+  if(LUdecompose(A,n,index,warn) == 0)
+    {Stack->push<realarray2*>(new array(0)); return;}
+
+  array *x=new array(n);
+  
+  for(size_t i=0; i < n; ++i) {
+    real *Ai=A+i*n;
+    real *Bi=B+i*m;
+    real *Bip=B+index[i]*m;
+    for(size_t k=0; k < m; ++k) {
+      real sum=Bip[k];
+      Bip[k]=Bi[k];
+      size_t jk=k;
+      for(size_t j=0; j < i; ++j, jk += m)
+        sum -= Ai[j]*B[jk];
+      Bi[k]=sum;
+    }
+  }
+  
+  for(size_t i=n; i > 0;) {
+    --i;
+    real *Ai=A+i*n;
+    real *Bi=B+i*m;
+    for(size_t k=0; k < m; ++k) {
+      real sum=Bi[k];
+      size_t jk=(i+1)*m+k;
+      for(size_t j=i+1; j < n; ++j, jk += m)
+        sum -= Ai[j]*B[jk];
+      Bi[k]=sum/Ai[i];
+    }
+  }
+  
+  for(size_t i=0; i < n; ++i) {
+    real *Bi=B+i*m;
+    array *xi=new array(m);
+    (*x)[i]=xi;
+    for(size_t j=0; j < m; ++j)
+      (*xi)[j]=Bi[j];
+  }
+    
+  delete[] index;
+  delete[] B;
+  delete[] A;
+  
+  {Stack->push<realarray2*>(x); return;}
+}
+
+// Compute the determinant of an n x n matrix.
+#line 1287 "runarray.in"
+// real determinant(realarray2 *a);
+void gen_runarray52(stack *Stack)
+{
+  realarray2 * a=vm::pop<realarray2 *>(Stack);
+#line 1288 "runarray.in"
+  real *A=copyArray2C(a);
+  size_t n=checkArray(a);
+  
+  real det=LUdecompose(A,n,NULL,false);
+  size_t n1=n+1;
+  for(size_t i=0; i < n; ++i)
+    det *= A[i*n1];
+  
+  delete[] A;
+  
+  {Stack->push<real>(det); return;}
+}
+
+#line 1302 "runarray.in"
+// realarray* *(realarray2 *a, realarray *b);
+void gen_runarray53(stack *Stack)
+{
+  realarray * b=vm::pop<realarray *>(Stack);
+  realarray2 * a=vm::pop<realarray2 *>(Stack);
+#line 1303 "runarray.in"
+  size_t n=checkArray(a);
+  size_t m=checkArray(b);
+  array *c=new array(n);
+  real *B=copyArrayC(b);
+  for(size_t i=0; i < n; ++i) {
+    array *ai=read<array*>(a,i);
+    if(checkArray(ai) != m) error(incommensurate);
+    real sum=0.0;
+    for(size_t j=0; j < m; ++j)
+      sum += read<real>(ai,j)*B[j];
+    (*c)[i]=sum;
+  }
+  delete[] B;
+  {Stack->push<realarray*>(c); return;}
+}
+
+#line 1320 "runarray.in"
+// realarray* *(realarray *a, realarray2 *b);
+void gen_runarray54(stack *Stack)
+{
+  realarray2 * b=vm::pop<realarray2 *>(Stack);
+  realarray * a=vm::pop<realarray *>(Stack);
+#line 1321 "runarray.in"
+  size_t n=checkArray(a);
+  if(n != checkArray(b)) error(incommensurate);
+  real *A=copyArrayC(a);
+
+  array **B=new array*[n];
+  array *bk=read<array *>(b,0);
+  B[0]=bk;
+  size_t m=bk->size();
+  for(size_t k=1; k < n; k++) {
+    array *bk=read<array *>(b,k);
+    if(bk->size() != m) error(incommensurate);
+    B[k]=bk;
+  }
+  array *c=new array(m);
+
+  for(size_t i=0; i < m; ++i) {
+    real sum=0.0;
+    for(size_t k=0; k < n; ++k)
+      sum += A[k]*read<real>(B[k],i);
+    (*c)[i]=sum;
+  }
+  delete[] B;
+  delete[] A;
+  {Stack->push<realarray*>(c); return;}
+}
+
+#line 1348 "runarray.in"
+// realarray2* *(realarray2 *a, realarray2 *b);
+void gen_runarray55(stack *Stack)
+{
+  realarray2 * b=vm::pop<realarray2 *>(Stack);
+  realarray2 * a=vm::pop<realarray2 *>(Stack);
+#line 1349 "runarray.in"
+  size_t n=checkArray(a);
+  
+  size_t nb=checkArray(b);
+  size_t na0=n == 0 ? 0 : checkArray(read<array*>(a,0));
+  if(na0 != nb) 
+    error(incommensurate);
+  
+  size_t nb0=nb == 0 ? 0 : checkArray(read<array*>(b,0));
+    
+  array *c=new array(n);
+
+  real *A=copyArray2C(a,false);
+  real *B=copyArray2C(b,false);
+
+  for(size_t i=0; i < n; ++i) {
+    real *Ai=A+i*nb;
+    array *ci=new array(nb0);
+    (*c)[i]=ci;
+    for(size_t j=0; j < nb0; ++j) {
+      real sum=0.0;
+      size_t kj=j;
+      for(size_t k=0; k < nb; ++k, kj += nb0)
+        sum += Ai[k]*B[kj];
+      (*ci)[j]=sum;
+    }
+  }
+  
+  delete[] B;
+  delete[] A;
+  
+  {Stack->push<realarray2*>(c); return;}
+}
+
+#line 1383 "runarray.in"
+// triple *(realarray2 *t, triple v);
+void gen_runarray56(stack *Stack)
+{
+  triple v=vm::pop<triple>(Stack);
+  realarray2 * t=vm::pop<realarray2 *>(Stack);
+#line 1384 "runarray.in"
+  {Stack->push<triple>(*t*v); return;}
+}
+
+#line 1388 "runarray.in"
+// pair project(triple v, realarray2 *t);
+void gen_runarray57(stack *Stack)
+{
+  realarray2 * t=vm::pop<realarray2 *>(Stack);
+  triple v=vm::pop<triple>(Stack);
+#line 1389 "runarray.in"
+  size_t n=checkArray(t);
+  if(n != 4) error(incommensurate);
+  array *t0=read<array*>(t,0);
+  array *t1=read<array*>(t,1);
+  array *t3=read<array*>(t,3);
+  if(checkArray(t0) != 4 || checkArray(t1) != 4 || checkArray(t3) != 4)
+    error(incommensurate);
+  
+  real x=v.getx();
+  real y=v.gety();
+  real z=v.getz();
+  
+  real f=read<real>(t3,0)*x+read<real>(t3,1)*y+read<real>(t3,2)*z+
+    read<real>(t3,3);
+  if(f == 0.0) dividebyzero();
+  f=1.0/f;
+  
+  {Stack->push<pair>(pair((read<real>(t0,0)*x+read<real>(t0,1)*y+read<real>(t0,2)*z+
+               read<real>(t0,3))*f,
+              (read<real>(t1,0)*x+read<real>(t1,1)*y+read<real>(t1,2)*z+
+               read<real>(t1,3))*f)); return;}
+}
+
+// Compute the dot product of vectors a and b.
+#line 1414 "runarray.in"
+// real dot(realarray *a, realarray *b);
+void gen_runarray58(stack *Stack)
+{
+  realarray * b=vm::pop<realarray *>(Stack);
+  realarray * a=vm::pop<realarray *>(Stack);
+#line 1415 "runarray.in"
+  size_t n=checkArrays(a,b);
+  real sum=0.0;
+  for(size_t i=0; i < n; ++i)
+    sum += read<real>(a,i)*read<real>(b,i);
+  {Stack->push<real>(sum); return;}
+}
+
+// Solve the problem L\inv f, where f is an n vector and L is the n x n matrix
+//
+// [ b[0] c[0]           a[0]   ]
+// [ a[1] b[1] c[1]             ]
+// [      a[2] b[2] c[2]        ]
+// [                ...         ]
+// [ c[n-1]       a[n-1] b[n-1] ]
+#line 1430 "runarray.in"
+// realarray* tridiagonal(realarray *a, realarray *b, realarray *c, realarray *f);
+void gen_runarray59(stack *Stack)
+{
+  realarray * f=vm::pop<realarray *>(Stack);
+  realarray * c=vm::pop<realarray *>(Stack);
+  realarray * b=vm::pop<realarray *>(Stack);
+  realarray * a=vm::pop<realarray *>(Stack);
+#line 1431 "runarray.in"
+  size_t n=checkArrays(a,b);
+  checkEqual(n,checkArray(c));
+  checkEqual(n,checkArray(f));
+  
+  array *up=new array(n);
+  array& u=*up;
+
+  if(n == 0) {Stack->push<realarray*>(up); return;}
+  
+  // Special case: zero Dirichlet boundary conditions
+  if(read<real>(a,0) == 0.0 && read<real>(c,n-1) == 0.0) {
+    real temp=read<real>(b,0);
+    if(temp == 0.0) dividebyzero();
+    temp=1.0/temp;
+    
+    real *work=new real[n];
+    u[0]=read<real>(f,0)*temp;
+    work[0]=-read<real>(c,0)*temp;
+        
+    for(size_t i=1; i < n; i++) {
+      real temp=(read<real>(b,i)+read<real>(a,i)*work[i-1]);
+      if(temp == 0.0) {delete[] work; dividebyzero();}
+      temp=1.0/temp;
+      u[i]=(read<real>(f,i)-read<real>(a,i)*read<real>(u,i-1))*temp;
+      work[i]=-read<real>(c,i)*temp;
+    }
+
+    for(size_t i=n-1; i >= 1; i--)
+      u[i-1]=read<real>(u,i-1)+work[i-1]*read<real>(u,i);
+    
+    delete[] work;
+    {Stack->push<realarray*>(up); return;}
+  }
+  
+  real binv=read<real>(b,0);
+  if(binv == 0.0) dividebyzero();
+  binv=1.0/binv;
+  
+  if(n == 1) {u[0]=read<real>(f,0)*binv; {Stack->push<realarray*>(up); return;}}
+  if(n == 2) {
+    real factor=(read<real>(b,0)*read<real>(b,1)-
+                 read<real>(a,0)*read<real>(c,1));
+    if(factor== 0.0) dividebyzero();
+    factor=1.0/factor;
+    real temp=(read<real>(b,0)*read<real>(f,1)-
+               read<real>(c,1)*read<real>(f,0))*factor;
+    u[0]=(read<real>(b,1)*read<real>(f,0)-
+          read<real>(a,0)*read<real>(f,1))*factor;
+    u[1]=temp;
+    {Stack->push<realarray*>(up); return;}
+  }
+        
+  real *gamma=new real[n-2];
+  real *delta=new real[n-2];
+  
+  gamma[0]=read<real>(c,0)*binv;
+  delta[0]=read<real>(a,0)*binv;
+  u[0]=read<real>(f,0)*binv;
+  real beta=read<real>(c,n-1);
+  real fn=read<real>(f,n-1)-beta*read<real>(u,0);
+  real alpha=read<real>(b,n-1)-beta*delta[0];
+
+  for(size_t i=1; i <= n-3; i++) {
+    real alphainv=read<real>(b,i)-read<real>(a,i)*gamma[i-1];
+    if(alphainv == 0.0) {delete[] gamma; delete[] delta; dividebyzero();}
+    alphainv=1.0/alphainv;
+    beta *= -gamma[i-1];
+    gamma[i]=read<real>(c,i)*alphainv;
+    u[i]=(read<real>(f,i)-read<real>(a,i)*read<real>(u,i-1))*alphainv;
+    fn -= beta*read<real>(u,i);
+    delta[i]=-read<real>(a,i)*delta[i-1]*alphainv;
+    alpha -= beta*delta[i];
+  }
+        
+  real alphainv=read<real>(b,n-2)-read<real>(a,n-2)*gamma[n-3];
+  if(alphainv == 0.0) {delete[] gamma; delete[] delta; dividebyzero();}
+  alphainv=1.0/alphainv;
+  u[n-2]=(read<real>(f,n-2)-read<real>(a,n-2)*read<real>(u,n-3))
+    *alphainv;
+  beta=read<real>(a,n-1)-beta*gamma[n-3];
+  real dnm1=(read<real>(c,n-2)-read<real>(a,n-2)*delta[n-3])*alphainv;
+  real temp=alpha-beta*dnm1;
+  if(temp == 0.0) {delete[] gamma; delete[] delta; dividebyzero();}
+  u[n-1]=temp=(fn-beta*read<real>(u,n-2))/temp;
+  u[n-2]=read<real>(u,n-2)-dnm1*temp;
+        
+  for(size_t i=n-2; i >= 1; i--)
+    u[i-1]=read<real>(u,i-1)-gamma[i-1]*read<real>(u,i)-delta[i-1]*temp;
+  
+  delete[] delta;
+  delete[] gamma;
+  
+  {Stack->push<realarray*>(up); return;}
+}
+
+// Root solve by Newton-Raphson
+#line 1528 "runarray.in"
+// real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x,            bool verbose=false);
+void gen_runarray60(stack *Stack)
+{
+  bool verbose=vm::pop<bool>(Stack,false);
+  real x=vm::pop<real>(Stack);
+  callableReal * fprime=vm::pop<callableReal *>(Stack);
+  callableReal * f=vm::pop<callableReal *>(Stack);
+  Int iterations=vm::pop<Int>(Stack,100);
+#line 1530 "runarray.in"
+  static const real fuzz=1000.0*DBL_EPSILON;
+  Int i=0;
+  size_t oldPrec=0;
+  if(verbose) 
+    oldPrec=cout.precision(DBL_DIG);
+
+  real diff=DBL_MAX;
+  real lastdiff;
+  do {
+    real x0=x;
+    
+    Stack->push(x);
+    fprime->call(Stack);
+    real dfdx=pop<real>(Stack);
+    
+    if(dfdx == 0.0) {
+      x=DBL_MAX;
+      break;
+    }
+
+    Stack->push(x);
+    f->call(Stack);
+    real fx=pop<real>(Stack);
+    
+    x -= fx/dfdx;
+
+    lastdiff=diff;
+    
+    if(verbose)
+      cout << "Newton-Raphson: " << x << endl;
+    
+    diff=fabs(x-x0);
+    if(++i == iterations) {
+      x=DBL_MAX;
+      break;
+    }
+  } while (diff != 0.0 && (diff < lastdiff || diff > fuzz*fabs(x)));
+
+  if(verbose)
+    cout.precision(oldPrec);
+  {Stack->push<real>(x); return;}
+}
+
+// Root solve by Newton-Raphson bisection
+// cf. routine rtsafe (Press et al.,  Numerical Recipes, 1991).
+#line 1576 "runarray.in"
+// real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x1,            real x2, bool verbose=false);
+void gen_runarray61(stack *Stack)
+{
+  bool verbose=vm::pop<bool>(Stack,false);
+  real x2=vm::pop<real>(Stack);
+  real x1=vm::pop<real>(Stack);
+  callableReal * fprime=vm::pop<callableReal *>(Stack);
+  callableReal * f=vm::pop<callableReal *>(Stack);
+  Int iterations=vm::pop<Int>(Stack,100);
+#line 1578 "runarray.in"
+  static const real fuzz=1000.0*DBL_EPSILON;
+  size_t oldPrec=0;
+  if(verbose) 
+    oldPrec=cout.precision(DBL_DIG);
+
+  Stack->push(x1);
+  f->call(Stack);
+  real f1=pop<real>(Stack);
+  if(f1 == 0.0) {Stack->push<real>(x1); return;}
+  
+  Stack->push(x2);
+  f->call(Stack);
+  real f2=pop<real>(Stack);
+  if(f2 == 0.0) {Stack->push<real>(x2); return;}
+        
+  if((f1 > 0.0 && f2 > 0.0) || (f1 < 0.0 && f2 < 0.0)) {
+    ostringstream buf;
+    buf << "root not bracketed, f(x1)=" << f1 << ", f(x2)=" << f2 << endl;
+    error(buf);
+  }
+
+  real x=0.5*(x1+x2);
+  real dxold=fabs(x2-x1);
+  if(f1 > 0.0) {
+    real temp=x1;
+    x1=x2;
+    x2=temp;
+  }
+        
+  if(verbose)
+    cout << "midpoint: " << x << endl;
+
+  real dx=dxold;
+  Stack->push(x);
+  f->call(Stack);
+  real y=pop<real>(Stack);
+  
+  Stack->push(x);
+  fprime->call(Stack);
+  real dy=pop<real>(Stack);
+
+  Int j;
+  for(j=0; j < iterations; j++) {
+    if(((x-x2)*dy-y)*((x-x1)*dy-y) >= 0.0 || fabs(2.0*y) > fabs(dxold*dy)) {
+      dxold=dx;
+      dx=0.5*(x2-x1);
+      x=x1+dx;
+      if(verbose)
+        cout << "bisection: " << x << endl;
+      if(x1 == x) {Stack->push<real>(x); return;}
+    } else {
+      dxold=dx;
+      dx=y/dy;
+      real temp=x;
+      x -= dx;
+      if(verbose)
+        cout << "Newton-Raphson: " << x << endl;
+      if(temp == x) {Stack->push<real>(x); return;}
+    }
+    if(fabs(dx) < fuzz*fabs(x)) {Stack->push<real>(x); return;}
+    
+    Stack->push(x);
+    f->call(Stack);
+    y=pop<real>(Stack);
+    
+    Stack->push(x);
+    fprime->call(Stack);
+    dy=pop<real>(Stack);
+
+    if(y < 0.0) x1=x;
+    else x2=x;
+  }
+  if(verbose)
+    cout.precision(oldPrec);
+  {Stack->push<real>((j == iterations) ? DBL_MAX : x); return;}
+}
+
+#line 1656 "runarray.in"
+// real simpson(callableReal *f, real a, real b, real acc=DBL_EPSILON,             real dxmax=0);
+void gen_runarray62(stack *Stack)
+{
+  real dxmax=vm::pop<real>(Stack,0);
+  real acc=vm::pop<real>(Stack,DBL_EPSILON);
+  real b=vm::pop<real>(Stack);
+  real a=vm::pop<real>(Stack);
+  callableReal * f=vm::pop<callableReal *>(Stack);
+#line 1658 "runarray.in"
+  real integral;
+  if(dxmax == 0) dxmax=b-a;
+  Func=f;
+  FuncStack=Stack;
+  if(!simpson(integral,wrapFunction,a,b,acc,dxmax))
+    error("nesting capacity exceeded in simpson");
+  {Stack->push<real>(integral); return;}
+}
+
+// Compute the fast Fourier transform of a pair array
+#line 1669 "runarray.in"
+void pairArrayFFT(stack *Stack)
+{
+  Int sign=vm::pop<Int>(Stack,1);
+  pairarray * a=vm::pop<pairarray *>(Stack);
+#line 1670 "runarray.in"
+  unsigned n=(unsigned) checkArray(a);
+#ifdef HAVE_LIBFFTW3
+  array *c=new array(n);
+  if(n) {
+    Complex *f=FFTWComplex(n);
+    fft1d Forward(n,intcast(sign),f);
+  
+    for(size_t i=0; i < n; i++) {
+      pair z=read<pair>(a,i);
+      f[i]=Complex(z.getx(),z.gety());
+    }
+    Forward.fft(f);
+  
+    for(size_t i=0; i < n; i++) {
+      Complex z=f[i];
+      (*c)[i]=pair(z.real(),z.imag());
+    }
+    FFTWdelete(f);
+  }
+#else
+  unused(&n);
+  unused(&sign);
+  array *c=new array(0);
+#endif //  HAVE_LIBFFTW3
+  {Stack->push<pairarray*>(c); return;}
+}
+
+#line 1698 "runarray.in"
+// Intarray2* triangulate(pairarray *z);
+void gen_runarray64(stack *Stack)
+{
+  pairarray * z=vm::pop<pairarray *>(Stack);
+#line 1699 "runarray.in"
+  size_t nv=checkArray(z);
+// Call robust version of Gilles Dumoulin's port of Paul Bourke's
+// triangulation code.
+
+  XYZ *pxyz=new XYZ[nv+3];
+  ITRIANGLE *V=new ITRIANGLE[4*nv];
+  
+  for(size_t i=0; i < nv; ++i) {
+    pair w=read<pair>(z,i);
+    pxyz[i].p[0]=w.getx();
+    pxyz[i].p[1]=w.gety();
+    pxyz[i].i=(Int) i;
+  }
+  
+  Int ntri;
+  Triangulate((Int) nv,pxyz,V,ntri,true,false);
+
+  size_t nt=(size_t) ntri;
+  array *t=new array(nt);
+  for(size_t i=0; i < nt; ++i) {
+    array *ti=new array(3);
+    (*t)[i]=ti;
+    ITRIANGLE *Vi=V+i;
+    (*ti)[0]=pxyz[Vi->p1].i;
+    (*ti)[1]=pxyz[Vi->p2].i;
+    (*ti)[2]=pxyz[Vi->p3].i;
+  }
+   
+  delete[] V;
+  delete[] pxyz;
+  {Stack->push<Intarray2*>(t); return;}
+}
+
+#line 1733 "runarray.in"
+// real norm(realarray *a);
+void gen_runarray65(stack *Stack)
+{
+  realarray * a=vm::pop<realarray *>(Stack);
+#line 1734 "runarray.in"
+  size_t n=checkArray(a);
+  real M=0.0;
+  for(size_t i=0; i < n; ++i) {
+    real x=fabs(vm::read<real>(a,i));
+    if(x > M) M=x;
+  }
+  {Stack->push<real>(M); return;}
+}
+
+#line 1744 "runarray.in"
+// real norm(realarray2 *a);
+void gen_runarray66(stack *Stack)
+{
+  realarray2 * a=vm::pop<realarray2 *>(Stack);
+#line 1745 "runarray.in"
+  size_t n=checkArray(a);
+  real M=0.0;
+  for(size_t i=0; i < n; ++i) {
+    vm::array *ai=vm::read<vm::array*>(a,i);
+    size_t m=checkArray(ai);
+    for(size_t j=0; j < m; ++j) {
+      real a=fabs(vm::read<real>(ai,j));
+      if(a > M) M=a;
+    }
+  }
+  {Stack->push<real>(M); return;}
+}
+
+#line 1759 "runarray.in"
+// real norm(triplearray2 *a);
+void gen_runarray67(stack *Stack)
+{
+  triplearray2 * a=vm::pop<triplearray2 *>(Stack);
+#line 1760 "runarray.in"
+  size_t n=checkArray(a);
+  real M=0.0;
+  for(size_t i=0; i < n; ++i) {
+    vm::array *ai=vm::read<vm::array*>(a,i);
+    size_t m=checkArray(ai);
+    for(size_t j=0; j < m; ++j) {
+      real a=vm::read<triple>(ai,j).abs2();
+      if(a > M) M=a;
+    }
+  }
+  {Stack->push<real>(sqrt(M)); return;}
+}
+
+#line 1774 "runarray.in"
+// real change2(triplearray2 *a);
+void gen_runarray68(stack *Stack)
+{
+  triplearray2 * a=vm::pop<triplearray2 *>(Stack);
+#line 1775 "runarray.in"
+  size_t n=checkArray(a);
+  if(n == 0) {Stack->push<real>(0.0); return;}
+  
+  vm::array *a0=vm::read<vm::array*>(a,0);
+  size_t m=checkArray(a0);
+  if(m == 0) {Stack->push<real>(0.0); return;}
+  triple a00=vm::read<triple>(a0,0);
+  real M=0.0;
+    
+  for(size_t i=0; i < n; ++i) {
+    vm::array *ai=vm::read<vm::array*>(a,i);
+    size_t m=checkArray(ai);
+    for(size_t j=0; j < m; ++j) {
+      real a=(vm::read<triple>(ai,j)-a00).abs2();
+      if(a > M) M=a;
+    }
+  }
+  {Stack->push<real>(M); return;}
+}
+
+#line 1796 "runarray.in"
+// triple minbezier(triplearray2 *P, triple b);
+void gen_runarray69(stack *Stack)
+{
+  triple b=vm::pop<triple>(Stack);
+  triplearray2 * P=vm::pop<triplearray2 *>(Stack);
+#line 1797 "runarray.in"
+  real *A=copyTripleArray2Components(P,true,4);
+  b=triple(bound(A,::min,b.getx(),sqrtFuzz*norm(A,16)),
+           bound(A+16,::min,b.gety(),sqrtFuzz*norm(A+16,16)),
+           bound(A+32,::min,b.getz(),sqrtFuzz*norm(A+32,16)));
+  delete[] A;
+  {Stack->push<triple>(b); return;}
+}
+
+#line 1806 "runarray.in"
+// triple maxbezier(triplearray2 *P, triple b);
+void gen_runarray70(stack *Stack)
+{
+  triple b=vm::pop<triple>(Stack);
+  triplearray2 * P=vm::pop<triplearray2 *>(Stack);
+#line 1807 "runarray.in"
+  real *A=copyTripleArray2Components(P,true,4);
+  b=triple(bound(A,::max,b.getx(),sqrtFuzz*norm(A,16)),
+           bound(A+16,::max,b.gety(),sqrtFuzz*norm(A+16,16)),
+           bound(A+32,::max,b.getz(),sqrtFuzz*norm(A+32,16)));
+  delete[] A;
+  {Stack->push<triple>(b); return;}
+}
+
+#line 1816 "runarray.in"
+// pair minratio(triplearray2 *P, pair b);
+void gen_runarray71(stack *Stack)
+{
+  pair b=vm::pop<pair>(Stack);
+  triplearray2 * P=vm::pop<triplearray2 *>(Stack);
+#line 1817 "runarray.in"
+  triple *A=copyTripleArray2C(P,true,4);
+  real fuzz=sqrtFuzz*norm(A,16);
+  b=pair(bound(A,::min,xratio,b.getx(),fuzz),
+         bound(A,::min,yratio,b.gety(),fuzz));
+  delete[] A;
+  {Stack->push<pair>(b); return;}
+}
+
+#line 1826 "runarray.in"
+// pair maxratio(triplearray2 *P, pair b);
+void gen_runarray72(stack *Stack)
+{
+  pair b=vm::pop<pair>(Stack);
+  triplearray2 * P=vm::pop<triplearray2 *>(Stack);
+#line 1827 "runarray.in"
+  triple *A=copyTripleArray2C(P,true,4);
+  real fuzz=sqrtFuzz*norm(A,16);
+  b=pair(bound(A,::max,xratio,b.getx(),fuzz),
+         bound(A,::max,yratio,b.gety(),fuzz));
+  delete[] A;
+  {Stack->push<pair>(b); return;}
+}
+
+} // namespace run
+
+namespace trans {
+
+void gen_runarray_venv(venv &ve)
+{
+#line 473 "runarray.in"
+  REGISTER_BLTIN(run::emptyArray,"emptyArray");
+#line 479 "runarray.in"
+  REGISTER_BLTIN(run::newDeepArray,"newDeepArray");
+#line 501 "runarray.in"
+  REGISTER_BLTIN(run::newInitializedArray,"newInitializedArray");
+#line 516 "runarray.in"
+  REGISTER_BLTIN(run::newAppendedArray,"newAppendedArray");
+#line 532 "runarray.in"
+  REGISTER_BLTIN(run::newDuplicateArray,"newDuplicateArray");
+#line 542 "runarray.in"
+  REGISTER_BLTIN(run::arrayRead,"arrayRead");
+#line 554 "runarray.in"
+  REGISTER_BLTIN(run::arraySliceRead,"arraySliceRead");
+#line 561 "runarray.in"
+  REGISTER_BLTIN(run::arraySliceReadToEnd,"arraySliceReadToEnd");
+#line 569 "runarray.in"
+  REGISTER_BLTIN(run::arrayArrayRead,"arrayArrayRead");
+#line 578 "runarray.in"
+  REGISTER_BLTIN(run::arrayWrite,"arrayWrite");
+#line 594 "runarray.in"
+  REGISTER_BLTIN(run::arraySliceWrite,"arraySliceWrite");
+#line 602 "runarray.in"
+  REGISTER_BLTIN(run::arraySliceWriteToEnd,"arraySliceWriteToEnd");
+#line 610 "runarray.in"
+  REGISTER_BLTIN(run::arrayLength,"arrayLength");
+#line 616 "runarray.in"
+  REGISTER_BLTIN(run::arrayKeys,"arrayKeys");
+#line 631 "runarray.in"
+  REGISTER_BLTIN(run::arrayCyclicFlag,"arrayCyclicFlag");
+#line 638 "runarray.in"
+  REGISTER_BLTIN(run::arraySetCyclicFlag,"arraySetCyclicFlag");
+#line 645 "runarray.in"
+  REGISTER_BLTIN(run::arrayInitializedHelper,"arrayInitializedHelper");
+#line 656 "runarray.in"
+  REGISTER_BLTIN(run::arrayInitialized,"arrayInitialized");
+#line 662 "runarray.in"
+  REGISTER_BLTIN(run::arrayCyclicHelper,"arrayCyclicHelper");
+#line 669 "runarray.in"
+  REGISTER_BLTIN(run::arrayCyclic,"arrayCyclic");
+#line 675 "runarray.in"
+  REGISTER_BLTIN(run::arrayPushHelper,"arrayPushHelper");
+#line 683 "runarray.in"
+  REGISTER_BLTIN(run::arrayPush,"arrayPush");
+#line 689 "runarray.in"
+  REGISTER_BLTIN(run::arrayAppendHelper,"arrayAppendHelper");
+#line 698 "runarray.in"
+  REGISTER_BLTIN(run::arrayAppend,"arrayAppend");
+#line 704 "runarray.in"
+  REGISTER_BLTIN(run::arrayPopHelper,"arrayPopHelper");
+#line 713 "runarray.in"
+  REGISTER_BLTIN(run::arrayPop,"arrayPop");
+#line 719 "runarray.in"
+  REGISTER_BLTIN(run::arrayInsertHelper,"arrayInsertHelper");
+#line 730 "runarray.in"
+  REGISTER_BLTIN(run::arrayInsert,"arrayInsert");
+#line 736 "runarray.in"
+  REGISTER_BLTIN(run::arrayDelete,"arrayDelete");
+#line 742 "runarray.in"
+  REGISTER_BLTIN(run::arrayAlias,"arrayAlias");
+#line 747 "runarray.in"
+  REGISTER_BLTIN(run::arrayIntArray,"arrayIntArray");
+#line 765 "runarray.in"
+  addFunc(ve, run::gen_runarray31, IntArray(), "complement", formal(IntArray(), "a", false, false), formal(primInt(), "n", false, false));
+#line 784 "runarray.in"
+  REGISTER_BLTIN(run::arraySequence,"arraySequence");
+#line 797 "runarray.in"
+  addFunc(ve, run::gen_runarray33, IntArray(), "sequence", formal(primInt(), "n", false, false));
+#line 808 "runarray.in"
+  REGISTER_BLTIN(run::arrayFunction,"arrayFunction");
+#line 821 "runarray.in"
+  REGISTER_BLTIN(run::arraySort,"arraySort");
+#line 830 "runarray.in"
+  addFunc(ve, run::gen_runarray36, primBoolean(), "all", formal(booleanArray(), "a", false, false));
+#line 839 "runarray.in"
+  addFunc(ve, run::gen_runarray37, booleanArray(), "!", formal(booleanArray(), "a", false, false));
+#line 848 "runarray.in"
+  addFunc(ve, run::gen_runarray38, primInt(), "sum", formal(booleanArray(), "a", false, false));
+#line 857 "runarray.in"
+  REGISTER_BLTIN(run::arrayCopy,"arrayCopy");
+#line 862 "runarray.in"
+  REGISTER_BLTIN(run::arrayConcat,"arrayConcat");
+#line 890 "runarray.in"
+  REGISTER_BLTIN(run::array2Copy,"array2Copy");
+#line 895 "runarray.in"
+  REGISTER_BLTIN(run::array3Copy,"array3Copy");
+#line 900 "runarray.in"
+  REGISTER_BLTIN(run::array2Transpose,"array2Transpose");
+#line 924 "runarray.in"
+  REGISTER_BLTIN(run::array3Transpose,"array3Transpose");
+#line 1011 "runarray.in"
+  addFunc(ve, run::gen_runarray45, primInt(), "find", formal(booleanArray(), "a", false, false), formal(primInt(), "n", true, false));
+#line 1030 "runarray.in"
+  REGISTER_BLTIN(run::arrayConditional,"arrayConditional");
+#line 1056 "runarray.in"
+  addFunc(ve, run::gen_runarray47, realArray2(), "identity", formal(primInt(), "n", false, false));
+#line 1062 "runarray.in"
+  REGISTER_BLTIN(run::diagonal,"diagonal");
+#line 1079 "runarray.in"
+  addFunc(ve, run::gen_runarray49, realArray2(), "inverse", formal(realArray2(), "a", false, false));
+#line 1171 "runarray.in"
+  addFunc(ve, run::gen_runarray50, realArray(), "solve", formal(realArray2(), "a", false, false), formal(realArray(), "b", false, false), formal(primBoolean(), "warn", true, false));
+#line 1222 "runarray.in"
+  addFunc(ve, run::gen_runarray51, realArray2(), "solve", formal(realArray2(), "a", false, false), formal(realArray2(), "b", false, false), formal(primBoolean(), "warn", true, false));
+#line 1286 "runarray.in"
+  addFunc(ve, run::gen_runarray52, primReal(), "determinant", formal(realArray2(), "a", false, false));
+#line 1302 "runarray.in"
+  addFunc(ve, run::gen_runarray53, realArray(), "*", formal(realArray2(), "a", false, false), formal(realArray(), "b", false, false));
+#line 1320 "runarray.in"
+  addFunc(ve, run::gen_runarray54, realArray(), "*", formal(realArray(), "a", false, false), formal(realArray2(), "b", false, false));
+#line 1348 "runarray.in"
+  addFunc(ve, run::gen_runarray55, realArray2(), "*", formal(realArray2(), "a", false, false), formal(realArray2(), "b", false, false));
+#line 1383 "runarray.in"
+  addFunc(ve, run::gen_runarray56, primTriple(), "*", formal(realArray2(), "t", false, false), formal(primTriple(), "v", false, false));
+#line 1388 "runarray.in"
+  addFunc(ve, run::gen_runarray57, primPair(), "project", formal(primTriple(), "v", false, false), formal(realArray2(), "t", false, false));
+#line 1413 "runarray.in"
+  addFunc(ve, run::gen_runarray58, primReal(), "dot", formal(realArray(), "a", false, false), formal(realArray(), "b", false, false));
+#line 1423 "runarray.in"
+  addFunc(ve, run::gen_runarray59, realArray(), "tridiagonal", formal(realArray(), "a", false, false), formal(realArray(), "b", false, false), formal(realArray(), "c", false, false), formal(realArray(), "f", false, false));
+#line 1527 "runarray.in"
+  addFunc(ve, run::gen_runarray60, primReal(), "newton", formal(primInt(), "iterations", true, false), formal(realRealFunction(), "f", false, false), formal(realRealFunction(), "fprime", false, false), formal(primReal(), "x", false, false), formal(primBoolean(), "verbose", true, false));
+#line 1574 "runarray.in"
+  addFunc(ve, run::gen_runarray61, primReal(), "newton", formal(primInt(), "iterations", true, false), formal(realRealFunction(), "f", false, false), formal(realRealFunction(), "fprime", false, false), formal(primReal(), "x1", false, false), formal(primReal(), "x2", false, false), formal(primBoolean(), "verbose", true, false));
+#line 1656 "runarray.in"
+  addFunc(ve, run::gen_runarray62, primReal(), "simpson", formal(realRealFunction(), "f", false, false), formal(primReal(), "a", false, false), formal(primReal(), "b", false, false), formal(primReal(), "acc", true, false), formal(primReal(), "dxmax", true, false));
+#line 1668 "runarray.in"
+  REGISTER_BLTIN(run::pairArrayFFT,"pairArrayFFT");
+#line 1698 "runarray.in"
+  addFunc(ve, run::gen_runarray64, IntArray2(), "triangulate", formal(pairArray(), "z", false, false));
+#line 1733 "runarray.in"
+  addFunc(ve, run::gen_runarray65, primReal(), "norm", formal(realArray(), "a", false, false));
+#line 1744 "runarray.in"
+  addFunc(ve, run::gen_runarray66, primReal(), "norm", formal(realArray2(), "a", false, false));
+#line 1759 "runarray.in"
+  addFunc(ve, run::gen_runarray67, primReal(), "norm", formal(tripleArray2(), "a", false, false));
+#line 1774 "runarray.in"
+  addFunc(ve, run::gen_runarray68, primReal(), "change2", formal(tripleArray2(), "a", false, false));
+#line 1796 "runarray.in"
+  addFunc(ve, run::gen_runarray69, primTriple(), "minbezier", formal(tripleArray2(), "p", false, false), formal(primTriple(), "b", false, false));
+#line 1806 "runarray.in"
+  addFunc(ve, run::gen_runarray70, primTriple(), "maxbezier", formal(tripleArray2(), "p", false, false), formal(primTriple(), "b", false, false));
+#line 1816 "runarray.in"
+  addFunc(ve, run::gen_runarray71, primPair(), "minratio", formal(tripleArray2(), "p", false, false), formal(primPair(), "b", false, false));
+#line 1826 "runarray.in"
+  addFunc(ve, run::gen_runarray72, primPair(), "maxratio", formal(tripleArray2(), "p", false, false), formal(primPair(), "b", false, false));
+}
+
+} // namespace trans