summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/runarray.cc
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/utils/asymptote/runarray.cc')
-rw-r--r--Build/source/utils/asymptote/runarray.cc2715
1 files changed, 2715 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..05c7f0a1a7c
--- /dev/null
+++ b/Build/source/utils/asymptote/runarray.cc
@@ -0,0 +1,2715 @@
+/***** 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"
+#include "opsymbols.h"
+
+using vm::stack;
+using vm::error;
+using vm::array;
+using vm::read;
+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 *copyTripleArray2Components(array *a, size_t &N,
+ GCPlacement placement=NoGC);
+triple *copyTripleArray2C(array *a, size_t &N,
+ GCPlacement placement=NoGC);
+}
+
+function *realRealFunction();
+
+#define CURRENTPEN processData().currentpen
+
+#line 22 "runarray.in"
+#include "array.h"
+#include "arrayop.h"
+#include "triple.h"
+#include "path3.h"
+#include "Delaunay.h"
+#include "glrender.h"
+
+#ifdef HAVE_LIBFFTW3
+#include "fftw++.h"
+ static const char *rectangular="matrix must be rectangular";
+#else
+static const char *installFFTW=
+ "Please install fftw3, run ./configure, and recompile";
+#endif
+
+using namespace camp;
+using namespace vm;
+
+namespace run {
+extern pair zero;
+}
+
+typedef array boolarray;
+typedef array Intarray;
+typedef array Intarray2;
+typedef array realarray;
+typedef array realarray2;
+typedef array pairarray;
+typedef array pairarray2;
+typedef array pairarray3;
+typedef array triplearray2;
+
+using types::booleanArray;
+using types::IntArray;
+using types::IntArray2;
+using types::realArray;
+using types::realArray2;
+using types::pairArray;
+using types::pairArray2;
+using types::pairArray3;
+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 const char *invalidarraylength="Invalid array length: ";
+
+static size_t *pivot,*Row,*Col;
+
+bound_double *bounddouble(int N)
+{
+ if(N == 16) return bound;
+ if(N == 10) return boundtri;
+ ostringstream buf;
+ buf << invalidarraylength << " " << N;
+ error(buf);
+ return NULL;
+}
+
+bound_triple *boundtriple(int N)
+{
+ if(N == 16) return bound;
+ if(N == 10) return boundtri;
+ ostringstream buf;
+ buf << invalidarraylength << " " << N;
+ error(buf);
+ return NULL;
+}
+
+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;
+}
+
+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;
+}
+
+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;
+}
+
+double *copyTripleArray2Components(array *a, size_t &N, GCPlacement placement)
+{
+ size_t n=checkArray(a);
+ N=0;
+ for(size_t i=0; i < n; i++)
+ N += checkArray(read<array*>(a,i));
+
+ double *A=(placement == NoGC) ? new double [3*N] :
+ new(placement) double[3*N];
+ double *p=A;
+
+ for(size_t i=0; i < n; i++) {
+ array *ai=read<array*>(a,i);
+ size_t m=checkArray(ai);
+ for(size_t j=0; j < m; j++) {
+ triple v=read<triple>(ai,j);
+ *p=v.getx();
+ *(p+N)=v.gety();
+ *(p+2*N)=v.getz();
+ ++p;
+ }
+ }
+ return A;
+}
+
+triple *copyTripleArray2C(array *a, size_t &N, GCPlacement placement)
+{
+ size_t n=checkArray(a);
+ N=0;
+ for(size_t i=0; i < n; i++)
+ N += checkArray(read<array*>(a,i));
+
+ triple *A=(placement == NoGC) ? new triple [N] :
+ new(placement) triple[N];
+ triple *p=A;
+
+ for(size_t i=0; i < n; i++) {
+ array *ai=read<array*>(a,i);
+ size_t m=checkArray(ai);
+ for(size_t j=0; j < m; j++)
+ *(p++)=read<triple>(ai,j);
+ }
+ return A;
+}
+
+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);
+}
+
+template<class T>
+array *mult(array *a, array *b)
+{
+ 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);
+
+ T *A,*B;
+ copyArray2C(A,a,false);
+ copyArray2C(B,b,false);
+
+ for(size_t i=0; i < n; ++i) {
+ T *Ai=A+i*nb;
+ array *ci=new array(nb0);
+ (*c)[i]=ci;
+ for(size_t j=0; j < nb0; ++j) {
+ T sum=T();
+ 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;
+
+ return c;
+}
+
+// Compute transpose(A)*A where A is an n x m matrix.
+template<class T>
+array *AtA(array *a)
+{
+ size_t n=checkArray(a);
+ size_t m=n == 0 ? 0 : checkArray(read<array*>(a,0));
+
+ array *c=new array(m);
+
+ T *A;
+ copyArray2C(A,a,false);
+
+ for(size_t i=0; i < m; ++i) {
+ array *ci=new array(m);
+ (*c)[i]=ci;
+ for(size_t j=0; j < m; ++j) {
+ T sum=T();
+ size_t kj=j;
+ size_t ki=i;
+ for(size_t k=0; k < n; ++k, kj += m, ki += m)
+ sum += A[ki]*A[kj];
+ (*ci)[j]=sum;
+ }
+ }
+
+ delete[] A;
+ return c;
+}
+
+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);
+}
+
+// Transpose an n x n matrix in place.
+void transpose(double *a, size_t n)
+{
+ for(size_t i=1; i < n; i++) {
+ for(size_t j=0; j < i; j++) {
+ size_t ij=n*i+j;
+ size_t ji=n*j+i;
+ double temp=a[ij];
+ a[ij]=a[ji];
+ a[ji]=temp;
+ }
+ }
+}
+
+// Invert an n x n array in place.
+void inverse(double *M, size_t n)
+{
+ if(n == 2) {
+ real a=M[0];
+ real b=M[1];
+ real c=M[2];
+ real d=M[3];
+ real det=a*d-b*c;
+ if(det == 0.0)
+ error(singular);
+ det=1.0/det;
+ M[0]=d*det;
+ M[1]=-b*det;
+ M[2]=-c*det;
+ M[3]=a*det;
+ return;
+ }
+
+ if(n == 3) {
+ real a=M[0], b=M[1], c=M[2];
+ real d=M[3], e=M[4], f=M[5];
+ real g=M[6], h=M[7], i=M[8];
+
+ real A=e*i-f*h;
+ real B=f*g-d*i;
+ real C=d*h-e*g;
+
+ real det=a*A+b*B+c*C;
+ if(det == 0.0)
+ error(singular);
+ det=1.0/det;
+
+ M[0]=A*det; M[1]=(c*h-b*i)*det; M[2]=(b*f-c*e)*det;
+ M[3]=B*det; M[4]=(a*i-c*g)*det; M[5]=(c*d-a*f)*det;
+ M[6]=C*det; M[7]=(b*g-a*h)*det; M[8]=(a*e-b*d)*det;
+ return;
+ }
+
+ 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++) {
+ double *aj=M+n*j;
+ if(pivot[j] != 1) {
+ for(size_t k=0; k < n; k++) {
+ if(pivot[k] == 0) {
+ real temp=fabs(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.
+ double *acol=M+n*col;
+ if(row != col) {
+ double *arow=M+n*row;
+ for(size_t k=0; k < n; k++) {
+ real temp=arow[k];
+ arow[k]=acol[k];
+ acol[k]=temp;
+ }
+ }
+
+ Row[i]=row;
+ Col[i]=col;
+
+ // Divide the pivot row by the pivot element.
+ real denom=acol[col];
+ if(denom == 0.0) {
+ inverseDeallocate();
+ error(singular);
+ }
+ real pivinv=1.0/denom;
+ acol[col]=1.0;
+ for(size_t k=0; k < n; k++)
+ acol[k]=acol[k]*pivinv;
+
+ // Reduce all rows except for the pivoted one.
+ for(size_t k=0; k < n; k++) {
+ if(k != col) {
+ double *ak=M+n*k;
+ real akcol=ak[col];
+ ak[col]=0.0;
+ for(size_t j=0; j < n; j++)
+ ak[j] -= acol[j]*akcol;
+ }
+ }
+ }
+
+ // Unscramble the inverse matrix in view of the column interchanges.
+ for(size_t k=n; k > 0;) {
+ k--;
+ size_t r=Row[k];
+ size_t c=Col[k];
+ if(r != c) {
+ for(size_t j=0; j < n; j++) {
+ double *aj=M+n*j;
+ real temp=aj[r];
+ aj[r]=aj[c];
+ aj[c]=temp;
+ }
+ }
+ }
+ inverseDeallocate();
+}
+}
+
+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);
+}
+
+// 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:
+
+
+
+#ifndef NOSYM
+#include "runarray.symbols.h"
+
+#endif
+namespace run {
+// Create an empty array.
+#line 594 "runarray.in"
+void emptyArray(stack *Stack)
+{
+#line 595 "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 604 "runarray.in"
+void newDeepArray(stack *Stack)
+{
+ Int depth=vm::pop<Int>(Stack);
+#line 605 "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 624 "runarray.in"
+void newInitializedArray(stack *Stack)
+{
+ Int n=vm::pop<Int>(Stack);
+#line 625 "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 638 "runarray.in"
+void newAppendedArray(stack *Stack)
+{
+ Int n=vm::pop<Int>(Stack);
+ array* tail=vm::pop<array*>(Stack);
+#line 639 "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;}
+}
+
+// Produce an array of n deep copies of value.
+// typeDepth is the true depth of the array determined at compile-time when the
+// operations for the array type are added. This typeDepth argument is
+// automatically pushed on the stack and is not visible to the user.
+#line 656 "runarray.in"
+void copyArrayValue(stack *Stack)
+{
+ Int typeDepth=vm::pop<Int>(Stack);
+ Int depth=vm::pop<Int>(Stack,Int_MAX);
+ item value=vm::pop(Stack);
+ Int n=vm::pop<Int>(Stack);
+#line 657 "runarray.in"
+ if(n < 0) error("cannot create a negative length array");
+ if(depth < 0) error("cannot copy to a negative depth");
+ if(depth > typeDepth) depth=typeDepth;
+ {Stack->push<array*>(new array((size_t) n, value, depth)); return;}
+}
+
+// Deep copy of array.
+// typeDepth is the true depth of the array determined at compile-time when the
+// operations for the array type are added. This typeDepth argument is
+// automatically pushed on the stack and is not visible to the user.
+#line 668 "runarray.in"
+void copyArray(stack *Stack)
+{
+ Int typeDepth=vm::pop<Int>(Stack);
+ Int depth=vm::pop<Int>(Stack,Int_MAX);
+ array * a=vm::pop<array *>(Stack);
+#line 669 "runarray.in"
+ if(a == 0) vm::error(dereferenceNullArray);
+ if(depth < 0) error("cannot copy to a negative depth");
+ if(depth > typeDepth) depth=typeDepth;
+ {Stack->push<array*>(a->copyToDepth(depth)); return;}
+}
+
+// Read an element from an array. Checks for initialization & bounds.
+#line 677 "runarray.in"
+void arrayRead(stack *Stack)
+{
+ Int n=vm::pop<Int>(Stack);
+ array * a=vm::pop<array *>(Stack);
+#line 678 "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 689 "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 690 "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 697 "runarray.in"
+void arraySliceReadToEnd(stack *Stack)
+{
+ Int left=vm::pop<Int>(Stack);
+ array * a=vm::pop<array *>(Stack);
+#line 698 "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 705 "runarray.in"
+void arrayArrayRead(stack *Stack)
+{
+ Int n=vm::pop<Int>(Stack);
+ array * a=vm::pop<array *>(Stack);
+#line 706 "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.
+// TODO: Add arrayWriteAndPop
+#line 714 "runarray.in"
+void arrayWrite(stack *Stack)
+{
+ item value=vm::pop(Stack);
+ Int n=vm::pop<Int>(Stack);
+ array * a=vm::pop<array *>(Stack);
+#line 715 "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 729 "runarray.in"
+void arraySliceWrite(stack *Stack)
+{
+ array * src=vm::pop<array *>(Stack);
+ Int right=vm::pop<Int>(Stack);
+ Int left=vm::pop<Int>(Stack);
+ array * dest=vm::pop<array *>(Stack);
+#line 730 "runarray.in"
+ checkArray(src);
+ checkArray(dest);
+ dest->setSlice(left, right, src);
+ {Stack->push<array*>(src); return;}
+}
+
+#line 737 "runarray.in"
+void arraySliceWriteToEnd(stack *Stack)
+{
+ array * src=vm::pop<array *>(Stack);
+ Int left=vm::pop<Int>(Stack);
+ array * dest=vm::pop<array *>(Stack);
+#line 738 "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 746 "runarray.in"
+void arrayLength(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+#line 747 "runarray.in"
+ {Stack->push<Int>((Int) checkArray(a)); return;}
+}
+
+// Returns an array of integers representing the keys of the array.
+#line 752 "runarray.in"
+void arrayKeys(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+#line 753 "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 767 "runarray.in"
+void arrayCyclicFlag(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+#line 768 "runarray.in"
+ checkArray(a);
+ {Stack->push<bool>(a->cyclic()); return;}
+}
+
+#line 773 "runarray.in"
+void arraySetCyclicFlag(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+ bool b=vm::pop<bool>(Stack);
+#line 774 "runarray.in"
+ checkArray(a);
+ a->cyclic(b);
+ {Stack->push<bool>(b); return;}
+}
+
+// Check to see if an array element is initialized.
+#line 781 "runarray.in"
+void arrayInitializedHelper(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+ Int n=vm::pop<Int>(Stack);
+#line 782 "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 792 "runarray.in"
+void arrayInitialized(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+#line 793 "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 798 "runarray.in"
+void arrayCyclicHelper(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+ bool b=vm::pop<bool>(Stack);
+#line 799 "runarray.in"
+ checkArray(a);
+ a->cyclic(b);
+}
+
+// Set the cyclic flag for an array.
+#line 805 "runarray.in"
+void arrayCyclic(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+#line 806 "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 811 "runarray.in"
+void arrayPushHelper(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+ item x=vm::pop(Stack);
+#line 812 "runarray.in"
+ checkArray(a);
+ a->push(x);
+ {Stack->push(x); return;}
+}
+
+// Returns the push method for an array.
+#line 819 "runarray.in"
+void arrayPush(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+#line 820 "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 825 "runarray.in"
+void arrayAppendHelper(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+ array * b=vm::pop<array *>(Stack);
+#line 826 "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 834 "runarray.in"
+void arrayAppend(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+#line 835 "runarray.in"
+ {Stack->push<callable*>(new thunk(new bfunc(arrayAppendHelper),a)); return;}
+}
+
+// The helper function for the pop method.
+#line 840 "runarray.in"
+void arrayPopHelper(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+#line 841 "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 849 "runarray.in"
+void arrayPop(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+#line 850 "runarray.in"
+ {Stack->push<callable*>(new thunk(new bfunc(arrayPopHelper),a)); return;}
+}
+
+// The helper function for the insert method.
+#line 855 "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 856 "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 866 "runarray.in"
+void arrayInsert(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+#line 867 "runarray.in"
+ {Stack->push<callable*>(new thunk(new bfunc(arrayInsertHelper),a)); return;}
+}
+
+// Returns the delete method for an array.
+#line 872 "runarray.in"
+void arrayDelete(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+#line 873 "runarray.in"
+ {Stack->push<callable*>(new thunk(new bfunc(arrayDeleteHelper),a)); return;}
+}
+
+#line 877 "runarray.in"
+void arrayAlias(stack *Stack)
+{
+ array * b=vm::pop<array *>(Stack);
+ array * a=vm::pop<array *>(Stack);
+#line 878 "runarray.in"
+ {Stack->push<bool>(a==b); return;}
+}
+
+// Return array formed by indexing array a with elements of integer array b
+#line 883 "runarray.in"
+void arrayIntArray(stack *Stack)
+{
+ array * b=vm::pop<array *>(Stack);
+ array * a=vm::pop<array *>(Stack);
+#line 884 "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 902 "runarray.in"
+// Intarray* complement(Intarray *a, Int n);
+void gen_runarray32(stack *Stack)
+{
+ Int n=vm::pop<Int>(Stack);
+ Intarray * a=vm::pop<Intarray *>(Stack);
+#line 903 "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 920 "runarray.in"
+void arraySequence(stack *Stack)
+{
+ Int n=vm::pop<Int>(Stack);
+ callable * f=vm::pop<callable *>(Stack);
+#line 921 "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 933 "runarray.in"
+// Intarray* sequence(Int n);
+void gen_runarray34(stack *Stack)
+{
+ Int n=vm::pop<Int>(Stack);
+#line 934 "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 944 "runarray.in"
+void arrayFunction(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+ callable * f=vm::pop<callable *>(Stack);
+#line 945 "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 956 "runarray.in"
+void arraySort(stack *Stack)
+{
+ bool stable=vm::pop<bool>(Stack,true);
+ callable * less=vm::pop<callable *>(Stack);
+ array * a=vm::pop<array *>(Stack);
+#line 957 "runarray.in"
+ array *c=copyArray(a);
+ compareFunc=less;
+ FuncStack=Stack;
+ if(stable) stable_sort(c->begin(),c->end(),compareFunction);
+ else sort(c->begin(),c->end(),compareFunction);
+ {Stack->push<array*>(c); return;}
+}
+
+#line 966 "runarray.in"
+void arraySearch(stack *Stack)
+{
+ callable * less=vm::pop<callable *>(Stack);
+ item key=vm::pop(Stack);
+ array * a=vm::pop<array *>(Stack);
+#line 967 "runarray.in"
+ size_t size=a->size();
+ compareFunc=less;
+ FuncStack=Stack;
+ if(size == 0 || compareFunction(key,(*a)[0])) {Stack->push<Int>(-1); return;}
+ size_t u=size-1;
+ if(!compareFunction(key,(*a)[u])) {Stack->push<Int>(Intcast(u)); return;}
+ size_t l=0;
+
+ while (l < u) {
+ size_t i=(l+u)/2;
+ if(compareFunction(key,(*a)[i])) u=i;
+ else if(compareFunction(key,(*a)[i+1])) {Stack->push<Int>(Intcast(i)); return;}
+ else l=i+1;
+ }
+ {Stack->push<Int>(0); return;}
+}
+
+#line 985 "runarray.in"
+// bool all(boolarray *a);
+void gen_runarray38(stack *Stack)
+{
+ boolarray * a=vm::pop<boolarray *>(Stack);
+#line 986 "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 994 "runarray.in"
+// boolarray* !(boolarray* a);
+void gen_runarray39(stack *Stack)
+{
+ boolarray* a=vm::pop<boolarray*>(Stack);
+#line 995 "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 1003 "runarray.in"
+// Int sum(boolarray *a);
+void gen_runarray40(stack *Stack)
+{
+ boolarray * a=vm::pop<boolarray *>(Stack);
+#line 1004 "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 1012 "runarray.in"
+void arrayConcat(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+#line 1013 "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 1040 "runarray.in"
+void array2Transpose(stack *Stack)
+{
+ array * a=vm::pop<array *>(Stack);
+#line 1041 "runarray.in"
+ size_t asize=checkArray(a);
+ array *c=new array(0);
+ size_t csize=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);
+ if(c->size() < aisize) {
+ c->resize(aisize);
+ for(size_t j=csize; j < aisize; j++)
+ (*c)[j]=new array(0);
+ csize=aisize;
+ }
+ for(size_t j=0; j < aisize; j++) {
+ if(!(*ai)[j].empty()) {
+ 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 1070 "runarray.in"
+void array3Transpose(stack *Stack)
+{
+ array * perm=vm::pop<array *>(Stack);
+ array * a=vm::pop<array *>(Stack);
+#line 1071 "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;}
+}
+
+// Find the index of the nth true value in a boolean array or -1 if not found.
+// If n is negative, search backwards.
+#line 1155 "runarray.in"
+// Int find(boolarray *a, Int n=1);
+void gen_runarray44(stack *Stack)
+{
+ Int n=vm::pop<Int>(Stack,1);
+ boolarray * a=vm::pop<boolarray *>(Stack);
+#line 1156 "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;}
+}
+
+// Find all indices of true values in a boolean array.
+#line 1173 "runarray.in"
+// Intarray* findall(boolarray *a);
+void gen_runarray45(stack *Stack)
+{
+ boolarray * a=vm::pop<boolarray *>(Stack);
+#line 1174 "runarray.in"
+ size_t size=checkArray(a);
+ array *b=new array(0);
+ for(size_t i=0; i < size; i++) {
+ if(read<bool>(a,i)) {
+ b->push((Int) i);
+ }
+ }
+ {Stack->push<Intarray*>(b); 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 1187 "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 1188 "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 1212 "runarray.in"
+// realarray2* identity(Int n);
+void gen_runarray47(stack *Stack)
+{
+ Int n=vm::pop<Int>(Stack);
+#line 1213 "runarray.in"
+ {Stack->push<realarray2*>(Identity(n)); return;}
+}
+
+// Return the inverse of an n x n matrix a using Gauss-Jordan elimination.
+#line 1218 "runarray.in"
+// realarray2* inverse(realarray2 *a);
+void gen_runarray48(stack *Stack)
+{
+ realarray2 * a=vm::pop<realarray2 *>(Stack);
+#line 1219 "runarray.in"
+ size_t n=checkArray(a);
+ double *A;
+ copyArray2C(A,a,true,0,NoGC);
+ inverse(A,n);
+ a=copyCArray2(n,n,A);
+ delete[] A;
+ {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 1232 "runarray.in"
+// realarray* solve(realarray2 *a, realarray *b, bool warn=true);
+void gen_runarray49(stack *Stack)
+{
+ bool warn=vm::pop<bool>(Stack,true);
+ realarray * b=vm::pop<realarray *>(Stack);
+ realarray2 * a=vm::pop<realarray2 *>(Stack);
+#line 1233 "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,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,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 1285 "runarray.in"
+// realarray2* solve(realarray2 *a, realarray2 *b, bool warn=true);
+void gen_runarray50(stack *Stack)
+{
+ bool warn=vm::pop<bool>(Stack,true);
+ realarray2 * b=vm::pop<realarray2 *>(Stack);
+ realarray2 * a=vm::pop<realarray2 *>(Stack);
+#line 1286 "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,*B;
+ copyArray2C(A,a);
+ copyArray2C(B,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 1348 "runarray.in"
+// real determinant(realarray2 *a);
+void gen_runarray51(stack *Stack)
+{
+ realarray2 * a=vm::pop<realarray2 *>(Stack);
+#line 1349 "runarray.in"
+ real *A;
+ copyArray2C(A,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 1364 "runarray.in"
+// realarray* *(realarray2 *a, realarray *b);
+void gen_runarray52(stack *Stack)
+{
+ realarray * b=vm::pop<realarray *>(Stack);
+ realarray2 * a=vm::pop<realarray2 *>(Stack);
+#line 1365 "runarray.in"
+ size_t n=checkArray(a);
+ size_t m=checkArray(b);
+ array *c=new array(n);
+ real *B;
+ copyArrayC(B,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 1383 "runarray.in"
+// realarray* *(realarray *a, realarray2 *b);
+void gen_runarray53(stack *Stack)
+{
+ realarray2 * b=vm::pop<realarray2 *>(Stack);
+ realarray * a=vm::pop<realarray *>(Stack);
+#line 1384 "runarray.in"
+ size_t n=checkArray(a);
+ if(n != checkArray(b)) error(incommensurate);
+ real *A;
+ copyArrayC(A,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 1412 "runarray.in"
+// Intarray2* *(Intarray2 *a, Intarray2 *b);
+void gen_runarray54(stack *Stack)
+{
+ Intarray2 * b=vm::pop<Intarray2 *>(Stack);
+ Intarray2 * a=vm::pop<Intarray2 *>(Stack);
+#line 1413 "runarray.in"
+ {Stack->push<Intarray2*>(mult<Int>(a,b)); return;}
+}
+
+#line 1417 "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 1418 "runarray.in"
+ {Stack->push<realarray2*>(mult<real>(a,b)); return;}
+}
+
+#line 1422 "runarray.in"
+// pairarray2* *(pairarray2 *a, pairarray2 *b);
+void gen_runarray56(stack *Stack)
+{
+ pairarray2 * b=vm::pop<pairarray2 *>(Stack);
+ pairarray2 * a=vm::pop<pairarray2 *>(Stack);
+#line 1423 "runarray.in"
+ {Stack->push<pairarray2*>(mult<pair>(a,b)); return;}
+}
+
+#line 1427 "runarray.in"
+// triple *(realarray2 *t, triple v);
+void gen_runarray57(stack *Stack)
+{
+ triple v=vm::pop<triple>(Stack);
+ realarray2 * t=vm::pop<realarray2 *>(Stack);
+#line 1428 "runarray.in"
+ {Stack->push<triple>(*t*v); return;}
+}
+
+#line 1432 "runarray.in"
+// realarray2* AtA(realarray2 *a);
+void gen_runarray58(stack *Stack)
+{
+ realarray2 * a=vm::pop<realarray2 *>(Stack);
+#line 1433 "runarray.in"
+ {Stack->push<realarray2*>(AtA<real>(a)); return;}
+}
+
+#line 1437 "runarray.in"
+// pair project(triple v, realarray2 *t);
+void gen_runarray59(stack *Stack)
+{
+ realarray2 * t=vm::pop<realarray2 *>(Stack);
+ triple v=vm::pop<triple>(Stack);
+#line 1438 "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 1463 "runarray.in"
+// real dot(realarray *a, realarray *b);
+void gen_runarray60(stack *Stack)
+{
+ realarray * b=vm::pop<realarray *>(Stack);
+ realarray * a=vm::pop<realarray *>(Stack);
+#line 1464 "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;}
+}
+
+// Compute the complex dot product of vectors a and b.
+#line 1473 "runarray.in"
+// pair dot(pairarray *a, pairarray *b);
+void gen_runarray61(stack *Stack)
+{
+ pairarray * b=vm::pop<pairarray *>(Stack);
+ pairarray * a=vm::pop<pairarray *>(Stack);
+#line 1474 "runarray.in"
+ size_t n=checkArrays(a,b);
+ pair sum=zero;
+ for(size_t i=0; i < n; ++i)
+ sum += read<pair>(a,i)*conj(read<pair>(b,i));
+ {Stack->push<pair>(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 1489 "runarray.in"
+// realarray* tridiagonal(realarray *a, realarray *b, realarray *c, realarray *f);
+void gen_runarray62(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 1490 "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 1587 "runarray.in"
+// real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x, bool verbose=false);
+void gen_runarray63(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 1589 "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 1635 "runarray.in"
+// real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x1, real x2, bool verbose=false);
+void gen_runarray64(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 1637 "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;}
+}
+
+// Find a root for the specified continuous (but not necessarily
+// differentiable) function. Whatever value t is returned, it is guaranteed
+// that t is within [a, b] and within tolerance of a sign change.
+// An error is thrown if fa and fb are both positive or both negative.
+//
+// In this implementation, the binary search is interleaved
+// with a modified version of quadratic interpolation.
+// This is a C++ port of the Asymptote routine written by Charles Staats III.
+#line 1723 "runarray.in"
+// real _findroot(callableReal *f, real a, real b, real tolerance, real fa, real fb);
+void gen_runarray65(stack *Stack)
+{
+ real fb=vm::pop<real>(Stack);
+ real fa=vm::pop<real>(Stack);
+ real tolerance=vm::pop<real>(Stack);
+ real b=vm::pop<real>(Stack);
+ real a=vm::pop<real>(Stack);
+ callableReal * f=vm::pop<callableReal *>(Stack);
+#line 1725 "runarray.in"
+ if(fa == 0.0) {Stack->push<real>(a); return;}
+ if(fb == 0.0) {Stack->push<real>(b); return;}
+
+ const char* oppsign="fa and fb must have opposite signs";
+ int sign;
+
+ if(fa < 0.0) {
+ if(fb < 0.0) error(oppsign);
+ sign=1;
+ } else {
+ if(fb > 0.0) error(oppsign);
+ fa=-fa;
+ fb=-fb;
+ sign=-1;
+ }
+
+ real t=a;
+ real ft=fa;
+ real twicetolerance=2.0*tolerance;
+
+ while(b-a > tolerance) {
+ t=(a+b)*0.5;
+
+ Stack->push(t);
+ f->call(Stack);
+ ft=sign*pop<double>(Stack);
+ if(ft == 0.0) {Stack->push<real>(t); return;}
+
+ // If halving the interval already puts us within tolerance,
+ // don't bother with the interpolation step.
+ if(b-a >= twicetolerance) {
+
+ real factor=1.0/(b-a);
+ real q_A=2.0*(fa-2.0*ft+fb)*factor*factor;
+ real q_B=(fb-fa)*factor;
+ quadraticroots Q=quadraticroots(q_A,q_B,ft);
+
+ // If the interpolation somehow failed, continue on to the next binary
+ // search step. This may or may not be possible, depending on what
+ // theoretical guarantees are provided by the quadraticroots function.
+
+ real root;
+ bool found=Q.roots > 0;
+ if(found) {
+ root=t+Q.t1;
+ if(root <= a || root >= b) {
+ if(Q.roots == 1) found=false;
+ else {
+ root=t+Q.t2;
+ if(root <= a || root >= b) found=false;
+ }
+ }
+ }
+
+ if(found) {
+ if(ft > 0.0) {
+ b=t;
+ fb=ft;
+ } else {
+ a=t;
+ fa=ft;
+ }
+
+ t=root;
+
+ // If the interpolated value is close to one edge of
+ // the interval, move it farther away from the edge in
+ // an effort to catch the root in the middle.
+ real margin=(b-a)*1.0e-3;
+ if(t-a < margin) t=a+2.0*(t-a);
+ else if(b-t < margin) t=b-2.0*(b-t);
+
+ Stack->push(t);
+ f->call(Stack);
+ ft=sign*pop<double>(Stack);
+
+ if(ft == 0.0) {Stack->push<real>(t); return;}
+ }
+ }
+
+ if(ft > 0.0) {
+ b=t;
+ fb=ft;
+ } else if(ft < 0.0) {
+ a=t;
+ fa=ft;
+ }
+ }
+ {Stack->push<real>(a-(b-a)/(fb-fa)*fa); return;}
+}
+
+#line 1817 "runarray.in"
+// real simpson(callableReal *f, real a, real b, real acc=DBL_EPSILON, real dxmax=0);
+void gen_runarray66(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 1819 "runarray.in"
+ real integral;
+ if(dxmax <= 0) dxmax=fabs(b-a);
+ callable *oldFunc=Func;
+ Func=f;
+ FuncStack=Stack;
+ if(!simpson(integral,wrapFunction,a,b,acc,dxmax))
+ error("nesting capacity exceeded in simpson");
+ Func=oldFunc;
+ {Stack->push<real>(integral); return;}
+}
+
+// Compute the fast Fourier transform of a pair array
+#line 1832 "runarray.in"
+// pairarray* fft(pairarray *a, Int sign=1);
+void gen_runarray67(stack *Stack)
+{
+ Int sign=vm::pop<Int>(Stack,1);
+ pairarray * a=vm::pop<pairarray *>(Stack);
+#line 1833 "runarray.in"
+#ifdef HAVE_LIBFFTW3
+ unsigned n=(unsigned) checkArray(a);
+ array *c=new array(n);
+ if(n) {
+ Complex *f=utils::ComplexAlign(n);
+ fftwpp::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());
+ }
+ utils::deleteAlign(f);
+ }
+#else
+ unused(a);
+ unused(&sign);
+ array *c=new array(0);
+ error(installFFTW);
+#endif // HAVE_LIBFFTW3
+ {Stack->push<pairarray*>(c); return;}
+}
+
+// Compute the fast Fourier transform of a 2D pair array
+#line 1863 "runarray.in"
+// pairarray2* fft(pairarray2 *a, Int sign=1);
+void gen_runarray68(stack *Stack)
+{
+ Int sign=vm::pop<Int>(Stack,1);
+ pairarray2 * a=vm::pop<pairarray2 *>(Stack);
+#line 1864 "runarray.in"
+#ifdef HAVE_LIBFFTW3
+ size_t n=checkArray(a);
+ size_t m=n == 0 ? 0 : checkArray(read<array*>(a,0));
+
+ array *c=new array(n);
+ Complex *f=utils::ComplexAlign(n*m);
+ fftwpp::fft2d Forward(n,m,intcast(sign),f);
+
+ if(n) {
+ for(size_t i=0; i < n; ++i) {
+ array *ai=read<array *>(a,i);
+ size_t aisize=checkArray(ai);
+ if(aisize != m) error(rectangular);
+ Complex *fi=f+m*i;
+ for(size_t j=0; j < m; ++j) {
+ pair z=read<pair>(ai,j);
+ fi[j]=Complex(z.getx(),z.gety());
+ }
+ }
+
+ Forward.fft(f);
+
+ for(size_t i=0; i < n; ++i) {
+ array *ci=new array(m);
+ (*c)[i]=ci;
+ Complex *fi=f+m*i;
+ for(size_t j=0; j < m; ++j) {
+ Complex z=fi[j];
+ (*ci)[j]=pair(z.real(),z.imag());
+ }
+ }
+
+ utils::deleteAlign(f);
+ }
+#else
+ unused(a);
+ unused(&sign);
+ array *c=new array(0);
+ error(installFFTW);
+#endif // HAVE_LIBFFTW3
+ {Stack->push<pairarray2*>(c); return;}
+}
+
+// Compute the fast Fourier transform of a 3D pair array
+#line 1909 "runarray.in"
+// pairarray3* fft(pairarray3 *a, Int sign=1);
+void gen_runarray69(stack *Stack)
+{
+ Int sign=vm::pop<Int>(Stack,1);
+ pairarray3 * a=vm::pop<pairarray3 *>(Stack);
+#line 1910 "runarray.in"
+#ifdef HAVE_LIBFFTW3
+ size_t n=checkArray(a);
+ array *a0=read<array*>(a,0);
+ size_t m=n == 0 ? 0 : checkArray(a0);
+ size_t l=m == 0 ? 0 : checkArray(read<array*>(a0,0));
+
+ array *c=new array(n);
+ Complex *f=utils::ComplexAlign(n*m*l);
+ fftwpp::fft3d Forward(n,m,l,intcast(sign),f);
+
+ if(n) {
+ for(size_t i=0; i < n; ++i) {
+ array *ai=read<array *>(a,i);
+ size_t aisize=checkArray(ai);
+ if(aisize != m) error(rectangular);
+ Complex *fi=f+m*l*i;
+ for(size_t j=0; j < m; ++j) {
+ array *aij=read<array *>(ai,j);
+ size_t aijsize=checkArray(aij);
+ if(aijsize != l) error(rectangular);
+ Complex *fij=fi+l*j;
+ for(size_t k=0; k < l; ++k) {
+ pair z=read<pair>(aij,k);
+ fij[k]=Complex(z.getx(),z.gety());
+ }
+ }
+ }
+
+ Forward.fft(f);
+
+ for(size_t i=0; i < n; ++i) {
+ array *ci=new array(m);
+ (*c)[i]=ci;
+ Complex *fi=f+m*l*i;
+ for(size_t j=0; j < m; ++j) {
+ array *cij=new array(l);
+ (*ci)[j]=cij;
+ Complex *fij=fi+l*j;
+ for(size_t k=0; k < l; ++k) {
+ Complex z=fij[k];
+ (*cij)[k]=pair(z.real(),z.imag());
+ }
+ }
+ }
+
+ utils::deleteAlign(f);
+ }
+#else
+ unused(a);
+ unused(&sign);
+ array *c=new array(0);
+ error(installFFTW);
+#endif // HAVE_LIBFFTW3
+ {Stack->push<pairarray3*>(c); return;}
+}
+
+#line 1967 "runarray.in"
+// Intarray2* triangulate(pairarray *z);
+void gen_runarray70(stack *Stack)
+{
+ pairarray * z=vm::pop<pairarray *>(Stack);
+#line 1968 "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 2002 "runarray.in"
+// real norm(realarray *a);
+void gen_runarray71(stack *Stack)
+{
+ realarray * a=vm::pop<realarray *>(Stack);
+#line 2003 "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 2013 "runarray.in"
+// real norm(realarray2 *a);
+void gen_runarray72(stack *Stack)
+{
+ realarray2 * a=vm::pop<realarray2 *>(Stack);
+#line 2014 "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 2028 "runarray.in"
+// real norm(triplearray2 *a);
+void gen_runarray73(stack *Stack)
+{
+ triplearray2 * a=vm::pop<triplearray2 *>(Stack);
+#line 2029 "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 2043 "runarray.in"
+// real change2(triplearray2 *a);
+void gen_runarray74(stack *Stack)
+{
+ triplearray2 * a=vm::pop<triplearray2 *>(Stack);
+#line 2044 "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 2065 "runarray.in"
+// triple minbezier(triplearray2 *P, triple b);
+void gen_runarray75(stack *Stack)
+{
+ triple b=vm::pop<triple>(Stack);
+ triplearray2 * P=vm::pop<triplearray2 *>(Stack);
+#line 2066 "runarray.in"
+ size_t N;
+ real *A=copyTripleArray2Components(P,N);
+ bound_double *B=bounddouble(N);
+ b=triple(B(A,::min,b.getx(),Fuzz*norm(A,N),maxdepth),
+ B(A+N,::min,b.gety(),Fuzz*norm(A+N,N),maxdepth),
+ B(A+2*N,::min,b.getz(),Fuzz*norm(A+2*N,N),maxdepth));
+ delete[] A;
+ {Stack->push<triple>(b); return;}
+}
+
+#line 2077 "runarray.in"
+// triple maxbezier(triplearray2 *P, triple b);
+void gen_runarray76(stack *Stack)
+{
+ triple b=vm::pop<triple>(Stack);
+ triplearray2 * P=vm::pop<triplearray2 *>(Stack);
+#line 2078 "runarray.in"
+ size_t N;
+ real *A=copyTripleArray2Components(P,N);
+ bound_double *B=bounddouble(N);
+ b=triple(B(A,::max,b.getx(),Fuzz*norm(A,N),maxdepth),
+ B(A+N,::max,b.gety(),Fuzz*norm(A+N,N),maxdepth),
+ B(A+2*N,::max,b.getz(),Fuzz*norm(A+2*N,N),maxdepth));
+ delete[] A;
+ {Stack->push<triple>(b); return;}
+}
+
+#line 2089 "runarray.in"
+// pair minratio(triplearray2 *P, pair b);
+void gen_runarray77(stack *Stack)
+{
+ pair b=vm::pop<pair>(Stack);
+ triplearray2 * P=vm::pop<triplearray2 *>(Stack);
+#line 2090 "runarray.in"
+ size_t N;
+ triple *A=copyTripleArray2C(P,N);
+ real fuzz=Fuzz*norm(A,N);
+ bound_triple *B=boundtriple(N);
+ b=pair(B(A,::min,xratio,b.getx(),fuzz,maxdepth),
+ B(A,::min,yratio,b.gety(),fuzz,maxdepth));
+ delete[] A;
+ {Stack->push<pair>(b); return;}
+}
+
+#line 2101 "runarray.in"
+// pair maxratio(triplearray2 *P, pair b);
+void gen_runarray78(stack *Stack)
+{
+ pair b=vm::pop<pair>(Stack);
+ triplearray2 * P=vm::pop<triplearray2 *>(Stack);
+#line 2102 "runarray.in"
+ size_t N;
+ triple *A=copyTripleArray2C(P,N);
+ bound_triple *B=boundtriple(N);
+ real fuzz=Fuzz*norm(A,N);
+ b=pair(B(A,::max,xratio,b.getx(),fuzz,maxdepth),
+ B(A,::max,yratio,b.gety(),fuzz,maxdepth));
+ delete[] A;
+ {Stack->push<pair>(b); return;}
+}
+
+#line 2113 "runarray.in"
+// realarray* _projection();
+void gen_runarray79(stack *Stack)
+{
+#line 2114 "runarray.in"
+#ifdef HAVE_GL
+ array *a=new array(14);
+ gl::projection P=gl::camera();
+ size_t k=0;
+ (*a)[k++]=P.orthographic ? 1.0 : 0.0;
+
+ triple camera=P.camera;
+ (*a)[k++]=camera.getx();
+ (*a)[k++]=camera.gety();
+ (*a)[k++]=camera.getz();
+
+ triple up=P.up;
+ (*a)[k++]=up.getx();
+ (*a)[k++]=up.gety();
+ (*a)[k++]=up.getz();
+
+ triple target=P.target;
+ (*a)[k++]=target.getx();
+ (*a)[k++]=target.gety();
+ (*a)[k++]=target.getz();
+
+ (*a)[k++]=P.zoom;
+ (*a)[k++]=P.angle;
+
+ (*a)[k++]=P.viewportshift.getx();
+ (*a)[k++]=P.viewportshift.gety();
+#endif
+ {Stack->push<realarray*>(new array(0)); return;}
+}
+
+} // namespace run
+
+namespace trans {
+
+void gen_runarray_venv(venv &ve)
+{
+#line 593 "runarray.in"
+ REGISTER_BLTIN(run::emptyArray,"emptyArray");
+#line 599 "runarray.in"
+ REGISTER_BLTIN(run::newDeepArray,"newDeepArray");
+#line 621 "runarray.in"
+ REGISTER_BLTIN(run::newInitializedArray,"newInitializedArray");
+#line 636 "runarray.in"
+ REGISTER_BLTIN(run::newAppendedArray,"newAppendedArray");
+#line 652 "runarray.in"
+ REGISTER_BLTIN(run::copyArrayValue,"copyArrayValue");
+#line 664 "runarray.in"
+ REGISTER_BLTIN(run::copyArray,"copyArray");
+#line 676 "runarray.in"
+ REGISTER_BLTIN(run::arrayRead,"arrayRead");
+#line 688 "runarray.in"
+ REGISTER_BLTIN(run::arraySliceRead,"arraySliceRead");
+#line 695 "runarray.in"
+ REGISTER_BLTIN(run::arraySliceReadToEnd,"arraySliceReadToEnd");
+#line 703 "runarray.in"
+ REGISTER_BLTIN(run::arrayArrayRead,"arrayArrayRead");
+#line 712 "runarray.in"
+ REGISTER_BLTIN(run::arrayWrite,"arrayWrite");
+#line 729 "runarray.in"
+ REGISTER_BLTIN(run::arraySliceWrite,"arraySliceWrite");
+#line 737 "runarray.in"
+ REGISTER_BLTIN(run::arraySliceWriteToEnd,"arraySliceWriteToEnd");
+#line 745 "runarray.in"
+ REGISTER_BLTIN(run::arrayLength,"arrayLength");
+#line 751 "runarray.in"
+ REGISTER_BLTIN(run::arrayKeys,"arrayKeys");
+#line 766 "runarray.in"
+ REGISTER_BLTIN(run::arrayCyclicFlag,"arrayCyclicFlag");
+#line 773 "runarray.in"
+ REGISTER_BLTIN(run::arraySetCyclicFlag,"arraySetCyclicFlag");
+#line 780 "runarray.in"
+ REGISTER_BLTIN(run::arrayInitializedHelper,"arrayInitializedHelper");
+#line 791 "runarray.in"
+ REGISTER_BLTIN(run::arrayInitialized,"arrayInitialized");
+#line 797 "runarray.in"
+ REGISTER_BLTIN(run::arrayCyclicHelper,"arrayCyclicHelper");
+#line 804 "runarray.in"
+ REGISTER_BLTIN(run::arrayCyclic,"arrayCyclic");
+#line 810 "runarray.in"
+ REGISTER_BLTIN(run::arrayPushHelper,"arrayPushHelper");
+#line 818 "runarray.in"
+ REGISTER_BLTIN(run::arrayPush,"arrayPush");
+#line 824 "runarray.in"
+ REGISTER_BLTIN(run::arrayAppendHelper,"arrayAppendHelper");
+#line 833 "runarray.in"
+ REGISTER_BLTIN(run::arrayAppend,"arrayAppend");
+#line 839 "runarray.in"
+ REGISTER_BLTIN(run::arrayPopHelper,"arrayPopHelper");
+#line 848 "runarray.in"
+ REGISTER_BLTIN(run::arrayPop,"arrayPop");
+#line 854 "runarray.in"
+ REGISTER_BLTIN(run::arrayInsertHelper,"arrayInsertHelper");
+#line 865 "runarray.in"
+ REGISTER_BLTIN(run::arrayInsert,"arrayInsert");
+#line 871 "runarray.in"
+ REGISTER_BLTIN(run::arrayDelete,"arrayDelete");
+#line 877 "runarray.in"
+ REGISTER_BLTIN(run::arrayAlias,"arrayAlias");
+#line 882 "runarray.in"
+ REGISTER_BLTIN(run::arrayIntArray,"arrayIntArray");
+#line 900 "runarray.in"
+ addFunc(ve, run::gen_runarray32, IntArray(), SYM(complement), formal(IntArray(), SYM(a), false, false), formal(primInt(), SYM(n), false, false));
+#line 919 "runarray.in"
+ REGISTER_BLTIN(run::arraySequence,"arraySequence");
+#line 932 "runarray.in"
+ addFunc(ve, run::gen_runarray34, IntArray(), SYM(sequence), formal(primInt(), SYM(n), false, false));
+#line 943 "runarray.in"
+ REGISTER_BLTIN(run::arrayFunction,"arrayFunction");
+#line 956 "runarray.in"
+ REGISTER_BLTIN(run::arraySort,"arraySort");
+#line 966 "runarray.in"
+ REGISTER_BLTIN(run::arraySearch,"arraySearch");
+#line 985 "runarray.in"
+ addFunc(ve, run::gen_runarray38, primBoolean(), SYM(all), formal(booleanArray(), SYM(a), false, false));
+#line 994 "runarray.in"
+ addFunc(ve, run::gen_runarray39, booleanArray(), SYM_LOGNOT, formal(booleanArray(), SYM(a), false, false));
+#line 1003 "runarray.in"
+ addFunc(ve, run::gen_runarray40, primInt(), SYM(sum), formal(booleanArray(), SYM(a), false, false));
+#line 1012 "runarray.in"
+ REGISTER_BLTIN(run::arrayConcat,"arrayConcat");
+#line 1040 "runarray.in"
+ REGISTER_BLTIN(run::array2Transpose,"array2Transpose");
+#line 1066 "runarray.in"
+ REGISTER_BLTIN(run::array3Transpose,"array3Transpose");
+#line 1153 "runarray.in"
+ addFunc(ve, run::gen_runarray44, primInt(), SYM(find), formal(booleanArray(), SYM(a), false, false), formal(primInt(), SYM(n), true, false));
+#line 1172 "runarray.in"
+ addFunc(ve, run::gen_runarray45, IntArray(), SYM(findall), formal(booleanArray(), SYM(a), false, false));
+#line 1185 "runarray.in"
+ REGISTER_BLTIN(run::arrayConditional,"arrayConditional");
+#line 1211 "runarray.in"
+ addFunc(ve, run::gen_runarray47, realArray2(), SYM(identity), formal(primInt(), SYM(n), false, false));
+#line 1217 "runarray.in"
+ addFunc(ve, run::gen_runarray48, realArray2(), SYM(inverse), formal(realArray2(), SYM(a), false, false));
+#line 1229 "runarray.in"
+ addFunc(ve, run::gen_runarray49, realArray(), SYM(solve), formal(realArray2(), SYM(a), false, false), formal(realArray(), SYM(b), false, false), formal(primBoolean(), SYM(warn), true, false));
+#line 1282 "runarray.in"
+ addFunc(ve, run::gen_runarray50, realArray2(), SYM(solve), formal(realArray2(), SYM(a), false, false), formal(realArray2(), SYM(b), false, false), formal(primBoolean(), SYM(warn), true, false));
+#line 1347 "runarray.in"
+ addFunc(ve, run::gen_runarray51, primReal(), SYM(determinant), formal(realArray2(), SYM(a), false, false));
+#line 1364 "runarray.in"
+ addFunc(ve, run::gen_runarray52, realArray(), SYM_TIMES, formal(realArray2(), SYM(a), false, false), formal(realArray(), SYM(b), false, false));
+#line 1383 "runarray.in"
+ addFunc(ve, run::gen_runarray53, realArray(), SYM_TIMES, formal(realArray(), SYM(a), false, false), formal(realArray2(), SYM(b), false, false));
+#line 1412 "runarray.in"
+ addFunc(ve, run::gen_runarray54, IntArray2(), SYM_TIMES, formal(IntArray2(), SYM(a), false, false), formal(IntArray2(), SYM(b), false, false));
+#line 1417 "runarray.in"
+ addFunc(ve, run::gen_runarray55, realArray2(), SYM_TIMES, formal(realArray2(), SYM(a), false, false), formal(realArray2(), SYM(b), false, false));
+#line 1422 "runarray.in"
+ addFunc(ve, run::gen_runarray56, pairArray2(), SYM_TIMES, formal(pairArray2(), SYM(a), false, false), formal(pairArray2(), SYM(b), false, false));
+#line 1427 "runarray.in"
+ addFunc(ve, run::gen_runarray57, primTriple(), SYM_TIMES, formal(realArray2(), SYM(t), false, false), formal(primTriple(), SYM(v), false, false));
+#line 1432 "runarray.in"
+ addFunc(ve, run::gen_runarray58, realArray2(), SYM(AtA), formal(realArray2(), SYM(a), false, false));
+#line 1437 "runarray.in"
+ addFunc(ve, run::gen_runarray59, primPair(), SYM(project), formal(primTriple(), SYM(v), false, false), formal(realArray2(), SYM(t), false, false));
+#line 1462 "runarray.in"
+ addFunc(ve, run::gen_runarray60, primReal(), SYM(dot), formal(realArray(), SYM(a), false, false), formal(realArray(), SYM(b), false, false));
+#line 1472 "runarray.in"
+ addFunc(ve, run::gen_runarray61, primPair(), SYM(dot), formal(pairArray(), SYM(a), false, false), formal(pairArray(), SYM(b), false, false));
+#line 1482 "runarray.in"
+ addFunc(ve, run::gen_runarray62, realArray(), SYM(tridiagonal), formal(realArray(), SYM(a), false, false), formal(realArray(), SYM(b), false, false), formal(realArray(), SYM(c), false, false), formal(realArray(), SYM(f), false, false));
+#line 1586 "runarray.in"
+ addFunc(ve, run::gen_runarray63, primReal(), SYM(newton), formal(primInt(), SYM(iterations), true, false), formal(realRealFunction(), SYM(f), false, false), formal(realRealFunction(), SYM(fprime), false, false), formal(primReal(), SYM(x), false, false), formal(primBoolean(), SYM(verbose), true, false));
+#line 1633 "runarray.in"
+ addFunc(ve, run::gen_runarray64, primReal(), SYM(newton), formal(primInt(), SYM(iterations), true, false), formal(realRealFunction(), SYM(f), false, false), formal(realRealFunction(), SYM(fprime), false, false), formal(primReal(), SYM(x1), false, false), formal(primReal(), SYM(x2), false, false), formal(primBoolean(), SYM(verbose), true, false));
+#line 1715 "runarray.in"
+ addFunc(ve, run::gen_runarray65, primReal(), SYM(_findroot), formal(realRealFunction(), SYM(f), false, false), formal(primReal(), SYM(a), false, false), formal(primReal(), SYM(b), false, false), formal(primReal(), SYM(tolerance), false, false), formal(primReal(), SYM(fa), false, false), formal(primReal(), SYM(fb), false, false));
+#line 1817 "runarray.in"
+ addFunc(ve, run::gen_runarray66, primReal(), SYM(simpson), formal(realRealFunction(), SYM(f), false, false), formal(primReal(), SYM(a), false, false), formal(primReal(), SYM(b), false, false), formal(primReal(), SYM(acc), true, false), formal(primReal(), SYM(dxmax), true, false));
+#line 1831 "runarray.in"
+ addFunc(ve, run::gen_runarray67, pairArray(), SYM(fft), formal(pairArray(), SYM(a), false, false), formal(primInt(), SYM(sign), true, false));
+#line 1862 "runarray.in"
+ addFunc(ve, run::gen_runarray68, pairArray2(), SYM(fft), formal(pairArray2(), SYM(a), false, false), formal(primInt(), SYM(sign), true, false));
+#line 1908 "runarray.in"
+ addFunc(ve, run::gen_runarray69, pairArray3(), SYM(fft), formal(pairArray3(), SYM(a), false, false), formal(primInt(), SYM(sign), true, false));
+#line 1967 "runarray.in"
+ addFunc(ve, run::gen_runarray70, IntArray2(), SYM(triangulate), formal(pairArray(), SYM(z), false, false));
+#line 2002 "runarray.in"
+ addFunc(ve, run::gen_runarray71, primReal(), SYM(norm), formal(realArray(), SYM(a), false, false));
+#line 2013 "runarray.in"
+ addFunc(ve, run::gen_runarray72, primReal(), SYM(norm), formal(realArray2(), SYM(a), false, false));
+#line 2028 "runarray.in"
+ addFunc(ve, run::gen_runarray73, primReal(), SYM(norm), formal(tripleArray2(), SYM(a), false, false));
+#line 2043 "runarray.in"
+ addFunc(ve, run::gen_runarray74, primReal(), SYM(change2), formal(tripleArray2(), SYM(a), false, false));
+#line 2065 "runarray.in"
+ addFunc(ve, run::gen_runarray75, primTriple(), SYM(minbezier), formal(tripleArray2(), SYM(p), false, false), formal(primTriple(), SYM(b), false, false));
+#line 2077 "runarray.in"
+ addFunc(ve, run::gen_runarray76, primTriple(), SYM(maxbezier), formal(tripleArray2(), SYM(p), false, false), formal(primTriple(), SYM(b), false, false));
+#line 2089 "runarray.in"
+ addFunc(ve, run::gen_runarray77, primPair(), SYM(minratio), formal(tripleArray2(), SYM(p), false, false), formal(primPair(), SYM(b), false, false));
+#line 2101 "runarray.in"
+ addFunc(ve, run::gen_runarray78, primPair(), SYM(maxratio), formal(tripleArray2(), SYM(p), false, false), formal(primPair(), SYM(b), false, false));
+#line 2113 "runarray.in"
+ addFunc(ve, run::gen_runarray79, realArray(), SYM(_projection));
+}
+
+} // namespace trans