/*****
 * runarray.in
 *
 * Runtime functions for array operations.
 *
 *****/

pair     => primPair()
triple   => primTriple()
boolarray* => booleanArray()
Intarray*  => IntArray()
Intarray2*  => IntArray2()
realarray* => realArray()
realarray2* => realArray2()
pairarray* => pairArray()
triplearray2* => tripleArray2()
callableReal* => realRealFunction()


#include "array.h"
#include "arrayop.h"
#include "triple.h"
#include "path3.h"
#include "Delaunay.h"
#include "glrender.h"

#ifdef HAVE_LIBFFTW3
#include "fftw++.h"
#endif

using namespace camp;
using namespace vm;

typedef array boolarray;
typedef array Intarray;
typedef array Intarray2;
typedef array realarray;
typedef array realarray2;
typedef array pairarray;
typedef array triplearray2;

using types::booleanArray;
using types::IntArray;
using types::IntArray2;
using types::realArray;
using types::realArray2;
using types::pairArray;
using types::tripleArray2;

typedef callable callableReal;

void outOfBounds(const char *op, size_t len, Int n)
{
  ostringstream buf;
  buf << op << " array of length " << len << " with out-of-bounds index " << n;
  error(buf);
}

inline item& arrayRead(array *a, Int n)  
{
  size_t len=checkArray(a);
  bool cyclic=a->cyclic();
  if(cyclic && len > 0) n=imod(n,len);
  else if(n < 0 || n >= (Int) len) outOfBounds("reading",len,n);
  return (*a)[(unsigned) n];
}

// Helper function to create deep arrays.
static array* deepArray(Int depth, Int *dims)
{
  assert(depth > 0);
  
  if (depth == 1) {
    return new array(dims[0]);
  } else {
    Int length = dims[0];
    depth--; dims++;

    array *a = new array(length);

    for (Int index = 0; index < length; index++) {
      (*a)[index] = deepArray(depth, dims);
    }
    return a;
  }
}

namespace run {
array *Identity(Int n)
{
  size_t N=(size_t) n;
  array *c=new array(N);
  for(size_t i=0; i < N; ++i) {
    array *ci=new array(N);
    (*c)[i]=ci;
    for(size_t j=0; j < N; ++j)
      (*ci)[j]=0.0;
    (*ci)[i]=1.0;
  }
  return c;
}
}

static const char *incommensurate="Incommensurate matrices";
static const char *singular="Singular matrix";
static size_t *pivot,*Row,*Col;

namespace run {

array *copyArray(array *a)
{
  size_t size=checkArray(a);
  array *c=new array(size);
  for(size_t i=0; i < size; i++) 
    (*c)[i]=(*a)[i];
  return c;
}

inline size_t checkdimension(const array *a, size_t dim)
{
  size_t size=checkArray(a);
  if(dim && size != dim) {
    ostringstream buf;
    buf << "array of length " << dim << " expected";
    error(buf);
  }
  return size;
}

double *copyArrayC(const array *a, size_t dim, GCPlacement placement)
{
  size_t size=checkdimension(a,dim);
  double *c=(placement == NoGC) ? new double [size] : 
    new(placement) double[size];
  for(size_t i=0; i < size; i++) 
    c[i]=read<double>(a,i);
  return c;
}

triple *copyTripleArrayC(const array *a, size_t dim)
{
  size_t size=checkdimension(a,dim);
  triple *c=new triple[size];
  for(size_t i=0; i < size; i++) 
    c[i]=read<triple>(a,i);
  return c;
}

array *copyArray2(array *a)
{
  size_t size=checkArray(a);
  array *c=new array(size);
  for(size_t i=0; i < size; i++) {
    array *ai=read<array*>(a,i);
    size_t aisize=checkArray(ai);
    array *ci=new array(aisize);
    (*c)[i]=ci;
    for(size_t j=0; j < aisize; j++) 
      (*ci)[j]=(*ai)[j];
  }
  return c;
}

array *copyArray3(array *a)
{
  size_t size=checkArray(a);
  array *c=new array(size);
  for(size_t i=0; i < size; i++) {
    array *ai=read<array*>(a,i);
    size_t aisize=checkArray(ai);
    array *ci=new array(aisize);
    (*c)[i]=ci;
    for(size_t j=0; j < aisize; j++) {
      array *aij=read<array*>(ai,j);
      size_t aijsize=checkArray(aij);
      array *cij=new array(aijsize);
      (*ci)[j]=cij;
      for(size_t k=0; k < aijsize; k++) 
        (*cij)[k]=(*aij)[k];
    }
  }
  return c;
}

double *copyArray2C(const array *a, bool square, size_t dim2,
                    GCPlacement placement)
{
  size_t n=checkArray(a);
  size_t m=(square || n == 0) ? n : checkArray(read<array*>(a,0));
  if(n > 0 && dim2 && m != dim2) {
    ostringstream buf;
    buf << "second matrix dimension must be " << dim2;
    error(buf);
  }
  
  double *c=(placement == NoGC) ? new double [n*m] : new(placement) double[n*m];
  for(size_t i=0; i < n; i++) {
    array *ai=read<array*>(a,i);
    size_t aisize=checkArray(ai);
    if(aisize == m) {
      double *ci=c+i*m;
      for(size_t j=0; j < m; j++) 
        ci[j]=read<double>(ai,j);
    } else
      error(square ? "matrix must be square" : "matrix must be rectangular");
  }
  return c;
}

triple *copyTripleArray2C(const array *a, bool square, size_t dim2)
{
  size_t n=checkArray(a);
  size_t m=(square || n == 0) ? n : checkArray(read<array*>(a,0));
  if(n > 0 && dim2 && m != dim2) {
    ostringstream buf;
    buf << "second matrix dimension must be " << dim2;
    error(buf);
  }
  
  triple *c=new triple[n*m];
  for(size_t i=0; i < n; i++) {
    array *ai=read<array*>(a,i);
    size_t aisize=checkArray(ai);
    if(aisize == m) {
      triple *ci=c+i*m;
      for(size_t j=0; j < m; j++) 
        ci[j]=read<triple>(ai,j);
    } else
      error(square ? "matrix must be square" : "matrix must be rectangular");
  }
  return c;
}

double *copyTripleArray2Components(array *a, bool square, size_t dim2,
                                   GCPlacement placement)
{
  size_t n=checkArray(a);
  size_t m=(square || n == 0) ? n : checkArray(read<array*>(a,0));
  if(n > 0 && dim2 && m != dim2) {
    ostringstream buf;
    buf << "second matrix dimension must be " << dim2;
    error(buf);
  }
  
  size_t nm=n*m;
  double *cx=(placement == NoGC) ? new double [3*nm] :
    new(placement) double[3*nm];
  double *cy=cx+nm;
  double *cz=cx+2*nm;
  for(size_t i=0; i < n; i++) {
    array *ai=read<array*>(a,i);
    size_t aisize=checkArray(ai);
    if(aisize == m) {
      double *xi=cx+i*m;
      double *yi=cy+i*m;
      double *zi=cz+i*m;
      for(size_t j=0; j < m; j++) {
        triple v=read<triple>(ai,j);
        xi[j]=v.getx();
        yi[j]=v.gety();
        zi[j]=v.getz();
      }
    } else
      error(square ? "matrix must be square" : "matrix must be rectangular");
  }
  return cx;
}

triple operator *(const array& t, const triple& v)
{
  size_t n=checkArray(&t);
  if(n != 4) error(incommensurate);
  array *t0=read<array*>(t,0);
  array *t1=read<array*>(t,1);
  array *t2=read<array*>(t,2);
  array *t3=read<array*>(t,3);
  
  if(checkArray(t0) != 4 || checkArray(t1) != 4 || 
     checkArray(t2) != 4 || checkArray(t3) != 4)
    error(incommensurate);

  double x=v.getx();
  double y=v.gety();
  double z=v.getz();
  
  double f=read<real>(t3,0)*x+read<real>(t3,1)*y+read<real>(t3,2)*z+
    read<real>(t3,3);
  if(f == 0.0) run::dividebyzero();
  f=1.0/f;
  
  return triple((read<real>(t0,0)*x+read<real>(t0,1)*y+read<real>(t0,2)*z+
                 read<real>(t0,3))*f,
                (read<real>(t1,0)*x+read<real>(t1,1)*y+read<real>(t1,2)*z+
                 read<real>(t1,3))*f,
                (read<real>(t2,0)*x+read<real>(t2,1)*y+read<real>(t2,2)*z+
                 read<real>(t2,3))*f);
}

triple multshiftless(const array& t, const triple& v)
{
  size_t n=checkArray(&t);
  if(n != 4) error(incommensurate);
  array *t0=read<array*>(t,0);
  array *t1=read<array*>(t,1);
  array *t2=read<array*>(t,2);
  array *t3=read<array*>(t,3);
  
  if(checkArray(t0) != 4 || checkArray(t1) != 4 || 
     checkArray(t2) != 4 || checkArray(t3) != 4)
    error(incommensurate);

  double x=v.getx();
  double y=v.gety();
  double z=v.getz();
  
  double f=read<real>(t3,0)*x+read<real>(t3,1)*y+read<real>(t3,2)*z+
    read<real>(t3,3);
  if(f == 0.0) run::dividebyzero();
  f=1.0/f;
  
  return triple((read<real>(t0,0)*x+read<real>(t0,1)*y+read<real>(t0,2)*z)*f,
                (read<real>(t1,0)*x+read<real>(t1,1)*y+read<real>(t1,2)*z)*f,
                (read<real>(t2,0)*x+read<real>(t2,1)*y+read<real>(t2,2)*z)*f);
}

double norm(double *a, size_t n) 
{
  if(n == 0) return 0.0;
  double M=fabs(a[0]);
  for(size_t i=1; i < n; ++i)
    M=::max(M,fabs(a[i]));
  return M;
}

double norm(triple *a, size_t n) 
{
  if(n == 0) return 0.0;
  double M=a[0].abs2();
  for(size_t i=1; i < n; ++i)
    M=::max(M,a[i].abs2());
  return sqrt(M);
}

}

static inline void inverseAllocate(size_t n)
{
  pivot=new size_t[n];
  Row=new size_t[n];
  Col=new size_t[n];
}

static inline void inverseDeallocate()
{
  delete[] pivot;
  delete[] Row;
  delete[] Col;
}

callable *Func;
stack *FuncStack;
double wrapFunction(double x)
{
  FuncStack->push(x);
  Func->call(FuncStack);
  return pop<double>(FuncStack);
}

callable *compareFunc;
bool compareFunction(const vm::item& i, const vm::item& j)
{
  FuncStack->push(i);
  FuncStack->push(j);
  compareFunc->call(FuncStack);
  return pop<bool>(FuncStack);
}

void checkSquare(array *a) 
{
  size_t n=checkArray(a);
  for(size_t i=0; i < n; i++)
    if(checkArray(read<array*>(a,i)) != n)
      error("matrix a must be square");
}

// Crout's algorithm for computing the LU decomposition of a square matrix.
// cf. routine ludcmp (Press et al.,  Numerical Recipes, 1991).
Int LUdecompose(double *a, size_t n, size_t* index, bool warn=true)
{
  double *vv=new double[n];
  Int swap=1;
  for(size_t i=0; i < n; ++i) {
    double big=0.0;
    double *ai=a+i*n;
    for(size_t j=0; j < n; ++j) {
      double temp=fabs(ai[j]);
      if(temp > big) big=temp;
    }
    if(big == 0.0) {
      delete[] vv;
      if(warn) error(singular);
      else return 0;
    }
    vv[i]=1.0/big;
  }
  for(size_t j=0; j < n; ++j) {
    for(size_t i=0; i < j; ++i) {
      double *ai=a+i*n;
      double sum=ai[j];
      for(size_t k=0; k < i; ++k) {
        sum -= ai[k]*a[k*n+j];
      }
      ai[j]=sum;
    }
    double big=0.0;
    size_t imax=j;
    for(size_t i=j; i < n; ++i) {
      double *ai=a+i*n;
      double sum=ai[j];
      for(size_t k=0; k < j; ++k)
        sum -= ai[k]*a[k*n+j];
      ai[j]=sum;
      double temp=vv[i]*fabs(sum);
      if(temp >= big) {
        big=temp;
        imax=i;
      }
    }
    double *aj=a+j*n;
    double *aimax=a+imax*n;
    if(j != imax) {
      for(size_t k=0; k < n; ++k) {
        double temp=aimax[k];
        aimax[k]=aj[k];
        aj[k]=temp;
      }
      swap *= -1;
      vv[imax]=vv[j];
    }
    if(index) 
      index[j]=imax;
    if(j != n) {
      double denom=aj[j];
      if(denom == 0.0) {
        delete[] vv;
        if(warn) error(singular);
        else return 0;
      }
      for(size_t i=j+1; i < n; ++i)
        a[i*n+j] /= denom;
    }
  }
  delete[] vv;
  return swap;
}

namespace run {
void dividebyzero(size_t i)
{
  ostringstream buf;
  if(i > 0) buf << "array element " << i << ": ";
  buf << "Divide by zero";
  error(buf);
}
  
void integeroverflow(size_t i)
{
  ostringstream buf;
  if(i > 0) buf << "array element " << i << ": ";
  buf << "Integer overflow";
  error(buf);
}
}

// Autogenerated routines:


// Create an empty array.
array* :emptyArray()
{
  return new array(0);
}

// 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.
array* :newDeepArray(Int depth)
{
  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;
  return a;
}

// Creates an array with elements already specified.  First, the number
// of elements is popped off the stack, followed by each element in
// reverse order.
array* :newInitializedArray(Int n)
{
  assert(n >= 0);

  array *a = new array(n);

  for (Int index = n-1; index >= 0; index--)
    (*a)[index] = pop(Stack);

  return a;
}

// Similar to newInitializedArray, but after the n elements, append another
// array to it.
array* :newAppendedArray(array* tail, Int n)
{
  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));

  return a;
}

// The function T[] array(int n, T value, int depth=0) produces a array of n
// copies of x, where each copy is copied up to depth.
array* :newDuplicateArray(Int n, item value, Int depth=Int_MAX)
{
  if(n < 0) error("cannot create a negative length array");
  if(depth < 0) error("cannot copy to a negative depth");

  return new array(n, value, depth);
}

// Read an element from an array. Checks for initialization & bounds.
item :arrayRead(array *a, Int n)
{
  item& i=arrayRead(a,n);
  if (i.empty()) {
    ostringstream buf;
    buf << "read uninitialized value from array at index " << n;
    error(buf);
  }
  return i;
}

// Slice a substring from an array.
item :arraySliceRead(array *a, Int left, Int right)
{
  checkArray(a);
  return a->slice(left, right);
}

// 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.
item :arraySliceReadToEnd(array *a, Int left)
{
  size_t len=checkArray(a);
  return a->slice(left, (Int)len);
}

// Read an element from an array of arrays. Check bounds and initialize
// as necessary.
item :arrayArrayRead(array *a, Int n)
{
  item& i=arrayRead(a,n);
  if (i.empty()) i=new array(0);
  return i;
}

// Write an element to an array.  Increase size if necessary.
item :arrayWrite(item value, array *a, Int n)
{
  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;
  return value;
}

array * :arraySliceWrite(array *src, array *dest, Int left, Int right)
{
  checkArray(src);
  checkArray(dest);
  dest->setSlice(left, right, src);
  return src;
}

array * :arraySliceWriteToEnd(array *src, array *dest, Int left)
{
  checkArray(src);
  size_t len=checkArray(dest);
  dest->setSlice(left, (Int) len, src);
  return src;
}

// Returns the length of an array.
Int :arrayLength(array *a)
{
  return (Int) checkArray(a);
}

// Returns an array of integers representing the keys of the array.
array * :arrayKeys(array *a)
{
  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);
  }

  return keys;
}

// Return the cyclic flag for an array.
bool :arrayCyclicFlag(array *a)
{
  checkArray(a);
  return a->cyclic();
}

bool :arraySetCyclicFlag(bool b, array *a)
{
  checkArray(a);
  a->cyclic(b);
  return b;
}

// Check to see if an array element is initialized.
bool :arrayInitializedHelper(Int n, array *a)
{
  size_t len=checkArray(a);
  bool cyclic=a->cyclic();
  if(cyclic && len > 0) n=imod(n,len);
  else if(n < 0 || n >= (Int) len) return false;
  item&i=(*a)[(unsigned) n];
  return !i.empty();
}

// Returns the initialize method for an array.
callable* :arrayInitialized(array *a)
{
  return new thunk(new bfunc(arrayInitializedHelper),a);
}

// The helper function for the cyclic method that sets the cyclic flag.
void :arrayCyclicHelper(bool b, array *a)
{
  checkArray(a);
  a->cyclic(b);
}

// Set the cyclic flag for an array.
callable* :arrayCyclic(array *a)
{
  return new thunk(new bfunc(arrayCyclicHelper),a);
}

// The helper function for the push method that does the actual operation.
item :arrayPushHelper(item x, array *a)
{
  checkArray(a);
  a->push(x);
  return x;
}

// Returns the push method for an array.
callable* :arrayPush(array *a)
{
  return new thunk(new bfunc(arrayPushHelper),a);
}

// The helper function for the append method that appends b to a.
void :arrayAppendHelper(array *b, array *a)
{
  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.
callable* :arrayAppend(array *a)
{
  return new thunk(new bfunc(arrayAppendHelper),a);
}

// The helper function for the pop method.
item :arrayPopHelper(array *a)
{
  size_t asize=checkArray(a);
  if(asize == 0) 
    error("cannot pop element from empty array");
  return a->pop();
}

// Returns the pop method for an array.
callable* :arrayPop(array *a)
{
  return new thunk(new bfunc(arrayPopHelper),a);
}

// The helper function for the insert method.
item :arrayInsertHelper(Int i, array *x, array *a)
{
  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.
callable* :arrayInsert(array *a)
{
  return new thunk(new bfunc(arrayInsertHelper),a);
}

// Returns the delete method for an array.
callable* :arrayDelete(array *a)
{
  return new thunk(new bfunc(arrayDeleteHelper),a);
}

bool :arrayAlias(array *a, array *b)
{
  return a==b;
}

// Return array formed by indexing array a with elements of integer array b
array* :arrayIntArray(array *a, array *b)
{
  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];
  }
  return r;
}

// 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].
Intarray* complement(Intarray *a, Int n)
{
  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;
  return r;
}

// Generate the sequence {f(i) : i=0,1,...n-1} given a function f and integer n
Intarray* :arraySequence(callable *f, Int n)
{
  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);
  }
  return a;
}

// Return the array {0,1,...n-1}
Intarray *sequence(Int n)
{
  if(n < 0) n=0;
  array *a=new array(n);
  for(Int i=0; i < n; ++i) {
    (*a)[i]=i;
  }
  return a;
}

// Apply a function to each element of an array
array* :arrayFunction(callable *f, array *a)
{
  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);
  }
  return b;
}

array* :arraySort(array *a, callable *f)
{
  array *c=copyArray(a);
  compareFunc=f;
  FuncStack=Stack;
  stable_sort(c->begin(),c->end(),compareFunction);
  return c;
}

bool all(boolarray *a)
{
  size_t size=checkArray(a);
  bool c=true;
  for(size_t i=0; i < size; i++)
    if(!get<bool>((*a)[i])) {c=false; break;}
  return c;
}

boolarray* !(boolarray* a)
{
  size_t size=checkArray(a);
  array *c=new array(size);
  for(size_t i=0; i < size; i++)
    (*c)[i]=!read<bool>(a,i);
  return c;
}

Int sum(boolarray *a)
{
  size_t size=checkArray(a);
  Int sum=0;
  for(size_t i=0; i < size; i++)
    sum += read<bool>(a,i) ? 1 : 0;
  return sum;
}

array* :arrayCopy(array *a)
{
  return copyArray(a);
}

array* :arrayConcat(array *a)
{
  // 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;
    }
  }

  return result;
}

array* :array2Copy(array *a)
{
  return copyArray2(a);
}

array* :array3Copy(array *a)
{
  return copyArray3(a);
}

array* :array2Transpose(array *a)
{
  size_t asize=checkArray(a);
  array *c=new array(0);
  for(size_t i=0; i < asize; i++) {
    size_t ip=i+1;
    array *ai=read<array*>(a,i);
    size_t aisize=checkArray(ai);
    size_t csize=checkArray(c);
    if(csize < aisize) {
      c->resize(aisize);
      for(size_t j=csize; j < aisize; j++) {
        (*c)[j]=new array(ip);
      }
    }
    for(size_t j=0; j < aisize; j++) {
      array *cj=read<array*>(c,j);
      if(checkArray(cj) < ip) cj->resize(ip);
      (*cj)[i]=(*ai)[j];
    }
  }
  return c;
}

// 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}. 
array* :array3Transpose(array *a, array *perm)
{
  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;  

  return c;
}

// In a boolean array, find the index of the nth true value or -1 if not found
// If n is negative, search backwards.
Int find(boolarray *a, Int n=1)
{
  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;}
      }
  return j;
}

// construct vector obtained by replacing those elements of b for which the
// corresponding elements of a are false by the corresponding element of c.
array* :arrayConditional(array *a, array *b, array *c)
{
  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]);
    }
  }
  return r;
}

// Return an n x n identity matrix.
realarray2 *identity(Int n)
{
  return Identity(n);
}

// Return the diagonal matrix with diagonal entries given by a.
realarray2* :diagonal(realarray *a)
{
  size_t n=checkArray(a);
  array *c=new array(n);
  for(size_t i=0; i < n; ++i) {
    array *ci=new array(n);
    (*c)[i]=ci;
    for(size_t j=0; j < i; ++j)
      (*ci)[j]=0.0;
    (*ci)[i]=read<real>(a,i);
    for(size_t j=i+1; j < n; ++j)
      (*ci)[j]=0.0;
  }
  return c;
}

// Return the inverse of an n x n matrix a using Gauss-Jordan elimination.
realarray2 *inverse(realarray2 *a)
{
  a=copyArray2(a);
  size_t n=checkArray(a);
  checkSquare(a);
  
  inverseAllocate(n);
  
  for(size_t i=0; i < n; i++)
    pivot[i]=0;
 
  size_t col=0, row=0;
  // This is the main loop over the columns to be reduced.
  for(size_t i=0; i < n; i++) {
    real big=0.0;
    // This is the outer loop of the search for a pivot element.
    for(size_t j=0; j < n; j++) {
      array *aj=read<array*>(a,j);
      if(pivot[j] != 1) {
        for(size_t k=0; k < n; k++) {
          if(pivot[k] == 0) {
            real temp=fabs(read<real>(aj,k));
            if(temp >= big) {
              big=temp;
              row=j;
              col=k;
            }
          } else if(pivot[k] > 1) {
            inverseDeallocate();
            error(singular);
          }
        }
      }
    }
    ++(pivot[col]);
    
    // Interchange rows, if needed, to put the pivot element on the diagonal.
    array *acol=read<array*>(a,col);
    if(row != col) {
      array *arow=read<array*>(a,row);
      for(size_t l=0; l < n; l++) {
        real temp=read<real>(arow,l);
        (*arow)[l]=read<real>(acol,l);
        (*acol)[l]=temp;
      }
    }
    
    Row[i]=row; 
    Col[i]=col;

    // Divide the pivot row by the pivot element.
    real denom=read<real>(acol,col);
    if(denom == 0.0) {
      inverseDeallocate();
      error(singular);
    }
    real pivinv=1.0/denom;
    (*acol)[col]=1.0;
    for(size_t l=0; l < n; l++) 
      (*acol)[l]=read<real>(acol,l)*pivinv;
    
    // Reduce all rows except for the pivoted one.
    for(size_t k=0; k < n; k++) {
      if(k != col) { 
        array *ak=read<array*>(a,k);
        real akcol=read<real>(ak,col);
        (*ak)[col]=0.0;
        for(size_t l=0; l < n; l++)
          (*ak)[l]=read<real>(ak,l)-read<real>(acol,l)*akcol;
      }
    }
  }
  
  // Unscramble the inverse matrix in view of the column interchanges.
  for(size_t l=n; l > 0;) {
    l--;
    size_t r=Row[l];
    size_t c=Col[l];
    if(r != c) {
      for(size_t k=0; k < n; k++) {
        array *ak=read<array*>(a,k);
        real temp=read<real>(ak,r);
        (*ak)[r]=read<real>(ak,c);
        (*ak)[c]=temp;
      }
    }
  }
  inverseDeallocate();
  return a;
}

// 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.
realarray *solve(realarray2 *a, realarray *b, bool warn=true)
{
  size_t n=checkArray(a);
  
  if(n == 0) return new array(0);
  
  size_t m=checkArray(b);
  if(m != n) error(incommensurate);
  
  real *A=copyArray2C(a);
  size_t *index=new size_t[n];
  
  if(LUdecompose(A,n,index,warn) == 0)
    return new array(0);

  array *x=new array(n);
  
  real *B=copyArrayC(b);
  
  for(size_t i=0; i < n; ++i) {
    size_t ip=index[i];
    real sum=B[ip];
    B[ip]=B[i];
    real *Ai=A+i*n;
    for(size_t j=0; j < i; ++j)
      sum -= Ai[j]*B[j];
    B[i]=sum;
  }
  
  for(size_t i=n; i > 0;) {
    --i;
    real sum=B[i];
    real *Ai=A+i*n;
    for(size_t j=i+1; j < n; ++j)
      sum -= Ai[j]*B[j];
    B[i]=sum/Ai[i];
  }
  
  for(size_t i=0; i < n; ++i)
    (*x)[i]=B[i];

  delete[] index;
  delete[] B;
  delete[] A;
  
  return x;
}

// 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.
realarray2 *solve(realarray2 *a, realarray2 *b, bool warn=true)
{
  size_t n=checkArray(a);
  
  if(n == 0) return new array(0);
  
  if(checkArray(b) != n) error(incommensurate);
  size_t m=checkArray(read<array*>(b,0));
  
  real *A=copyArray2C(a);
  real *B=copyArray2C(b,false);
  
  size_t *index=new size_t[n];
  
  if(LUdecompose(A,n,index,warn) == 0)
    return new array(0);

  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;
  
  return x;
}

// Compute the determinant of an n x n matrix.
real determinant(realarray2 *a)
{
  real *A=copyArray2C(a);
  size_t n=checkArray(a);
  
  real det=LUdecompose(A,n,NULL,false);
  size_t n1=n+1;
  for(size_t i=0; i < n; ++i)
    det *= A[i*n1];
  
  delete[] A;
  
  return det;
}

realarray *Operator *(realarray2 *a, realarray *b)
{
  size_t n=checkArray(a);
  size_t m=checkArray(b);
  array *c=new array(n);
  real *B=copyArrayC(b);
  for(size_t i=0; i < n; ++i) {
    array *ai=read<array*>(a,i);
    if(checkArray(ai) != m) error(incommensurate);
    real sum=0.0;
    for(size_t j=0; j < m; ++j)
      sum += read<real>(ai,j)*B[j];
    (*c)[i]=sum;
  }
  delete[] B;
  return c;
}

realarray *Operator *(realarray *a, realarray2 *b)
{
  size_t n=checkArray(a);
  if(n != checkArray(b)) error(incommensurate);
  real *A=copyArrayC(a);

  array **B=new array*[n];
  array *bk=read<array *>(b,0);
  B[0]=bk;
  size_t m=bk->size();
  for(size_t k=1; k < n; k++) {
    array *bk=read<array *>(b,k);
    if(bk->size() != m) error(incommensurate);
    B[k]=bk;
  }
  array *c=new array(m);

  for(size_t i=0; i < m; ++i) {
    real sum=0.0;
    for(size_t k=0; k < n; ++k)
      sum += A[k]*read<real>(B[k],i);
    (*c)[i]=sum;
  }
  delete[] B;
  delete[] A;
  return c;
}

realarray2 *Operator *(realarray2 *a, realarray2 *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);

  real *A=copyArray2C(a,false);
  real *B=copyArray2C(b,false);

  for(size_t i=0; i < n; ++i) {
    real *Ai=A+i*nb;
    array *ci=new array(nb0);
    (*c)[i]=ci;
    for(size_t j=0; j < nb0; ++j) {
      real sum=0.0;
      size_t kj=j;
      for(size_t k=0; k < nb; ++k, kj += nb0)
        sum += Ai[k]*B[kj];
      (*ci)[j]=sum;
    }
  }
  
  delete[] B;
  delete[] A;
  
  return c;
}

triple Operator *(realarray2 *t, triple v)
{
  return *t*v;
}

pair project(triple v, realarray2 *t)
{
  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;
  
  return 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);
}

// Compute the dot product of vectors a and b.
real dot(realarray *a, realarray *b) 
{
  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);
  return sum;
}

// 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] ]
realarray *tridiagonal(realarray *a, realarray *b, realarray *c, realarray *f)
{
  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) return up;
  
  // 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;
    return up;
  }
  
  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; return up;}
  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;
    return up;
  }
        
  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;
  
  return up;
}

// Root solve by Newton-Raphson
real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x,
            bool verbose=false)
{
  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);
  return x;
}

// Root solve by Newton-Raphson bisection
// cf. routine rtsafe (Press et al.,  Numerical Recipes, 1991).
real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x1,
            real x2, bool verbose=false)
{
  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) return x1;
  
  Stack->push(x2);
  f->call(Stack);
  real f2=pop<real>(Stack);
  if(f2 == 0.0) return x2;
        
  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) return x;
    } else {
      dxold=dx;
      dx=y/dy;
      real temp=x;
      x -= dx;
      if(verbose)
        cout << "Newton-Raphson: " << x << endl;
      if(temp == x) return x;
    }
    if(fabs(dx) < fuzz*fabs(x)) return x;
    
    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);
  return (j == iterations) ? DBL_MAX : x;
}

real simpson(callableReal *f, real a, real b, real acc=DBL_EPSILON,
             real dxmax=0)
{
  real integral;
  if(dxmax == 0) dxmax=b-a;
  Func=f;
  FuncStack=Stack;
  if(!simpson(integral,wrapFunction,a,b,acc,dxmax))
    error("nesting capacity exceeded in simpson");
  return integral;
}

// Compute the fast Fourier transform of a pair array
pairarray* fft(pairarray *a, Int sign=1)
{
  unsigned n=(unsigned) checkArray(a);
#ifdef HAVE_LIBFFTW3
  array *c=new array(n);
  if(n) {
    Complex *f=FFTWComplex(n);
    fft1d Forward(n,intcast(sign),f);
  
    for(size_t i=0; i < n; i++) {
      pair z=read<pair>(a,i);
      f[i]=Complex(z.getx(),z.gety());
    }
    Forward.fft(f);
  
    for(size_t i=0; i < n; i++) {
      Complex z=f[i];
      (*c)[i]=pair(z.real(),z.imag());
    }
    FFTWdelete(f);
  }
#else
  unused(&n);
  unused(&sign);
  array *c=new array(0);
  error("Please install fftw3, run ./configure, and recompile");
#endif //  HAVE_LIBFFTW3
  return c;
}

Intarray2 *triangulate(pairarray *z)
{
  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;
  return t;
}

real norm(realarray *a)
{
  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;
  }
  return M;
}

real norm(realarray2 *a)
{
  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;
    }
  }
  return M;
}

real norm(triplearray2 *a)
{
  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;
    }
  }
  return sqrt(M);
}

real change2(triplearray2 *a)
{
  size_t n=checkArray(a);
  if(n == 0) return 0.0;
  
  vm::array *a0=vm::read<vm::array*>(a,0);
  size_t m=checkArray(a0);
  if(m == 0) return 0.0;
  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;
    }
  }
  return M;
}

triple minbezier(triplearray2 *P, triple b)
{
  real *A=copyTripleArray2Components(P,true,4);
  b=triple(bound(A,::min,b.getx(),sqrtFuzz*norm(A,16)),
           bound(A+16,::min,b.gety(),sqrtFuzz*norm(A+16,16)),
           bound(A+32,::min,b.getz(),sqrtFuzz*norm(A+32,16)));
  delete[] A;
  return b;
}

triple maxbezier(triplearray2 *P, triple b)
{
  real *A=copyTripleArray2Components(P,true,4);
  b=triple(bound(A,::max,b.getx(),sqrtFuzz*norm(A,16)),
           bound(A+16,::max,b.gety(),sqrtFuzz*norm(A+16,16)),
           bound(A+32,::max,b.getz(),sqrtFuzz*norm(A+32,16)));
  delete[] A;
  return b;
}

pair minratio(triplearray2 *P, pair b)
{
  triple *A=copyTripleArray2C(P,true,4);
  real fuzz=sqrtFuzz*norm(A,16);
  b=pair(bound(A,::min,xratio,b.getx(),fuzz),
         bound(A,::min,yratio,b.gety(),fuzz));
  delete[] A;
  return b;
}

pair maxratio(triplearray2 *P, pair b)
{
  triple *A=copyTripleArray2C(P,true,4);
  real fuzz=sqrtFuzz*norm(A,16);
  b=pair(bound(A,::max,xratio,b.getx(),fuzz),
         bound(A,::max,yratio,b.gety(),fuzz));
  delete[] A;
  return b;
}

realarray *_projection()
{
#ifdef HAVE_LIBGL
  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
  return new array(0);
}