/***** Autogenerated from runarray.in; changes will be overwritten *****/

/***** 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 23 "./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

void unused(void *);

namespace run {
typedef double real;

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 23 "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, then ./configure; make";
#endif

#ifdef HAVE_EIGEN_DENSE
#include <Eigen/Dense>
typedef std::complex<double> Complex;
static const char *square="matrix must be square";
using Eigen::MatrixXd;
using Eigen::MatrixXcd;
using Eigen::RealSchur;
using Eigen::ComplexSchur;
#else
static const char *installEIGEN=
"Please install eigen3, then ./configure; make";
#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 realarray3;
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::realArray3;
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 610 "./runarray.in"
void emptyArray(stack *Stack)
{
#line 611 "./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 620 "./runarray.in"
void newDeepArray(stack *Stack)
{
Int depth=vm::pop<Int>(Stack);
#line 621 "./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 640 "./runarray.in"
void newInitializedArray(stack *Stack)
{
Int n=vm::pop<Int>(Stack);
#line 641 "./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 654 "./runarray.in"
void newAppendedArray(stack *Stack)
{
Int n=vm::pop<Int>(Stack);
array* tail=vm::pop<array*>(Stack);
#line 655 "./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 672 "./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 673 "./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 684 "./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 685 "./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 693 "./runarray.in"
void arrayRead(stack *Stack)
{
Int n=vm::pop<Int>(Stack);
array * a=vm::pop<array *>(Stack);
#line 694 "./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 705 "./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 706 "./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 713 "./runarray.in"
void arraySliceReadToEnd(stack *Stack)
{
Int left=vm::pop<Int>(Stack);
array * a=vm::pop<array *>(Stack);
#line 714 "./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 721 "./runarray.in"
void arrayArrayRead(stack *Stack)
{
Int n=vm::pop<Int>(Stack);
array * a=vm::pop<array *>(Stack);
#line 722 "./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 730 "./runarray.in"
void arrayWrite(stack *Stack)
{
item value=vm::pop(Stack);
Int n=vm::pop<Int>(Stack);
array * a=vm::pop<array *>(Stack);
#line 731 "./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 745 "./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 746 "./runarray.in"
checkArray(src);
checkArray(dest);
dest->setSlice(left, right, src);
{Stack->push<array*>(src); return;}
}

#line 753 "./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 754 "./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 762 "./runarray.in"
void arrayLength(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
#line 763 "./runarray.in"
{Stack->push<Int>((Int) checkArray(a)); return;}
}

// Returns an array of integers representing the keys of the array.
#line 768 "./runarray.in"
void arrayKeys(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
#line 769 "./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 783 "./runarray.in"
void arrayCyclicFlag(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
#line 784 "./runarray.in"
checkArray(a);
{Stack->push<bool>(a->cyclic()); return;}
}

#line 789 "./runarray.in"
void arraySetCyclicFlag(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
bool b=vm::pop<bool>(Stack);
#line 790 "./runarray.in"
checkArray(a);
a->cyclic(b);
{Stack->push<bool>(b); return;}
}

// Check to see if an array element is initialized.
#line 797 "./runarray.in"
void arrayInitializedHelper(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
Int n=vm::pop<Int>(Stack);
#line 798 "./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 808 "./runarray.in"
void arrayInitialized(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
#line 809 "./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 814 "./runarray.in"
void arrayCyclicHelper(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
bool b=vm::pop<bool>(Stack);
#line 815 "./runarray.in"
checkArray(a);
a->cyclic(b);
}

// Set the cyclic flag for an array.
#line 821 "./runarray.in"
void arrayCyclic(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
#line 822 "./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 827 "./runarray.in"
void arrayPushHelper(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
item x=vm::pop(Stack);
#line 828 "./runarray.in"
checkArray(a);
a->push(x);
{Stack->push(x); return;}
}

// Returns the push method for an array.
#line 835 "./runarray.in"
void arrayPush(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
#line 836 "./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 841 "./runarray.in"
void arrayAppendHelper(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
array * b=vm::pop<array *>(Stack);
#line 842 "./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 850 "./runarray.in"
void arrayAppend(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
#line 851 "./runarray.in"
{Stack->push<callable*>(new thunk(new bfunc(arrayAppendHelper),a)); return;}
}

// The helper function for the pop method.
#line 856 "./runarray.in"
void arrayPopHelper(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
#line 857 "./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 865 "./runarray.in"
void arrayPop(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
#line 866 "./runarray.in"
{Stack->push<callable*>(new thunk(new bfunc(arrayPopHelper),a)); return;}
}

// The helper function for the insert method.
#line 871 "./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 872 "./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 882 "./runarray.in"
void arrayInsert(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
#line 883 "./runarray.in"
{Stack->push<callable*>(new thunk(new bfunc(arrayInsertHelper),a)); return;}
}

// Returns the delete method for an array.
#line 888 "./runarray.in"
void arrayDelete(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
#line 889 "./runarray.in"
{Stack->push<callable*>(new thunk(new bfunc(arrayDeleteHelper),a)); return;}
}

#line 893 "./runarray.in"
void arrayAlias(stack *Stack)
{
array * b=vm::pop<array *>(Stack);
array * a=vm::pop<array *>(Stack);
#line 894 "./runarray.in"
{Stack->push<bool>(a==b); return;}
}

// Return array formed by indexing array a with elements of integer array b
#line 899 "./runarray.in"
void arrayIntArray(stack *Stack)
{
array * b=vm::pop<array *>(Stack);
array * a=vm::pop<array *>(Stack);
#line 900 "./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 918 "./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 919 "./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 936 "./runarray.in"
void arraySequence(stack *Stack)
{
Int n=vm::pop<Int>(Stack);
callable * f=vm::pop<callable *>(Stack);
#line 937 "./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 949 "./runarray.in"
// Intarray* sequence(Int n);
void gen_runarray34(stack *Stack)
{
Int n=vm::pop<Int>(Stack);
#line 950 "./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 960 "./runarray.in"
void arrayFunction(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
callable * f=vm::pop<callable *>(Stack);
#line 961 "./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 972 "./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 973 "./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 982 "./runarray.in"
void arraySearch(stack *Stack)
{
callable * less=vm::pop<callable *>(Stack);
item key=vm::pop(Stack);
array * a=vm::pop<array *>(Stack);
#line 983 "./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 1001 "./runarray.in"
// bool all(boolarray *a);
void gen_runarray38(stack *Stack)
{
boolarray * a=vm::pop<boolarray *>(Stack);
#line 1002 "./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 1010 "./runarray.in"
// boolarray* !(boolarray* a);
void gen_runarray39(stack *Stack)
{
boolarray* a=vm::pop<boolarray*>(Stack);
#line 1011 "./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 1019 "./runarray.in"
// Int sum(boolarray *a);
void gen_runarray40(stack *Stack)
{
boolarray * a=vm::pop<boolarray *>(Stack);
#line 1020 "./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 1028 "./runarray.in"
void arrayConcat(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
#line 1029 "./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 1056 "./runarray.in"
void array2Transpose(stack *Stack)
{
array * a=vm::pop<array *>(Stack);
#line 1057 "./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 1086 "./runarray.in"
void array3Transpose(stack *Stack)
{
array * perm=vm::pop<array *>(Stack);
array * a=vm::pop<array *>(Stack);
#line 1087 "./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 1171 "./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 1172 "./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 1189 "./runarray.in"
// Intarray* findall(boolarray *a);
void gen_runarray45(stack *Stack)
{
boolarray * a=vm::pop<boolarray *>(Stack);
#line 1190 "./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 1203 "./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 1204 "./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 1228 "./runarray.in"
// realarray2* identity(Int n);
void gen_runarray47(stack *Stack)
{
Int n=vm::pop<Int>(Stack);
#line 1229 "./runarray.in"
{Stack->push<realarray2*>(Identity(n)); return;}
}

// Return the inverse of an n x n matrix a using Gauss-Jordan elimination.
#line 1234 "./runarray.in"
// realarray2* inverse(realarray2 *a);
void gen_runarray48(stack *Stack)
{
realarray2 * a=vm::pop<realarray2 *>(Stack);
#line 1235 "./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 1248 "./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 1249 "./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 1301 "./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 1302 "./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 1364 "./runarray.in"
// real determinant(realarray2 *a);
void gen_runarray51(stack *Stack)
{
realarray2 * a=vm::pop<realarray2 *>(Stack);
#line 1365 "./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 1380 "./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 1381 "./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 1399 "./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 1400 "./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 1428 "./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 1429 "./runarray.in"
{Stack->push<Intarray2*>(mult<Int>(a,b)); return;}
}

#line 1433 "./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 1434 "./runarray.in"
{Stack->push<realarray2*>(mult<real>(a,b)); return;}
}

#line 1438 "./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 1439 "./runarray.in"
{Stack->push<pairarray2*>(mult<pair>(a,b)); return;}
}

#line 1443 "./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 1444 "./runarray.in"
{Stack->push<triple>(*t*v); return;}
}

#line 1448 "./runarray.in"
// realarray2* AtA(realarray2 *a);
void gen_runarray58(stack *Stack)
{
realarray2 * a=vm::pop<realarray2 *>(Stack);
#line 1449 "./runarray.in"
{Stack->push<realarray2*>(AtA<real>(a)); return;}
}

#line 1453 "./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 1454 "./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 1479 "./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 1480 "./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 1489 "./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 1490 "./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 1505 "./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 1506 "./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 1603 "./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 1605 "./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 1651 "./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 1653 "./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 1739 "./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 1741 "./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 1833 "./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 1835 "./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 1848 "./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 1849 "./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 1879 "./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 1880 "./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 1925 "./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 1926 "./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;}
}

// Compute the real Schur decomposition of a 2D pair array
#line 1984 "./runarray.in"
// realarray3* _schur(realarray2 *a);
void gen_runarray70(stack *Stack)
{
realarray2 * a=vm::pop<realarray2 *>(Stack);
#line 1985 "./runarray.in"
#ifdef HAVE_EIGEN_DENSE
size_t n=checkArray(a);

MatrixXd A(n,n);
RealSchur<MatrixXd> schur(n);

array *S=new array(2);

if(n) {
for(size_t i=0; i < n; ++i) {
array *ai=read<array *>(a,i);
size_t aisize=checkArray(ai);
if(aisize != n) error(square);
for(size_t j=0; j < n; ++j)
A(i,j)=read<double>(ai,j);
}

schur.compute(A);
MatrixXd U=schur.matrixU();
MatrixXd T=schur.matrixT();

array *u=new array(n);
array *t=new array(n);
(*S)[0]=u;
(*S)[1]=t;

for(size_t i=0; i < n; ++i) {
array *ui=new array(n);
array *ti=new array(n);
(*u)[i]=ui;
(*t)[i]=ti;
for(size_t j=0; j < n; ++j) {
(*ui)[j]=U(i,j);
(*ti)[j]=T(i,j);
}
}
}
#else
unused(a);
array *S=new array(0);
error(installEIGEN);
#endif // HAVE_EIGEN_DENSE
{Stack->push<realarray3*>(S); return;}
}

// Compute the Schur decomposition of a 2D pair array
#line 2032 "./runarray.in"
// pairarray3* _schur(pairarray2 *a);
void gen_runarray71(stack *Stack)
{
pairarray2 * a=vm::pop<pairarray2 *>(Stack);
#line 2033 "./runarray.in"
#ifdef HAVE_EIGEN_DENSE
size_t n=checkArray(a);

MatrixXcd A(n,n);
ComplexSchur<MatrixXcd> schur(n);

array *S=new array(2);

if(n) {
for(size_t i=0; i < n; ++i) {
array *ai=read<array *>(a,i);
size_t aisize=checkArray(ai);
if(aisize != n) error(square);
for(size_t j=0; j < n; ++j) {
pair z=read<pair>(ai,j);
A(i,j)=Complex(z.getx(),z.gety());
}
}

schur.compute(A);
MatrixXcd U=schur.matrixU();
MatrixXcd T=schur.matrixT();

array *u=new array(n);
array *t=new array(n);
(*S)[0]=u;
(*S)[1]=t;

for(size_t i=0; i < n; ++i) {
array *ui=new array(n);
array *ti=new array(n);
(*u)[i]=ui;
(*t)[i]=ti;
for(size_t j=0; j < n; ++j) {
Complex z=U(i,j);
Complex w=T(i,j);
(*ui)[j]=pair(z.real(),z.imag());
(*ti)[j]=pair(w.real(),w.imag());
}
}
}
#else
unused(a);
array *S=new array(0);
error(installEIGEN);
#endif // HAVE_EIGEN_DENSE
{Stack->push<pairarray3*>(S); return;}
}

#line 2083 "./runarray.in"
// Intarray2* triangulate(pairarray *z);
void gen_runarray72(stack *Stack)
{
pairarray * z=vm::pop<pairarray *>(Stack);
#line 2084 "./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 2118 "./runarray.in"
// real norm(realarray *a);
void gen_runarray73(stack *Stack)
{
realarray * a=vm::pop<realarray *>(Stack);
#line 2119 "./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 2129 "./runarray.in"
// real norm(realarray2 *a);
void gen_runarray74(stack *Stack)
{
realarray2 * a=vm::pop<realarray2 *>(Stack);
#line 2130 "./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 2144 "./runarray.in"
// real norm(triplearray2 *a);
void gen_runarray75(stack *Stack)
{
triplearray2 * a=vm::pop<triplearray2 *>(Stack);
#line 2145 "./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 2159 "./runarray.in"
// real change2(triplearray2 *a);
void gen_runarray76(stack *Stack)
{
triplearray2 * a=vm::pop<triplearray2 *>(Stack);
#line 2160 "./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 2181 "./runarray.in"
// triple minbezier(triplearray2 *P, triple b);
void gen_runarray77(stack *Stack)
{
triple b=vm::pop<triple>(Stack);
triplearray2 * P=vm::pop<triplearray2 *>(Stack);
#line 2182 "./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 2193 "./runarray.in"
// triple maxbezier(triplearray2 *P, triple b);
void gen_runarray78(stack *Stack)
{
triple b=vm::pop<triple>(Stack);
triplearray2 * P=vm::pop<triplearray2 *>(Stack);
#line 2194 "./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 2205 "./runarray.in"
// pair minratio(triplearray2 *P, pair b);
void gen_runarray79(stack *Stack)
{
pair b=vm::pop<pair>(Stack);
triplearray2 * P=vm::pop<triplearray2 *>(Stack);
#line 2206 "./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 2217 "./runarray.in"
// pair maxratio(triplearray2 *P, pair b);
void gen_runarray80(stack *Stack)
{
pair b=vm::pop<pair>(Stack);
triplearray2 * P=vm::pop<triplearray2 *>(Stack);
#line 2218 "./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 2229 "./runarray.in"
// realarray* _projection();
void gen_runarray81(stack *Stack)
{
#line 2230 "./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 609 "./runarray.in"
REGISTER_BLTIN(run::emptyArray,"emptyArray");
#line 615 "./runarray.in"
REGISTER_BLTIN(run::newDeepArray,"newDeepArray");
#line 637 "./runarray.in"
REGISTER_BLTIN(run::newInitializedArray,"newInitializedArray");
#line 652 "./runarray.in"
REGISTER_BLTIN(run::newAppendedArray,"newAppendedArray");
#line 668 "./runarray.in"
REGISTER_BLTIN(run::copyArrayValue,"copyArrayValue");
#line 680 "./runarray.in"
REGISTER_BLTIN(run::copyArray,"copyArray");
#line 692 "./runarray.in"
REGISTER_BLTIN(run::arrayRead,"arrayRead");
#line 704 "./runarray.in"
REGISTER_BLTIN(run::arraySliceRead,"arraySliceRead");
#line 711 "./runarray.in"
REGISTER_BLTIN(run::arraySliceReadToEnd,"arraySliceReadToEnd");
#line 719 "./runarray.in"
REGISTER_BLTIN(run::arrayArrayRead,"arrayArrayRead");
#line 728 "./runarray.in"
REGISTER_BLTIN(run::arrayWrite,"arrayWrite");
#line 745 "./runarray.in"
REGISTER_BLTIN(run::arraySliceWrite,"arraySliceWrite");
#line 753 "./runarray.in"
REGISTER_BLTIN(run::arraySliceWriteToEnd,"arraySliceWriteToEnd");
#line 761 "./runarray.in"
REGISTER_BLTIN(run::arrayLength,"arrayLength");
#line 767 "./runarray.in"
REGISTER_BLTIN(run::arrayKeys,"arrayKeys");
#line 782 "./runarray.in"
REGISTER_BLTIN(run::arrayCyclicFlag,"arrayCyclicFlag");
#line 789 "./runarray.in"
REGISTER_BLTIN(run::arraySetCyclicFlag,"arraySetCyclicFlag");
#line 796 "./runarray.in"
REGISTER_BLTIN(run::arrayInitializedHelper,"arrayInitializedHelper");
#line 807 "./runarray.in"
REGISTER_BLTIN(run::arrayInitialized,"arrayInitialized");
#line 813 "./runarray.in"
REGISTER_BLTIN(run::arrayCyclicHelper,"arrayCyclicHelper");
#line 820 "./runarray.in"
REGISTER_BLTIN(run::arrayCyclic,"arrayCyclic");
#line 826 "./runarray.in"
REGISTER_BLTIN(run::arrayPushHelper,"arrayPushHelper");
#line 834 "./runarray.in"
REGISTER_BLTIN(run::arrayPush,"arrayPush");
#line 840 "./runarray.in"
REGISTER_BLTIN(run::arrayAppendHelper,"arrayAppendHelper");
#line 849 "./runarray.in"
REGISTER_BLTIN(run::arrayAppend,"arrayAppend");
#line 855 "./runarray.in"
REGISTER_BLTIN(run::arrayPopHelper,"arrayPopHelper");
#line 864 "./runarray.in"
REGISTER_BLTIN(run::arrayPop,"arrayPop");
#line 870 "./runarray.in"
REGISTER_BLTIN(run::arrayInsertHelper,"arrayInsertHelper");
#line 881 "./runarray.in"
REGISTER_BLTIN(run::arrayInsert,"arrayInsert");
#line 887 "./runarray.in"
REGISTER_BLTIN(run::arrayDelete,"arrayDelete");
#line 893 "./runarray.in"
REGISTER_BLTIN(run::arrayAlias,"arrayAlias");
#line 898 "./runarray.in"
REGISTER_BLTIN(run::arrayIntArray,"arrayIntArray");
#line 916 "./runarray.in"
addFunc(ve, run::gen_runarray32, IntArray(), SYM(complement), formal(IntArray(), SYM(a), false, false), formal(primInt(), SYM(n), false, false));
#line 935 "./runarray.in"
REGISTER_BLTIN(run::arraySequence,"arraySequence");
#line 948 "./runarray.in"
addFunc(ve, run::gen_runarray34, IntArray(), SYM(sequence), formal(primInt(), SYM(n), false, false));
#line 959 "./runarray.in"
REGISTER_BLTIN(run::arrayFunction,"arrayFunction");
#line 972 "./runarray.in"
REGISTER_BLTIN(run::arraySort,"arraySort");
#line 982 "./runarray.in"
REGISTER_BLTIN(run::arraySearch,"arraySearch");
#line 1001 "./runarray.in"
addFunc(ve, run::gen_runarray38, primBoolean(), SYM(all), formal(booleanArray(), SYM(a), false, false));
#line 1010 "./runarray.in"
addFunc(ve, run::gen_runarray39, booleanArray(), SYM_LOGNOT, formal(booleanArray(), SYM(a), false, false));
#line 1019 "./runarray.in"
addFunc(ve, run::gen_runarray40, primInt(), SYM(sum), formal(booleanArray(), SYM(a), false, false));
#line 1028 "./runarray.in"
REGISTER_BLTIN(run::arrayConcat,"arrayConcat");
#line 1056 "./runarray.in"
REGISTER_BLTIN(run::array2Transpose,"array2Transpose");
#line 1082 "./runarray.in"
REGISTER_BLTIN(run::array3Transpose,"array3Transpose");
#line 1169 "./runarray.in"
addFunc(ve, run::gen_runarray44, primInt(), SYM(find), formal(booleanArray(), SYM(a), false, false), formal(primInt(), SYM(n), true, false));
#line 1188 "./runarray.in"
addFunc(ve, run::gen_runarray45, IntArray(), SYM(findall), formal(booleanArray(), SYM(a), false, false));
#line 1201 "./runarray.in"
REGISTER_BLTIN(run::arrayConditional,"arrayConditional");
#line 1227 "./runarray.in"
addFunc(ve, run::gen_runarray47, realArray2(), SYM(identity), formal(primInt(), SYM(n), false, false));
#line 1233 "./runarray.in"
addFunc(ve, run::gen_runarray48, realArray2(), SYM(inverse), formal(realArray2(), SYM(a), false, false));
#line 1245 "./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 1298 "./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 1363 "./runarray.in"
addFunc(ve, run::gen_runarray51, primReal(), SYM(determinant), formal(realArray2(), SYM(a), false, false));
#line 1380 "./runarray.in"
addFunc(ve, run::gen_runarray52, realArray(), SYM_TIMES, formal(realArray2(), SYM(a), false, false), formal(realArray(), SYM(b), false, false));
#line 1399 "./runarray.in"
addFunc(ve, run::gen_runarray53, realArray(), SYM_TIMES, formal(realArray(), SYM(a), false, false), formal(realArray2(), SYM(b), false, false));
#line 1428 "./runarray.in"
addFunc(ve, run::gen_runarray54, IntArray2(), SYM_TIMES, formal(IntArray2(), SYM(a), false, false), formal(IntArray2(), SYM(b), false, false));
#line 1433 "./runarray.in"
addFunc(ve, run::gen_runarray55, realArray2(), SYM_TIMES, formal(realArray2(), SYM(a), false, false), formal(realArray2(), SYM(b), false, false));
#line 1438 "./runarray.in"
addFunc(ve, run::gen_runarray56, pairArray2(), SYM_TIMES, formal(pairArray2(), SYM(a), false, false), formal(pairArray2(), SYM(b), false, false));
#line 1443 "./runarray.in"
addFunc(ve, run::gen_runarray57, primTriple(), SYM_TIMES, formal(realArray2(), SYM(t), false, false), formal(primTriple(), SYM(v), false, false));
#line 1448 "./runarray.in"
addFunc(ve, run::gen_runarray58, realArray2(), SYM(AtA), formal(realArray2(), SYM(a), false, false));
#line 1453 "./runarray.in"
addFunc(ve, run::gen_runarray59, primPair(), SYM(project), formal(primTriple(), SYM(v), false, false), formal(realArray2(), SYM(t), false, false));
#line 1478 "./runarray.in"
addFunc(ve, run::gen_runarray60, primReal(), SYM(dot), formal(realArray(), SYM(a), false, false), formal(realArray(), SYM(b), false, false));
#line 1488 "./runarray.in"
addFunc(ve, run::gen_runarray61, primPair(), SYM(dot), formal(pairArray(), SYM(a), false, false), formal(pairArray(), SYM(b), false, false));
#line 1498 "./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 1602 "./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 1649 "./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 1731 "./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 1833 "./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 1847 "./runarray.in"
addFunc(ve, run::gen_runarray67, pairArray(), SYM(fft), formal(pairArray(), SYM(a), false, false), formal(primInt(), SYM(sign), true, false));
#line 1878 "./runarray.in"
addFunc(ve, run::gen_runarray68, pairArray2(), SYM(fft), formal(pairArray2(), SYM(a), false, false), formal(primInt(), SYM(sign), true, false));
#line 1924 "./runarray.in"
addFunc(ve, run::gen_runarray69, pairArray3(), SYM(fft), formal(pairArray3(), SYM(a), false, false), formal(primInt(), SYM(sign), true, false));
#line 1983 "./runarray.in"
addFunc(ve, run::gen_runarray70, realArray3(), SYM(_schur), formal(realArray2(), SYM(a), false, false));
#line 2031 "./runarray.in"
addFunc(ve, run::gen_runarray71, pairArray3(), SYM(_schur), formal(pairArray2(), SYM(a), false, false));
#line 2083 "./runarray.in"
addFunc(ve, run::gen_runarray72, IntArray2(), SYM(triangulate), formal(pairArray(), SYM(z), false, false));
#line 2118 "./runarray.in"
addFunc(ve, run::gen_runarray73, primReal(), SYM(norm), formal(realArray(), SYM(a), false, false));
#line 2129 "./runarray.in"
addFunc(ve, run::gen_runarray74, primReal(), SYM(norm), formal(realArray2(), SYM(a), false, false));
#line 2144 "./runarray.in"
addFunc(ve, run::gen_runarray75, primReal(), SYM(norm), formal(tripleArray2(), SYM(a), false, false));
#line 2159 "./runarray.in"
addFunc(ve, run::gen_runarray76, primReal(), SYM(change2), formal(tripleArray2(), SYM(a), false, false));
#line 2181 "./runarray.in"
addFunc(ve, run::gen_runarray77, primTriple(), SYM(minbezier), formal(tripleArray2(), SYM(p), false, false), formal(primTriple(), SYM(b), false, false));
#line 2193 "./runarray.in"
addFunc(ve, run::gen_runarray78, primTriple(), SYM(maxbezier), formal(tripleArray2(), SYM(p), false, false), formal(primTriple(), SYM(b), false, false));
#line 2205 "./runarray.in"
addFunc(ve, run::gen_runarray79, primPair(), SYM(minratio), formal(tripleArray2(), SYM(p), false, false), formal(primPair(), SYM(b), false, false));
#line 2217 "./runarray.in"
addFunc(ve, run::gen_runarray80, primPair(), SYM(maxratio), formal(tripleArray2(), SYM(p), false, false), formal(primPair(), SYM(b), false, false));
#line 2229 "./runarray.in"
addFunc(ve, run::gen_runarray81, realArray(), SYM(_projection));
}

} // namespace trans