From ad547a6b5986815fda458221149728d9d9ab1d87 Mon Sep 17 00:00:00 2001 From: Karl Berry Date: Thu, 25 Feb 2021 19:22:25 +0000 Subject: restore Build,TODO from r57911 git-svn-id: svn://tug.org/texlive/trunk@57915 c570f23f-e606-0410-a88d-b1316a301751 --- Build/source/utils/asymptote/runarray.cc | 2715 ++++++++++++++++++++++++++++++ 1 file changed, 2715 insertions(+) create mode 100644 Build/source/utils/asymptote/runarray.cc (limited to 'Build/source/utils/asymptote/runarray.cc') diff --git a/Build/source/utils/asymptote/runarray.cc b/Build/source/utils/asymptote/runarray.cc new file mode 100644 index 00000000000..05c7f0a1a7c --- /dev/null +++ b/Build/source/utils/asymptote/runarray.cc @@ -0,0 +1,2715 @@ +/***** Autogenerated from runarray.in; changes will be overwritten *****/ + +#line 1 "runtimebase.in" +/***** + * runtimebase.in + * Andy Hammerlindl 2009/07/28 + * + * Common declarations needed for all code-generating .in files. + * + *****/ + + +#line 1 "runarray.in" +/***** + * runarray.in + * + * Runtime functions for array operations. + * + *****/ + +#line 1 "runtimebase.in" +#include "stack.h" +#include "types.h" +#include "builtin.h" +#include "entry.h" +#include "errormsg.h" +#include "array.h" +#include "triple.h" +#include "callable.h" +#include "opsymbols.h" + +using vm::stack; +using vm::error; +using vm::array; +using vm::read; +using vm::callable; +using types::formal; +using types::function; +using camp::triple; + +#define PRIMITIVE(name,Name,asyName) using types::prim##Name; +#include +#undef PRIMITIVE + +typedef double real; + +void unused(void *); + +namespace run { +array *copyArray(array *a); +array *copyArray2(array *a); +array *copyArray3(array *a); + +double *copyTripleArray2Components(array *a, size_t &N, + GCPlacement placement=NoGC); +triple *copyTripleArray2C(array *a, size_t &N, + GCPlacement placement=NoGC); +} + +function *realRealFunction(); + +#define CURRENTPEN processData().currentpen + +#line 22 "runarray.in" +#include "array.h" +#include "arrayop.h" +#include "triple.h" +#include "path3.h" +#include "Delaunay.h" +#include "glrender.h" + +#ifdef HAVE_LIBFFTW3 +#include "fftw++.h" + static const char *rectangular="matrix must be rectangular"; +#else +static const char *installFFTW= + "Please install fftw3, run ./configure, and recompile"; +#endif + +using namespace camp; +using namespace vm; + +namespace run { +extern pair zero; +} + +typedef array boolarray; +typedef array Intarray; +typedef array Intarray2; +typedef array realarray; +typedef array realarray2; +typedef array pairarray; +typedef array pairarray2; +typedef array pairarray3; +typedef array triplearray2; + +using types::booleanArray; +using types::IntArray; +using types::IntArray2; +using types::realArray; +using types::realArray2; +using types::pairArray; +using types::pairArray2; +using types::pairArray3; +using types::tripleArray2; + +typedef callable callableReal; + +void outOfBounds(const char *op, size_t len, Int n) +{ + ostringstream buf; + buf << op << " array of length " << len << " with out-of-bounds index " << n; + error(buf); +} + +inline item& arrayRead(array *a, Int n) +{ + size_t len=checkArray(a); + bool cyclic=a->cyclic(); + if(cyclic && len > 0) n=imod(n,len); + else if(n < 0 || n >= (Int) len) outOfBounds("reading",len,n); + return (*a)[(unsigned) n]; +} + +// Helper function to create deep arrays. +static array* deepArray(Int depth, Int *dims) +{ + assert(depth > 0); + + if (depth == 1) { + return new array(dims[0]); + } else { + Int length = dims[0]; + depth--; dims++; + + array *a = new array(length); + + for (Int index = 0; index < length; index++) { + (*a)[index] = deepArray(depth, dims); + } + return a; + } +} + +namespace run { +array *Identity(Int n) +{ + size_t N=(size_t) n; + array *c=new array(N); + for(size_t i=0; i < N; ++i) { + array *ci=new array(N); + (*c)[i]=ci; + for(size_t j=0; j < N; ++j) + (*ci)[j]=0.0; + (*ci)[i]=1.0; + } + return c; +} +} + +static const char *incommensurate="Incommensurate matrices"; +static const char *singular="Singular matrix"; +static const char *invalidarraylength="Invalid array length: "; + +static size_t *pivot,*Row,*Col; + +bound_double *bounddouble(int N) +{ + if(N == 16) return bound; + if(N == 10) return boundtri; + ostringstream buf; + buf << invalidarraylength << " " << N; + error(buf); + return NULL; +} + +bound_triple *boundtriple(int N) +{ + if(N == 16) return bound; + if(N == 10) return boundtri; + ostringstream buf; + buf << invalidarraylength << " " << N; + error(buf); + return NULL; +} + +static inline void inverseAllocate(size_t n) +{ + pivot=new size_t[n]; + Row=new size_t[n]; + Col=new size_t[n]; +} + +static inline void inverseDeallocate() +{ + delete[] pivot; + delete[] Row; + delete[] Col; +} + +namespace run { + +array *copyArray(array *a) +{ + size_t size=checkArray(a); + array *c=new array(size); + for(size_t i=0; i < size; i++) + (*c)[i]=(*a)[i]; + return c; +} + +array *copyArray2(array *a) +{ + size_t size=checkArray(a); + array *c=new array(size); + for(size_t i=0; i < size; i++) { + array *ai=read(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(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(a,i); + size_t m=checkArray(ai); + for(size_t j=0; j < m; j++) { + triple v=read(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(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(a,i); + size_t m=checkArray(ai); + for(size_t j=0; j < m; j++) + *(p++)=read(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(t,0); + array *t1=read(t,1); + array *t2=read(t,2); + array *t3=read(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(t3,0)*x+read(t3,1)*y+read(t3,2)*z+ + read(t3,3); + if(f == 0.0) run::dividebyzero(); + f=1.0/f; + + return triple((read(t0,0)*x+read(t0,1)*y+read(t0,2)*z+ + read(t0,3))*f, + (read(t1,0)*x+read(t1,1)*y+read(t1,2)*z+ + read(t1,3))*f, + (read(t2,0)*x+read(t2,1)*y+read(t2,2)*z+ + read(t2,3))*f); +} + +template +array *mult(array *a, array *b) +{ + size_t n=checkArray(a); + + size_t nb=checkArray(b); + size_t na0=n == 0 ? 0 : checkArray(read(a,0)); + if(na0 != nb) + error(incommensurate); + + size_t nb0=nb == 0 ? 0 : checkArray(read(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 +array *AtA(array *a) +{ + size_t n=checkArray(a); + size_t m=n == 0 ? 0 : checkArray(read(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(FuncStack); +} + +callable *compareFunc; +bool compareFunction(const vm::item& i, const vm::item& j) +{ + FuncStack->push(i); + FuncStack->push(j); + compareFunc->call(FuncStack); + return pop(FuncStack); +} + +// Crout's algorithm for computing the LU decomposition of a square matrix. +// cf. routine ludcmp (Press et al., Numerical Recipes, 1991). +Int LUdecompose(double *a, size_t n, size_t* index, bool warn=true) +{ + double *vv=new double[n]; + Int swap=1; + for(size_t i=0; i < n; ++i) { + double big=0.0; + double *ai=a+i*n; + for(size_t j=0; j < n; ++j) { + double temp=fabs(ai[j]); + if(temp > big) big=temp; + } + if(big == 0.0) { + delete[] vv; + if(warn) error(singular); + else return 0; + } + vv[i]=1.0/big; + } + for(size_t j=0; j < n; ++j) { + for(size_t i=0; i < j; ++i) { + double *ai=a+i*n; + double sum=ai[j]; + for(size_t k=0; k < i; ++k) { + sum -= ai[k]*a[k*n+j]; + } + ai[j]=sum; + } + double big=0.0; + size_t imax=j; + for(size_t i=j; i < n; ++i) { + double *ai=a+i*n; + double sum=ai[j]; + for(size_t k=0; k < j; ++k) + sum -= ai[k]*a[k*n+j]; + ai[j]=sum; + double temp=vv[i]*fabs(sum); + if(temp >= big) { + big=temp; + imax=i; + } + } + double *aj=a+j*n; + double *aimax=a+imax*n; + if(j != imax) { + for(size_t k=0; k < n; ++k) { + double temp=aimax[k]; + aimax[k]=aj[k]; + aj[k]=temp; + } + swap *= -1; + vv[imax]=vv[j]; + } + if(index) + index[j]=imax; + if(j != n) { + double denom=aj[j]; + if(denom == 0.0) { + delete[] vv; + if(warn) error(singular); + else return 0; + } + for(size_t i=j+1; i < n; ++i) + a[i*n+j] /= denom; + } + } + delete[] vv; + return swap; +} + +namespace run { + +void dividebyzero(size_t i) +{ + ostringstream buf; + if(i > 0) buf << "array element " << i << ": "; + buf << "Divide by zero"; + error(buf); +} + +void integeroverflow(size_t i) +{ + ostringstream buf; + if(i > 0) buf << "array element " << i << ": "; + buf << "Integer overflow"; + error(buf); +} + +} + +// Autogenerated routines: + + + +#ifndef NOSYM +#include "runarray.symbols.h" + +#endif +namespace run { +// Create an empty array. +#line 594 "runarray.in" +void emptyArray(stack *Stack) +{ +#line 595 "runarray.in" + {Stack->push(new array(0)); return;} +} + +// Create a new array (technically a vector). +// This array will be multidimensional. First the number of dimensions +// is popped off the stack, followed by each dimension in reverse order. +// The array itself is technically a one dimensional array of one +// dimension arrays and so on. +#line 604 "runarray.in" +void newDeepArray(stack *Stack) +{ + Int depth=vm::pop(Stack); +#line 605 "runarray.in" + assert(depth > 0); + + Int *dims = new Int[depth]; + + for (Int index = depth-1; index >= 0; index--) { + Int i=pop(Stack); + if(i < 0) error("cannot create a negative length array"); + dims[index]=i; + } + + array *a=deepArray(depth, dims); + delete[] dims; + {Stack->push(a); return;} +} + +// Creates an array with elements already specified. First, the number +// of elements is popped off the stack, followed by each element in +// reverse order. +#line 624 "runarray.in" +void newInitializedArray(stack *Stack) +{ + Int n=vm::pop(Stack); +#line 625 "runarray.in" + assert(n >= 0); + + array *a = new array(n); + + for (Int index = n-1; index >= 0; index--) + (*a)[index] = pop(Stack); + + {Stack->push(a); return;} +} + +// Similar to newInitializedArray, but after the n elements, append another +// array to it. +#line 638 "runarray.in" +void newAppendedArray(stack *Stack) +{ + Int n=vm::pop(Stack); + array* tail=vm::pop(Stack); +#line 639 "runarray.in" + assert(n >= 0); + + array *a = new array(n); + + for (Int index = n-1; index >= 0; index--) + (*a)[index] = pop(Stack); + + copy(tail->begin(), tail->end(), back_inserter(*a)); + + {Stack->push(a); return;} +} + +// Produce an array of n deep copies of value. +// typeDepth is the true depth of the array determined at compile-time when the +// operations for the array type are added. This typeDepth argument is +// automatically pushed on the stack and is not visible to the user. +#line 656 "runarray.in" +void copyArrayValue(stack *Stack) +{ + Int typeDepth=vm::pop(Stack); + Int depth=vm::pop(Stack,Int_MAX); + item value=vm::pop(Stack); + Int n=vm::pop(Stack); +#line 657 "runarray.in" + if(n < 0) error("cannot create a negative length array"); + if(depth < 0) error("cannot copy to a negative depth"); + if(depth > typeDepth) depth=typeDepth; + {Stack->push(new array((size_t) n, value, depth)); return;} +} + +// Deep copy of array. +// typeDepth is the true depth of the array determined at compile-time when the +// operations for the array type are added. This typeDepth argument is +// automatically pushed on the stack and is not visible to the user. +#line 668 "runarray.in" +void copyArray(stack *Stack) +{ + Int typeDepth=vm::pop(Stack); + Int depth=vm::pop(Stack,Int_MAX); + array * a=vm::pop(Stack); +#line 669 "runarray.in" + if(a == 0) vm::error(dereferenceNullArray); + if(depth < 0) error("cannot copy to a negative depth"); + if(depth > typeDepth) depth=typeDepth; + {Stack->push(a->copyToDepth(depth)); return;} +} + +// Read an element from an array. Checks for initialization & bounds. +#line 677 "runarray.in" +void arrayRead(stack *Stack) +{ + Int n=vm::pop(Stack); + array * a=vm::pop(Stack); +#line 678 "runarray.in" + item& i=arrayRead(a,n); + if (i.empty()) { + ostringstream buf; + buf << "read uninitialized value from array at index " << n; + error(buf); + } + {Stack->push(i); return;} +} + +// Slice a substring from an array. +#line 689 "runarray.in" +void arraySliceRead(stack *Stack) +{ + Int right=vm::pop(Stack); + Int left=vm::pop(Stack); + array * a=vm::pop(Stack); +#line 690 "runarray.in" + checkArray(a); + {Stack->push(a->slice(left, right)); return;} +} + +// Slice a substring from an array. This implements the cases a[i:] and a[:] +// where the endpoint is not given, and assumed to be the length of the array. +#line 697 "runarray.in" +void arraySliceReadToEnd(stack *Stack) +{ + Int left=vm::pop(Stack); + array * a=vm::pop(Stack); +#line 698 "runarray.in" + size_t len=checkArray(a); + {Stack->push(a->slice(left, (Int)len)); return;} +} + +// Read an element from an array of arrays. Check bounds and initialize +// as necessary. +#line 705 "runarray.in" +void arrayArrayRead(stack *Stack) +{ + Int n=vm::pop(Stack); + array * a=vm::pop(Stack); +#line 706 "runarray.in" + item& i=arrayRead(a,n); + if (i.empty()) i=new array(0); + {Stack->push(i); return;} +} + +// Write an element to an array. Increase size if necessary. +// TODO: Add arrayWriteAndPop +#line 714 "runarray.in" +void arrayWrite(stack *Stack) +{ + item value=vm::pop(Stack); + Int n=vm::pop(Stack); + array * a=vm::pop(Stack); +#line 715 "runarray.in" + size_t len=checkArray(a); + bool cyclic=a->cyclic(); + if(cyclic && len > 0) n=imod(n,len); + else { + if(cyclic) outOfBounds("writing cyclic",len,n); + if(n < 0) outOfBounds("writing",len,n); + if(len <= (size_t) n) + a->resize(n+1); + } + (*a)[n] = value; + {Stack->push(value); return;} +} + +#line 729 "runarray.in" +void arraySliceWrite(stack *Stack) +{ + array * src=vm::pop(Stack); + Int right=vm::pop(Stack); + Int left=vm::pop(Stack); + array * dest=vm::pop(Stack); +#line 730 "runarray.in" + checkArray(src); + checkArray(dest); + dest->setSlice(left, right, src); + {Stack->push(src); return;} +} + +#line 737 "runarray.in" +void arraySliceWriteToEnd(stack *Stack) +{ + array * src=vm::pop(Stack); + Int left=vm::pop(Stack); + array * dest=vm::pop(Stack); +#line 738 "runarray.in" + checkArray(src); + size_t len=checkArray(dest); + dest->setSlice(left, (Int) len, src); + {Stack->push(src); return;} +} + +// Returns the length of an array. +#line 746 "runarray.in" +void arrayLength(stack *Stack) +{ + array * a=vm::pop(Stack); +#line 747 "runarray.in" + {Stack->push((Int) checkArray(a)); return;} +} + +// Returns an array of integers representing the keys of the array. +#line 752 "runarray.in" +void arrayKeys(stack *Stack) +{ + array * a=vm::pop(Stack); +#line 753 "runarray.in" + size_t size=checkArray(a); + + array *keys=new array(); + for (size_t i=0; ipush((Int)i); + } + + {Stack->push(keys); return;} +} + +// Return the cyclic flag for an array. +#line 767 "runarray.in" +void arrayCyclicFlag(stack *Stack) +{ + array * a=vm::pop(Stack); +#line 768 "runarray.in" + checkArray(a); + {Stack->push(a->cyclic()); return;} +} + +#line 773 "runarray.in" +void arraySetCyclicFlag(stack *Stack) +{ + array * a=vm::pop(Stack); + bool b=vm::pop(Stack); +#line 774 "runarray.in" + checkArray(a); + a->cyclic(b); + {Stack->push(b); return;} +} + +// Check to see if an array element is initialized. +#line 781 "runarray.in" +void arrayInitializedHelper(stack *Stack) +{ + array * a=vm::pop(Stack); + Int n=vm::pop(Stack); +#line 782 "runarray.in" + size_t len=checkArray(a); + bool cyclic=a->cyclic(); + if(cyclic && len > 0) n=imod(n,len); + else if(n < 0 || n >= (Int) len) {Stack->push(false); return;} + item&i=(*a)[(unsigned) n]; + {Stack->push(!i.empty()); return;} +} + +// Returns the initialize method for an array. +#line 792 "runarray.in" +void arrayInitialized(stack *Stack) +{ + array * a=vm::pop(Stack); +#line 793 "runarray.in" + {Stack->push(new thunk(new bfunc(arrayInitializedHelper),a)); return;} +} + +// The helper function for the cyclic method that sets the cyclic flag. +#line 798 "runarray.in" +void arrayCyclicHelper(stack *Stack) +{ + array * a=vm::pop(Stack); + bool b=vm::pop(Stack); +#line 799 "runarray.in" + checkArray(a); + a->cyclic(b); +} + +// Set the cyclic flag for an array. +#line 805 "runarray.in" +void arrayCyclic(stack *Stack) +{ + array * a=vm::pop(Stack); +#line 806 "runarray.in" + {Stack->push(new thunk(new bfunc(arrayCyclicHelper),a)); return;} +} + +// The helper function for the push method that does the actual operation. +#line 811 "runarray.in" +void arrayPushHelper(stack *Stack) +{ + array * a=vm::pop(Stack); + item x=vm::pop(Stack); +#line 812 "runarray.in" + checkArray(a); + a->push(x); + {Stack->push(x); return;} +} + +// Returns the push method for an array. +#line 819 "runarray.in" +void arrayPush(stack *Stack) +{ + array * a=vm::pop(Stack); +#line 820 "runarray.in" + {Stack->push(new thunk(new bfunc(arrayPushHelper),a)); return;} +} + +// The helper function for the append method that appends b to a. +#line 825 "runarray.in" +void arrayAppendHelper(stack *Stack) +{ + array * a=vm::pop(Stack); + array * b=vm::pop(Stack); +#line 826 "runarray.in" + checkArray(a); + size_t size=checkArray(b); + for(size_t i=0; i < size; i++) + a->push((*b)[i]); +} + +// Returns the append method for an array. +#line 834 "runarray.in" +void arrayAppend(stack *Stack) +{ + array * a=vm::pop(Stack); +#line 835 "runarray.in" + {Stack->push(new thunk(new bfunc(arrayAppendHelper),a)); return;} +} + +// The helper function for the pop method. +#line 840 "runarray.in" +void arrayPopHelper(stack *Stack) +{ + array * a=vm::pop(Stack); +#line 841 "runarray.in" + size_t asize=checkArray(a); + if(asize == 0) + error("cannot pop element from empty array"); + {Stack->push(a->pop()); return;} +} + +// Returns the pop method for an array. +#line 849 "runarray.in" +void arrayPop(stack *Stack) +{ + array * a=vm::pop(Stack); +#line 850 "runarray.in" + {Stack->push(new thunk(new bfunc(arrayPopHelper),a)); return;} +} + +// The helper function for the insert method. +#line 855 "runarray.in" +void arrayInsertHelper(stack *Stack) +{ + array * a=vm::pop(Stack); + array * x=vm::pop(Stack); + Int i=vm::pop(Stack); +#line 856 "runarray.in" + size_t asize=checkArray(a); + checkArray(x); + if(a->cyclic() && asize > 0) i=imod(i,asize); + if(i < 0 || i > (Int) asize) + outOfBounds("inserting",asize,i); + (*a).insert((*a).begin()+i,(*x).begin(),(*x).end()); +} + +// Returns the insert method for an array. +#line 866 "runarray.in" +void arrayInsert(stack *Stack) +{ + array * a=vm::pop(Stack); +#line 867 "runarray.in" + {Stack->push(new thunk(new bfunc(arrayInsertHelper),a)); return;} +} + +// Returns the delete method for an array. +#line 872 "runarray.in" +void arrayDelete(stack *Stack) +{ + array * a=vm::pop(Stack); +#line 873 "runarray.in" + {Stack->push(new thunk(new bfunc(arrayDeleteHelper),a)); return;} +} + +#line 877 "runarray.in" +void arrayAlias(stack *Stack) +{ + array * b=vm::pop(Stack); + array * a=vm::pop(Stack); +#line 878 "runarray.in" + {Stack->push(a==b); return;} +} + +// Return array formed by indexing array a with elements of integer array b +#line 883 "runarray.in" +void arrayIntArray(stack *Stack) +{ + array * b=vm::pop(Stack); + array * a=vm::pop(Stack); +#line 884 "runarray.in" + size_t asize=checkArray(a); + size_t bsize=checkArray(b); + array *r=new array(bsize); + bool cyclic=a->cyclic(); + for(size_t i=0; i < bsize; i++) { + Int index=read(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(r); return;} +} + +// returns the complement of the integer array a in {0,2,...,n-1}, +// so that b[complement(a,b.length)] yields the complement of b[a]. +#line 902 "runarray.in" +// Intarray* complement(Intarray *a, Int n); +void gen_runarray32(stack *Stack) +{ + Int n=vm::pop(Stack); + Intarray * a=vm::pop(Stack); +#line 903 "runarray.in" + size_t asize=checkArray(a); + array *r=new array(0); + bool *keep=new bool[n]; + for(Int i=0; i < n; ++i) keep[i]=true; + for(size_t i=0; i < asize; ++i) { + Int j=read(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(r); return;} +} + +// Generate the sequence {f(i) : i=0,1,...n-1} given a function f and integer n +#line 920 "runarray.in" +void arraySequence(stack *Stack) +{ + Int n=vm::pop(Stack); + callable * f=vm::pop(Stack); +#line 921 "runarray.in" + if(n < 0) n=0; + array *a=new array(n); + for(Int i=0; i < n; ++i) { + Stack->push(i); + f->call(Stack); + (*a)[i]=pop(Stack); + } + {Stack->push(a); return;} +} + +// Return the array {0,1,...n-1} +#line 933 "runarray.in" +// Intarray* sequence(Int n); +void gen_runarray34(stack *Stack) +{ + Int n=vm::pop(Stack); +#line 934 "runarray.in" + if(n < 0) n=0; + array *a=new array(n); + for(Int i=0; i < n; ++i) { + (*a)[i]=i; + } + {Stack->push(a); return;} +} + +// Apply a function to each element of an array +#line 944 "runarray.in" +void arrayFunction(stack *Stack) +{ + array * a=vm::pop(Stack); + callable * f=vm::pop(Stack); +#line 945 "runarray.in" + size_t size=checkArray(a); + array *b=new array(size); + for(size_t i=0; i < size; ++i) { + Stack->push((*a)[i]); + f->call(Stack); + (*b)[i]=pop(Stack); + } + {Stack->push(b); return;} +} + +#line 956 "runarray.in" +void arraySort(stack *Stack) +{ + bool stable=vm::pop(Stack,true); + callable * less=vm::pop(Stack); + array * a=vm::pop(Stack); +#line 957 "runarray.in" + array *c=copyArray(a); + compareFunc=less; + FuncStack=Stack; + if(stable) stable_sort(c->begin(),c->end(),compareFunction); + else sort(c->begin(),c->end(),compareFunction); + {Stack->push(c); return;} +} + +#line 966 "runarray.in" +void arraySearch(stack *Stack) +{ + callable * less=vm::pop(Stack); + item key=vm::pop(Stack); + array * a=vm::pop(Stack); +#line 967 "runarray.in" + size_t size=a->size(); + compareFunc=less; + FuncStack=Stack; + if(size == 0 || compareFunction(key,(*a)[0])) {Stack->push(-1); return;} + size_t u=size-1; + if(!compareFunction(key,(*a)[u])) {Stack->push(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(Intcast(i)); return;} + else l=i+1; + } + {Stack->push(0); return;} +} + +#line 985 "runarray.in" +// bool all(boolarray *a); +void gen_runarray38(stack *Stack) +{ + boolarray * a=vm::pop(Stack); +#line 986 "runarray.in" + size_t size=checkArray(a); + bool c=true; + for(size_t i=0; i < size; i++) + if(!get((*a)[i])) {c=false; break;} + {Stack->push(c); return;} +} + +#line 994 "runarray.in" +// boolarray* !(boolarray* a); +void gen_runarray39(stack *Stack) +{ + boolarray* a=vm::pop(Stack); +#line 995 "runarray.in" + size_t size=checkArray(a); + array *c=new array(size); + for(size_t i=0; i < size; i++) + (*c)[i]=!read(a,i); + {Stack->push(c); return;} +} + +#line 1003 "runarray.in" +// Int sum(boolarray *a); +void gen_runarray40(stack *Stack) +{ + boolarray * a=vm::pop(Stack); +#line 1004 "runarray.in" + size_t size=checkArray(a); + Int sum=0; + for(size_t i=0; i < size; i++) + sum += read(a,i) ? 1 : 0; + {Stack->push(sum); return;} +} + +#line 1012 "runarray.in" +void arrayConcat(stack *Stack) +{ + array * a=vm::pop(Stack); +#line 1013 "runarray.in" + // a is an array of arrays to be concatenated together. + // The signature is + // T[] concat(... T[][] a); + + size_t numArgs=checkArray(a); + size_t resultSize=0; + for (size_t i=0; i < numArgs; ++i) { + resultSize += checkArray(a->read(i)); + } + + array *result=new array(resultSize); + + size_t ri=0; + for (size_t i=0; i < numArgs; ++i) { + array *arg=a->read(i); + size_t size=checkArray(arg); + + for (size_t j=0; j < size; ++j) { + (*result)[ri]=(*arg)[j]; + ++ri; + } + } + + {Stack->push(result); return;} +} + +#line 1040 "runarray.in" +void array2Transpose(stack *Stack) +{ + array * a=vm::pop(Stack); +#line 1041 "runarray.in" + size_t asize=checkArray(a); + array *c=new array(0); + size_t csize=0; + for(size_t i=0; i < asize; i++) { + size_t ip=i+1; + array *ai=read(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(c,j); + if(checkArray(cj) < ip) cj->resize(ip); + (*cj)[i]=(*ai)[j]; + } + } + } + {Stack->push(c); return;} +} + +// a is a rectangular 3D array; perm is an Int array indicating the type of +// permutation (021 or 120, etc; original is 012). +// Transpose by sending respective members to the permutated locations: +// return the array obtained by putting a[i][j][k] into position perm{ijk}. +#line 1070 "runarray.in" +void array3Transpose(stack *Stack) +{ + array * perm=vm::pop(Stack); + array * a=vm::pop(Stack); +#line 1071 "runarray.in" + const size_t DIM=3; + + if(checkArray(perm) != DIM) { + ostringstream buf; + buf << "permutation array must have length " << DIM; + error(buf); + } + + size_t* size=new size_t[DIM]; + for(size_t i=0; i < DIM; ++i) size[i]=DIM; + + for(size_t i=0; i < DIM; ++i) { + Int p=read(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(a,0); + size[1]=checkArray(a0); + array *a00=read(a0,0); + size[2]=checkArray(a00); + for(size_t i=0; i < isize; i++) { + array *ai=read(a,i); + size_t jsize=checkArray(ai); + if(jsize != size[1]) error(rectangular); + for(size_t j=0; j < jsize; j++) { + array *aij=read(ai,j); + if(checkArray(aij) != size[2]) error(rectangular); + } + } + + size_t perm0=(size_t) read(perm,0); + size_t perm1=(size_t) read(perm,1); + size_t perm2=(size_t) read(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(a,i[0]); + for(i[1]=0; i[1] < size[1]; ++i[1]) { + array *a1=read(a0,i[1]); + for(i[2]=0; i[2] < size[2]; ++i[2]) { + array *c0=read(c,i[perm0]); + array *c1=read(c0,i[perm1]); + (*c1)[i[perm2]]=read(a1,i[2]); + } + } + } + + delete[] i; + delete[] size; + + {Stack->push(c); return;} +} + +// Find the index of the nth true value in a boolean array or -1 if not found. +// If n is negative, search backwards. +#line 1155 "runarray.in" +// Int find(boolarray *a, Int n=1); +void gen_runarray44(stack *Stack) +{ + Int n=vm::pop(Stack,1); + boolarray * a=vm::pop(Stack); +#line 1156 "runarray.in" + size_t size=checkArray(a); + Int j=-1; + if(n > 0) + for(size_t i=0; i < size; i++) + if(read(a,i)) { + n--; if(n == 0) {j=(Int) i; break;} + } + if(n < 0) + for(size_t i=size; i > 0;) + if(read(a,--i)) { + n++; if(n == 0) {j=(Int) i; break;} + } + {Stack->push(j); return;} +} + +// Find all indices of true values in a boolean array. +#line 1173 "runarray.in" +// Intarray* findall(boolarray *a); +void gen_runarray45(stack *Stack) +{ + boolarray * a=vm::pop(Stack); +#line 1174 "runarray.in" + size_t size=checkArray(a); + array *b=new array(0); + for(size_t i=0; i < size; i++) { + if(read(a,i)) { + b->push((Int) i); + } + } + {Stack->push(b); return;} +} + +// construct vector obtained by replacing those elements of b for which the +// corresponding elements of a are false by the corresponding element of c. +#line 1187 "runarray.in" +void arrayConditional(stack *Stack) +{ + array * c=vm::pop(Stack); + array * b=vm::pop(Stack); + array * a=vm::pop(Stack); +#line 1188 "runarray.in" + size_t size=checkArray(a); + array *r=new array(size); + if(b && c) { + checkArrays(a,b); + checkArrays(b,c); + for(size_t i=0; i < size; i++) + (*r)[i]=read(a,i) ? (*b)[i] : (*c)[i]; + } else { + r->clear(); + if(b) { + checkArrays(a,b); + for(size_t i=0; i < size; i++) + if(read(a,i)) r->push((*b)[i]); + } else if(c) { + checkArrays(a,c); + for(size_t i=0; i < size; i++) + if(!read(a,i)) r->push((*c)[i]); + } + } + {Stack->push(r); return;} +} + +// Return an n x n identity matrix. +#line 1212 "runarray.in" +// realarray2* identity(Int n); +void gen_runarray47(stack *Stack) +{ + Int n=vm::pop(Stack); +#line 1213 "runarray.in" + {Stack->push(Identity(n)); return;} +} + +// Return the inverse of an n x n matrix a using Gauss-Jordan elimination. +#line 1218 "runarray.in" +// realarray2* inverse(realarray2 *a); +void gen_runarray48(stack *Stack) +{ + realarray2 * a=vm::pop(Stack); +#line 1219 "runarray.in" + size_t n=checkArray(a); + double *A; + copyArray2C(A,a,true,0,NoGC); + inverse(A,n); + a=copyCArray2(n,n,A); + delete[] A; + {Stack->push(a); return;} +} + +// Solve the linear equation ax=b by LU decomposition, returning the +// solution x, where a is an n x n matrix and b is an array of length n. +// If no solution exists, return an empty array. +#line 1232 "runarray.in" +// realarray* solve(realarray2 *a, realarray *b, bool warn=true); +void gen_runarray49(stack *Stack) +{ + bool warn=vm::pop(Stack,true); + realarray * b=vm::pop(Stack); + realarray2 * a=vm::pop(Stack); +#line 1233 "runarray.in" + size_t n=checkArray(a); + + if(n == 0) {Stack->push(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(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(x); return;} +} + +// Solve the linear equation ax=b by LU decomposition, returning the +// solution x, where a is an n x n matrix and b is an n x m matrix. +// If no solution exists, return an empty array. +#line 1285 "runarray.in" +// realarray2* solve(realarray2 *a, realarray2 *b, bool warn=true); +void gen_runarray50(stack *Stack) +{ + bool warn=vm::pop(Stack,true); + realarray2 * b=vm::pop(Stack); + realarray2 * a=vm::pop(Stack); +#line 1286 "runarray.in" + size_t n=checkArray(a); + + if(n == 0) {Stack->push(new array(0)); return;} + + if(checkArray(b) != n) error(incommensurate); + size_t m=checkArray(read(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(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(x); return;} +} + +// Compute the determinant of an n x n matrix. +#line 1348 "runarray.in" +// real determinant(realarray2 *a); +void gen_runarray51(stack *Stack) +{ + realarray2 * a=vm::pop(Stack); +#line 1349 "runarray.in" + real *A; + copyArray2C(A,a); + size_t n=checkArray(a); + + real det=LUdecompose(A,n,NULL,false); + size_t n1=n+1; + for(size_t i=0; i < n; ++i) + det *= A[i*n1]; + + delete[] A; + + {Stack->push(det); return;} +} + +#line 1364 "runarray.in" +// realarray* *(realarray2 *a, realarray *b); +void gen_runarray52(stack *Stack) +{ + realarray * b=vm::pop(Stack); + realarray2 * a=vm::pop(Stack); +#line 1365 "runarray.in" + size_t n=checkArray(a); + size_t m=checkArray(b); + array *c=new array(n); + real *B; + copyArrayC(B,b); + for(size_t i=0; i < n; ++i) { + array *ai=read(a,i); + if(checkArray(ai) != m) error(incommensurate); + real sum=0.0; + for(size_t j=0; j < m; ++j) + sum += read(ai,j)*B[j]; + (*c)[i]=sum; + } + delete[] B; + {Stack->push(c); return;} +} + +#line 1383 "runarray.in" +// realarray* *(realarray *a, realarray2 *b); +void gen_runarray53(stack *Stack) +{ + realarray2 * b=vm::pop(Stack); + realarray * a=vm::pop(Stack); +#line 1384 "runarray.in" + size_t n=checkArray(a); + if(n != checkArray(b)) error(incommensurate); + real *A; + copyArrayC(A,a); + + array **B=new array*[n]; + array *bk=read(b,0); + B[0]=bk; + size_t m=bk->size(); + for(size_t k=1; k < n; k++) { + array *bk=read(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(B[k],i); + (*c)[i]=sum; + } + delete[] B; + delete[] A; + {Stack->push(c); return;} +} + +#line 1412 "runarray.in" +// Intarray2* *(Intarray2 *a, Intarray2 *b); +void gen_runarray54(stack *Stack) +{ + Intarray2 * b=vm::pop(Stack); + Intarray2 * a=vm::pop(Stack); +#line 1413 "runarray.in" + {Stack->push(mult(a,b)); return;} +} + +#line 1417 "runarray.in" +// realarray2* *(realarray2 *a, realarray2 *b); +void gen_runarray55(stack *Stack) +{ + realarray2 * b=vm::pop(Stack); + realarray2 * a=vm::pop(Stack); +#line 1418 "runarray.in" + {Stack->push(mult(a,b)); return;} +} + +#line 1422 "runarray.in" +// pairarray2* *(pairarray2 *a, pairarray2 *b); +void gen_runarray56(stack *Stack) +{ + pairarray2 * b=vm::pop(Stack); + pairarray2 * a=vm::pop(Stack); +#line 1423 "runarray.in" + {Stack->push(mult(a,b)); return;} +} + +#line 1427 "runarray.in" +// triple *(realarray2 *t, triple v); +void gen_runarray57(stack *Stack) +{ + triple v=vm::pop(Stack); + realarray2 * t=vm::pop(Stack); +#line 1428 "runarray.in" + {Stack->push(*t*v); return;} +} + +#line 1432 "runarray.in" +// realarray2* AtA(realarray2 *a); +void gen_runarray58(stack *Stack) +{ + realarray2 * a=vm::pop(Stack); +#line 1433 "runarray.in" + {Stack->push(AtA(a)); return;} +} + +#line 1437 "runarray.in" +// pair project(triple v, realarray2 *t); +void gen_runarray59(stack *Stack) +{ + realarray2 * t=vm::pop(Stack); + triple v=vm::pop(Stack); +#line 1438 "runarray.in" + size_t n=checkArray(t); + if(n != 4) error(incommensurate); + array *t0=read(t,0); + array *t1=read(t,1); + array *t3=read(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(t3,0)*x+read(t3,1)*y+read(t3,2)*z+ + read(t3,3); + if(f == 0.0) dividebyzero(); + f=1.0/f; + + {Stack->push(pair((read(t0,0)*x+read(t0,1)*y+read(t0,2)*z+ + read(t0,3))*f, + (read(t1,0)*x+read(t1,1)*y+read(t1,2)*z+ + read(t1,3))*f)); return;} +} + +// Compute the dot product of vectors a and b. +#line 1463 "runarray.in" +// real dot(realarray *a, realarray *b); +void gen_runarray60(stack *Stack) +{ + realarray * b=vm::pop(Stack); + realarray * a=vm::pop(Stack); +#line 1464 "runarray.in" + size_t n=checkArrays(a,b); + real sum=0.0; + for(size_t i=0; i < n; ++i) + sum += read(a,i)*read(b,i); + {Stack->push(sum); return;} +} + +// Compute the complex dot product of vectors a and b. +#line 1473 "runarray.in" +// pair dot(pairarray *a, pairarray *b); +void gen_runarray61(stack *Stack) +{ + pairarray * b=vm::pop(Stack); + pairarray * a=vm::pop(Stack); +#line 1474 "runarray.in" + size_t n=checkArrays(a,b); + pair sum=zero; + for(size_t i=0; i < n; ++i) + sum += read(a,i)*conj(read(b,i)); + {Stack->push(sum); return;} +} + +// Solve the problem L\inv f, where f is an n vector and L is the n x n matrix +// +// [ b[0] c[0] a[0] ] +// [ a[1] b[1] c[1] ] +// [ a[2] b[2] c[2] ] +// [ ... ] +// [ c[n-1] a[n-1] b[n-1] ] +#line 1489 "runarray.in" +// realarray* tridiagonal(realarray *a, realarray *b, realarray *c, realarray *f); +void gen_runarray62(stack *Stack) +{ + realarray * f=vm::pop(Stack); + realarray * c=vm::pop(Stack); + realarray * b=vm::pop(Stack); + realarray * a=vm::pop(Stack); +#line 1490 "runarray.in" + size_t n=checkArrays(a,b); + checkEqual(n,checkArray(c)); + checkEqual(n,checkArray(f)); + + array *up=new array(n); + array& u=*up; + + if(n == 0) {Stack->push(up); return;} + + // Special case: zero Dirichlet boundary conditions + if(read(a,0) == 0.0 && read(c,n-1) == 0.0) { + real temp=read(b,0); + if(temp == 0.0) dividebyzero(); + temp=1.0/temp; + + real *work=new real[n]; + u[0]=read(f,0)*temp; + work[0]=-read(c,0)*temp; + + for(size_t i=1; i < n; i++) { + real temp=(read(b,i)+read(a,i)*work[i-1]); + if(temp == 0.0) {delete[] work; dividebyzero();} + temp=1.0/temp; + u[i]=(read(f,i)-read(a,i)*read(u,i-1))*temp; + work[i]=-read(c,i)*temp; + } + + for(size_t i=n-1; i >= 1; i--) + u[i-1]=read(u,i-1)+work[i-1]*read(u,i); + + delete[] work; + {Stack->push(up); return;} + } + + real binv=read(b,0); + if(binv == 0.0) dividebyzero(); + binv=1.0/binv; + + if(n == 1) {u[0]=read(f,0)*binv; {Stack->push(up); return;}} + if(n == 2) { + real factor=(read(b,0)*read(b,1)- + read(a,0)*read(c,1)); + if(factor== 0.0) dividebyzero(); + factor=1.0/factor; + real temp=(read(b,0)*read(f,1)- + read(c,1)*read(f,0))*factor; + u[0]=(read(b,1)*read(f,0)- + read(a,0)*read(f,1))*factor; + u[1]=temp; + {Stack->push(up); return;} + } + + real *gamma=new real[n-2]; + real *delta=new real[n-2]; + + gamma[0]=read(c,0)*binv; + delta[0]=read(a,0)*binv; + u[0]=read(f,0)*binv; + real beta=read(c,n-1); + real fn=read(f,n-1)-beta*read(u,0); + real alpha=read(b,n-1)-beta*delta[0]; + + for(size_t i=1; i <= n-3; i++) { + real alphainv=read(b,i)-read(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(c,i)*alphainv; + u[i]=(read(f,i)-read(a,i)*read(u,i-1))*alphainv; + fn -= beta*read(u,i); + delta[i]=-read(a,i)*delta[i-1]*alphainv; + alpha -= beta*delta[i]; + } + + real alphainv=read(b,n-2)-read(a,n-2)*gamma[n-3]; + if(alphainv == 0.0) {delete[] gamma; delete[] delta; dividebyzero();} + alphainv=1.0/alphainv; + u[n-2]=(read(f,n-2)-read(a,n-2)*read(u,n-3)) + *alphainv; + beta=read(a,n-1)-beta*gamma[n-3]; + real dnm1=(read(c,n-2)-read(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(u,n-2))/temp; + u[n-2]=read(u,n-2)-dnm1*temp; + + for(size_t i=n-2; i >= 1; i--) + u[i-1]=read(u,i-1)-gamma[i-1]*read(u,i)-delta[i-1]*temp; + + delete[] delta; + delete[] gamma; + + {Stack->push(up); return;} +} + +// Root solve by Newton-Raphson +#line 1587 "runarray.in" +// real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x, bool verbose=false); +void gen_runarray63(stack *Stack) +{ + bool verbose=vm::pop(Stack,false); + real x=vm::pop(Stack); + callableReal * fprime=vm::pop(Stack); + callableReal * f=vm::pop(Stack); + Int iterations=vm::pop(Stack,100); +#line 1589 "runarray.in" + static const real fuzz=1000.0*DBL_EPSILON; + Int i=0; + size_t oldPrec=0; + if(verbose) + oldPrec=cout.precision(DBL_DIG); + + real diff=DBL_MAX; + real lastdiff; + do { + real x0=x; + + Stack->push(x); + fprime->call(Stack); + real dfdx=pop(Stack); + + if(dfdx == 0.0) { + x=DBL_MAX; + break; + } + + Stack->push(x); + f->call(Stack); + real fx=pop(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(x); return;} +} + +// Root solve by Newton-Raphson bisection +// cf. routine rtsafe (Press et al., Numerical Recipes, 1991). +#line 1635 "runarray.in" +// real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x1, real x2, bool verbose=false); +void gen_runarray64(stack *Stack) +{ + bool verbose=vm::pop(Stack,false); + real x2=vm::pop(Stack); + real x1=vm::pop(Stack); + callableReal * fprime=vm::pop(Stack); + callableReal * f=vm::pop(Stack); + Int iterations=vm::pop(Stack,100); +#line 1637 "runarray.in" + static const real fuzz=1000.0*DBL_EPSILON; + size_t oldPrec=0; + if(verbose) + oldPrec=cout.precision(DBL_DIG); + + Stack->push(x1); + f->call(Stack); + real f1=pop(Stack); + if(f1 == 0.0) {Stack->push(x1); return;} + + Stack->push(x2); + f->call(Stack); + real f2=pop(Stack); + if(f2 == 0.0) {Stack->push(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(Stack); + + Stack->push(x); + fprime->call(Stack); + real dy=pop(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(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(x); return;} + } + if(fabs(dx) < fuzz*fabs(x)) {Stack->push(x); return;} + + Stack->push(x); + f->call(Stack); + y=pop(Stack); + + Stack->push(x); + fprime->call(Stack); + dy=pop(Stack); + + if(y < 0.0) x1=x; + else x2=x; + } + if(verbose) + cout.precision(oldPrec); + {Stack->push((j == iterations) ? DBL_MAX : x); return;} +} + +// Find a root for the specified continuous (but not necessarily +// differentiable) function. Whatever value t is returned, it is guaranteed +// that t is within [a, b] and within tolerance of a sign change. +// An error is thrown if fa and fb are both positive or both negative. +// +// In this implementation, the binary search is interleaved +// with a modified version of quadratic interpolation. +// This is a C++ port of the Asymptote routine written by Charles Staats III. +#line 1723 "runarray.in" +// real _findroot(callableReal *f, real a, real b, real tolerance, real fa, real fb); +void gen_runarray65(stack *Stack) +{ + real fb=vm::pop(Stack); + real fa=vm::pop(Stack); + real tolerance=vm::pop(Stack); + real b=vm::pop(Stack); + real a=vm::pop(Stack); + callableReal * f=vm::pop(Stack); +#line 1725 "runarray.in" + if(fa == 0.0) {Stack->push(a); return;} + if(fb == 0.0) {Stack->push(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(Stack); + if(ft == 0.0) {Stack->push(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(Stack); + + if(ft == 0.0) {Stack->push(t); return;} + } + } + + if(ft > 0.0) { + b=t; + fb=ft; + } else if(ft < 0.0) { + a=t; + fa=ft; + } + } + {Stack->push(a-(b-a)/(fb-fa)*fa); return;} +} + +#line 1817 "runarray.in" +// real simpson(callableReal *f, real a, real b, real acc=DBL_EPSILON, real dxmax=0); +void gen_runarray66(stack *Stack) +{ + real dxmax=vm::pop(Stack,0); + real acc=vm::pop(Stack,DBL_EPSILON); + real b=vm::pop(Stack); + real a=vm::pop(Stack); + callableReal * f=vm::pop(Stack); +#line 1819 "runarray.in" + real integral; + if(dxmax <= 0) dxmax=fabs(b-a); + callable *oldFunc=Func; + Func=f; + FuncStack=Stack; + if(!simpson(integral,wrapFunction,a,b,acc,dxmax)) + error("nesting capacity exceeded in simpson"); + Func=oldFunc; + {Stack->push(integral); return;} +} + +// Compute the fast Fourier transform of a pair array +#line 1832 "runarray.in" +// pairarray* fft(pairarray *a, Int sign=1); +void gen_runarray67(stack *Stack) +{ + Int sign=vm::pop(Stack,1); + pairarray * a=vm::pop(Stack); +#line 1833 "runarray.in" +#ifdef HAVE_LIBFFTW3 + unsigned n=(unsigned) checkArray(a); + array *c=new array(n); + if(n) { + Complex *f=utils::ComplexAlign(n); + fftwpp::fft1d Forward(n,intcast(sign),f); + + for(size_t i=0; i < n; i++) { + pair z=read(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(c); return;} +} + +// Compute the fast Fourier transform of a 2D pair array +#line 1863 "runarray.in" +// pairarray2* fft(pairarray2 *a, Int sign=1); +void gen_runarray68(stack *Stack) +{ + Int sign=vm::pop(Stack,1); + pairarray2 * a=vm::pop(Stack); +#line 1864 "runarray.in" +#ifdef HAVE_LIBFFTW3 + size_t n=checkArray(a); + size_t m=n == 0 ? 0 : checkArray(read(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(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(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(c); return;} +} + +// Compute the fast Fourier transform of a 3D pair array +#line 1909 "runarray.in" +// pairarray3* fft(pairarray3 *a, Int sign=1); +void gen_runarray69(stack *Stack) +{ + Int sign=vm::pop(Stack,1); + pairarray3 * a=vm::pop(Stack); +#line 1910 "runarray.in" +#ifdef HAVE_LIBFFTW3 + size_t n=checkArray(a); + array *a0=read(a,0); + size_t m=n == 0 ? 0 : checkArray(a0); + size_t l=m == 0 ? 0 : checkArray(read(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(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(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(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(c); return;} +} + +#line 1967 "runarray.in" +// Intarray2* triangulate(pairarray *z); +void gen_runarray70(stack *Stack) +{ + pairarray * z=vm::pop(Stack); +#line 1968 "runarray.in" + size_t nv=checkArray(z); +// Call robust version of Gilles Dumoulin's port of Paul Bourke's +// triangulation code. + + XYZ *pxyz=new XYZ[nv+3]; + ITRIANGLE *V=new ITRIANGLE[4*nv]; + + for(size_t i=0; i < nv; ++i) { + pair w=read(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(t); return;} +} + +#line 2002 "runarray.in" +// real norm(realarray *a); +void gen_runarray71(stack *Stack) +{ + realarray * a=vm::pop(Stack); +#line 2003 "runarray.in" + size_t n=checkArray(a); + real M=0.0; + for(size_t i=0; i < n; ++i) { + real x=fabs(vm::read(a,i)); + if(x > M) M=x; + } + {Stack->push(M); return;} +} + +#line 2013 "runarray.in" +// real norm(realarray2 *a); +void gen_runarray72(stack *Stack) +{ + realarray2 * a=vm::pop(Stack); +#line 2014 "runarray.in" + size_t n=checkArray(a); + real M=0.0; + for(size_t i=0; i < n; ++i) { + vm::array *ai=vm::read(a,i); + size_t m=checkArray(ai); + for(size_t j=0; j < m; ++j) { + real a=fabs(vm::read(ai,j)); + if(a > M) M=a; + } + } + {Stack->push(M); return;} +} + +#line 2028 "runarray.in" +// real norm(triplearray2 *a); +void gen_runarray73(stack *Stack) +{ + triplearray2 * a=vm::pop(Stack); +#line 2029 "runarray.in" + size_t n=checkArray(a); + real M=0.0; + for(size_t i=0; i < n; ++i) { + vm::array *ai=vm::read(a,i); + size_t m=checkArray(ai); + for(size_t j=0; j < m; ++j) { + real a=vm::read(ai,j).abs2(); + if(a > M) M=a; + } + } + {Stack->push(sqrt(M)); return;} +} + +#line 2043 "runarray.in" +// real change2(triplearray2 *a); +void gen_runarray74(stack *Stack) +{ + triplearray2 * a=vm::pop(Stack); +#line 2044 "runarray.in" + size_t n=checkArray(a); + if(n == 0) {Stack->push(0.0); return;} + + vm::array *a0=vm::read(a,0); + size_t m=checkArray(a0); + if(m == 0) {Stack->push(0.0); return;} + triple a00=vm::read(a0,0); + real M=0.0; + + for(size_t i=0; i < n; ++i) { + vm::array *ai=vm::read(a,i); + size_t m=checkArray(ai); + for(size_t j=0; j < m; ++j) { + real a=(vm::read(ai,j)-a00).abs2(); + if(a > M) M=a; + } + } + {Stack->push(M); return;} +} + +#line 2065 "runarray.in" +// triple minbezier(triplearray2 *P, triple b); +void gen_runarray75(stack *Stack) +{ + triple b=vm::pop(Stack); + triplearray2 * P=vm::pop(Stack); +#line 2066 "runarray.in" + size_t N; + real *A=copyTripleArray2Components(P,N); + bound_double *B=bounddouble(N); + b=triple(B(A,::min,b.getx(),Fuzz*norm(A,N),maxdepth), + B(A+N,::min,b.gety(),Fuzz*norm(A+N,N),maxdepth), + B(A+2*N,::min,b.getz(),Fuzz*norm(A+2*N,N),maxdepth)); + delete[] A; + {Stack->push(b); return;} +} + +#line 2077 "runarray.in" +// triple maxbezier(triplearray2 *P, triple b); +void gen_runarray76(stack *Stack) +{ + triple b=vm::pop(Stack); + triplearray2 * P=vm::pop(Stack); +#line 2078 "runarray.in" + size_t N; + real *A=copyTripleArray2Components(P,N); + bound_double *B=bounddouble(N); + b=triple(B(A,::max,b.getx(),Fuzz*norm(A,N),maxdepth), + B(A+N,::max,b.gety(),Fuzz*norm(A+N,N),maxdepth), + B(A+2*N,::max,b.getz(),Fuzz*norm(A+2*N,N),maxdepth)); + delete[] A; + {Stack->push(b); return;} +} + +#line 2089 "runarray.in" +// pair minratio(triplearray2 *P, pair b); +void gen_runarray77(stack *Stack) +{ + pair b=vm::pop(Stack); + triplearray2 * P=vm::pop(Stack); +#line 2090 "runarray.in" + size_t N; + triple *A=copyTripleArray2C(P,N); + real fuzz=Fuzz*norm(A,N); + bound_triple *B=boundtriple(N); + b=pair(B(A,::min,xratio,b.getx(),fuzz,maxdepth), + B(A,::min,yratio,b.gety(),fuzz,maxdepth)); + delete[] A; + {Stack->push(b); return;} +} + +#line 2101 "runarray.in" +// pair maxratio(triplearray2 *P, pair b); +void gen_runarray78(stack *Stack) +{ + pair b=vm::pop(Stack); + triplearray2 * P=vm::pop(Stack); +#line 2102 "runarray.in" + size_t N; + triple *A=copyTripleArray2C(P,N); + bound_triple *B=boundtriple(N); + real fuzz=Fuzz*norm(A,N); + b=pair(B(A,::max,xratio,b.getx(),fuzz,maxdepth), + B(A,::max,yratio,b.gety(),fuzz,maxdepth)); + delete[] A; + {Stack->push(b); return;} +} + +#line 2113 "runarray.in" +// realarray* _projection(); +void gen_runarray79(stack *Stack) +{ +#line 2114 "runarray.in" +#ifdef HAVE_GL + array *a=new array(14); + gl::projection P=gl::camera(); + size_t k=0; + (*a)[k++]=P.orthographic ? 1.0 : 0.0; + + triple camera=P.camera; + (*a)[k++]=camera.getx(); + (*a)[k++]=camera.gety(); + (*a)[k++]=camera.getz(); + + triple up=P.up; + (*a)[k++]=up.getx(); + (*a)[k++]=up.gety(); + (*a)[k++]=up.getz(); + + triple target=P.target; + (*a)[k++]=target.getx(); + (*a)[k++]=target.gety(); + (*a)[k++]=target.getz(); + + (*a)[k++]=P.zoom; + (*a)[k++]=P.angle; + + (*a)[k++]=P.viewportshift.getx(); + (*a)[k++]=P.viewportshift.gety(); +#endif + {Stack->push(new array(0)); return;} +} + +} // namespace run + +namespace trans { + +void gen_runarray_venv(venv &ve) +{ +#line 593 "runarray.in" + REGISTER_BLTIN(run::emptyArray,"emptyArray"); +#line 599 "runarray.in" + REGISTER_BLTIN(run::newDeepArray,"newDeepArray"); +#line 621 "runarray.in" + REGISTER_BLTIN(run::newInitializedArray,"newInitializedArray"); +#line 636 "runarray.in" + REGISTER_BLTIN(run::newAppendedArray,"newAppendedArray"); +#line 652 "runarray.in" + REGISTER_BLTIN(run::copyArrayValue,"copyArrayValue"); +#line 664 "runarray.in" + REGISTER_BLTIN(run::copyArray,"copyArray"); +#line 676 "runarray.in" + REGISTER_BLTIN(run::arrayRead,"arrayRead"); +#line 688 "runarray.in" + REGISTER_BLTIN(run::arraySliceRead,"arraySliceRead"); +#line 695 "runarray.in" + REGISTER_BLTIN(run::arraySliceReadToEnd,"arraySliceReadToEnd"); +#line 703 "runarray.in" + REGISTER_BLTIN(run::arrayArrayRead,"arrayArrayRead"); +#line 712 "runarray.in" + REGISTER_BLTIN(run::arrayWrite,"arrayWrite"); +#line 729 "runarray.in" + REGISTER_BLTIN(run::arraySliceWrite,"arraySliceWrite"); +#line 737 "runarray.in" + REGISTER_BLTIN(run::arraySliceWriteToEnd,"arraySliceWriteToEnd"); +#line 745 "runarray.in" + REGISTER_BLTIN(run::arrayLength,"arrayLength"); +#line 751 "runarray.in" + REGISTER_BLTIN(run::arrayKeys,"arrayKeys"); +#line 766 "runarray.in" + REGISTER_BLTIN(run::arrayCyclicFlag,"arrayCyclicFlag"); +#line 773 "runarray.in" + REGISTER_BLTIN(run::arraySetCyclicFlag,"arraySetCyclicFlag"); +#line 780 "runarray.in" + REGISTER_BLTIN(run::arrayInitializedHelper,"arrayInitializedHelper"); +#line 791 "runarray.in" + REGISTER_BLTIN(run::arrayInitialized,"arrayInitialized"); +#line 797 "runarray.in" + REGISTER_BLTIN(run::arrayCyclicHelper,"arrayCyclicHelper"); +#line 804 "runarray.in" + REGISTER_BLTIN(run::arrayCyclic,"arrayCyclic"); +#line 810 "runarray.in" + REGISTER_BLTIN(run::arrayPushHelper,"arrayPushHelper"); +#line 818 "runarray.in" + REGISTER_BLTIN(run::arrayPush,"arrayPush"); +#line 824 "runarray.in" + REGISTER_BLTIN(run::arrayAppendHelper,"arrayAppendHelper"); +#line 833 "runarray.in" + REGISTER_BLTIN(run::arrayAppend,"arrayAppend"); +#line 839 "runarray.in" + REGISTER_BLTIN(run::arrayPopHelper,"arrayPopHelper"); +#line 848 "runarray.in" + REGISTER_BLTIN(run::arrayPop,"arrayPop"); +#line 854 "runarray.in" + REGISTER_BLTIN(run::arrayInsertHelper,"arrayInsertHelper"); +#line 865 "runarray.in" + REGISTER_BLTIN(run::arrayInsert,"arrayInsert"); +#line 871 "runarray.in" + REGISTER_BLTIN(run::arrayDelete,"arrayDelete"); +#line 877 "runarray.in" + REGISTER_BLTIN(run::arrayAlias,"arrayAlias"); +#line 882 "runarray.in" + REGISTER_BLTIN(run::arrayIntArray,"arrayIntArray"); +#line 900 "runarray.in" + addFunc(ve, run::gen_runarray32, IntArray(), SYM(complement), formal(IntArray(), SYM(a), false, false), formal(primInt(), SYM(n), false, false)); +#line 919 "runarray.in" + REGISTER_BLTIN(run::arraySequence,"arraySequence"); +#line 932 "runarray.in" + addFunc(ve, run::gen_runarray34, IntArray(), SYM(sequence), formal(primInt(), SYM(n), false, false)); +#line 943 "runarray.in" + REGISTER_BLTIN(run::arrayFunction,"arrayFunction"); +#line 956 "runarray.in" + REGISTER_BLTIN(run::arraySort,"arraySort"); +#line 966 "runarray.in" + REGISTER_BLTIN(run::arraySearch,"arraySearch"); +#line 985 "runarray.in" + addFunc(ve, run::gen_runarray38, primBoolean(), SYM(all), formal(booleanArray(), SYM(a), false, false)); +#line 994 "runarray.in" + addFunc(ve, run::gen_runarray39, booleanArray(), SYM_LOGNOT, formal(booleanArray(), SYM(a), false, false)); +#line 1003 "runarray.in" + addFunc(ve, run::gen_runarray40, primInt(), SYM(sum), formal(booleanArray(), SYM(a), false, false)); +#line 1012 "runarray.in" + REGISTER_BLTIN(run::arrayConcat,"arrayConcat"); +#line 1040 "runarray.in" + REGISTER_BLTIN(run::array2Transpose,"array2Transpose"); +#line 1066 "runarray.in" + REGISTER_BLTIN(run::array3Transpose,"array3Transpose"); +#line 1153 "runarray.in" + addFunc(ve, run::gen_runarray44, primInt(), SYM(find), formal(booleanArray(), SYM(a), false, false), formal(primInt(), SYM(n), true, false)); +#line 1172 "runarray.in" + addFunc(ve, run::gen_runarray45, IntArray(), SYM(findall), formal(booleanArray(), SYM(a), false, false)); +#line 1185 "runarray.in" + REGISTER_BLTIN(run::arrayConditional,"arrayConditional"); +#line 1211 "runarray.in" + addFunc(ve, run::gen_runarray47, realArray2(), SYM(identity), formal(primInt(), SYM(n), false, false)); +#line 1217 "runarray.in" + addFunc(ve, run::gen_runarray48, realArray2(), SYM(inverse), formal(realArray2(), SYM(a), false, false)); +#line 1229 "runarray.in" + addFunc(ve, run::gen_runarray49, realArray(), SYM(solve), formal(realArray2(), SYM(a), false, false), formal(realArray(), SYM(b), false, false), formal(primBoolean(), SYM(warn), true, false)); +#line 1282 "runarray.in" + addFunc(ve, run::gen_runarray50, realArray2(), SYM(solve), formal(realArray2(), SYM(a), false, false), formal(realArray2(), SYM(b), false, false), formal(primBoolean(), SYM(warn), true, false)); +#line 1347 "runarray.in" + addFunc(ve, run::gen_runarray51, primReal(), SYM(determinant), formal(realArray2(), SYM(a), false, false)); +#line 1364 "runarray.in" + addFunc(ve, run::gen_runarray52, realArray(), SYM_TIMES, formal(realArray2(), SYM(a), false, false), formal(realArray(), SYM(b), false, false)); +#line 1383 "runarray.in" + addFunc(ve, run::gen_runarray53, realArray(), SYM_TIMES, formal(realArray(), SYM(a), false, false), formal(realArray2(), SYM(b), false, false)); +#line 1412 "runarray.in" + addFunc(ve, run::gen_runarray54, IntArray2(), SYM_TIMES, formal(IntArray2(), SYM(a), false, false), formal(IntArray2(), SYM(b), false, false)); +#line 1417 "runarray.in" + addFunc(ve, run::gen_runarray55, realArray2(), SYM_TIMES, formal(realArray2(), SYM(a), false, false), formal(realArray2(), SYM(b), false, false)); +#line 1422 "runarray.in" + addFunc(ve, run::gen_runarray56, pairArray2(), SYM_TIMES, formal(pairArray2(), SYM(a), false, false), formal(pairArray2(), SYM(b), false, false)); +#line 1427 "runarray.in" + addFunc(ve, run::gen_runarray57, primTriple(), SYM_TIMES, formal(realArray2(), SYM(t), false, false), formal(primTriple(), SYM(v), false, false)); +#line 1432 "runarray.in" + addFunc(ve, run::gen_runarray58, realArray2(), SYM(AtA), formal(realArray2(), SYM(a), false, false)); +#line 1437 "runarray.in" + addFunc(ve, run::gen_runarray59, primPair(), SYM(project), formal(primTriple(), SYM(v), false, false), formal(realArray2(), SYM(t), false, false)); +#line 1462 "runarray.in" + addFunc(ve, run::gen_runarray60, primReal(), SYM(dot), formal(realArray(), SYM(a), false, false), formal(realArray(), SYM(b), false, false)); +#line 1472 "runarray.in" + addFunc(ve, run::gen_runarray61, primPair(), SYM(dot), formal(pairArray(), SYM(a), false, false), formal(pairArray(), SYM(b), false, false)); +#line 1482 "runarray.in" + addFunc(ve, run::gen_runarray62, realArray(), SYM(tridiagonal), formal(realArray(), SYM(a), false, false), formal(realArray(), SYM(b), false, false), formal(realArray(), SYM(c), false, false), formal(realArray(), SYM(f), false, false)); +#line 1586 "runarray.in" + addFunc(ve, run::gen_runarray63, primReal(), SYM(newton), formal(primInt(), SYM(iterations), true, false), formal(realRealFunction(), SYM(f), false, false), formal(realRealFunction(), SYM(fprime), false, false), formal(primReal(), SYM(x), false, false), formal(primBoolean(), SYM(verbose), true, false)); +#line 1633 "runarray.in" + addFunc(ve, run::gen_runarray64, primReal(), SYM(newton), formal(primInt(), SYM(iterations), true, false), formal(realRealFunction(), SYM(f), false, false), formal(realRealFunction(), SYM(fprime), false, false), formal(primReal(), SYM(x1), false, false), formal(primReal(), SYM(x2), false, false), formal(primBoolean(), SYM(verbose), true, false)); +#line 1715 "runarray.in" + addFunc(ve, run::gen_runarray65, primReal(), SYM(_findroot), formal(realRealFunction(), SYM(f), false, false), formal(primReal(), SYM(a), false, false), formal(primReal(), SYM(b), false, false), formal(primReal(), SYM(tolerance), false, false), formal(primReal(), SYM(fa), false, false), formal(primReal(), SYM(fb), false, false)); +#line 1817 "runarray.in" + addFunc(ve, run::gen_runarray66, primReal(), SYM(simpson), formal(realRealFunction(), SYM(f), false, false), formal(primReal(), SYM(a), false, false), formal(primReal(), SYM(b), false, false), formal(primReal(), SYM(acc), true, false), formal(primReal(), SYM(dxmax), true, false)); +#line 1831 "runarray.in" + addFunc(ve, run::gen_runarray67, pairArray(), SYM(fft), formal(pairArray(), SYM(a), false, false), formal(primInt(), SYM(sign), true, false)); +#line 1862 "runarray.in" + addFunc(ve, run::gen_runarray68, pairArray2(), SYM(fft), formal(pairArray2(), SYM(a), false, false), formal(primInt(), SYM(sign), true, false)); +#line 1908 "runarray.in" + addFunc(ve, run::gen_runarray69, pairArray3(), SYM(fft), formal(pairArray3(), SYM(a), false, false), formal(primInt(), SYM(sign), true, false)); +#line 1967 "runarray.in" + addFunc(ve, run::gen_runarray70, IntArray2(), SYM(triangulate), formal(pairArray(), SYM(z), false, false)); +#line 2002 "runarray.in" + addFunc(ve, run::gen_runarray71, primReal(), SYM(norm), formal(realArray(), SYM(a), false, false)); +#line 2013 "runarray.in" + addFunc(ve, run::gen_runarray72, primReal(), SYM(norm), formal(realArray2(), SYM(a), false, false)); +#line 2028 "runarray.in" + addFunc(ve, run::gen_runarray73, primReal(), SYM(norm), formal(tripleArray2(), SYM(a), false, false)); +#line 2043 "runarray.in" + addFunc(ve, run::gen_runarray74, primReal(), SYM(change2), formal(tripleArray2(), SYM(a), false, false)); +#line 2065 "runarray.in" + addFunc(ve, run::gen_runarray75, primTriple(), SYM(minbezier), formal(tripleArray2(), SYM(p), false, false), formal(primTriple(), SYM(b), false, false)); +#line 2077 "runarray.in" + addFunc(ve, run::gen_runarray76, primTriple(), SYM(maxbezier), formal(tripleArray2(), SYM(p), false, false), formal(primTriple(), SYM(b), false, false)); +#line 2089 "runarray.in" + addFunc(ve, run::gen_runarray77, primPair(), SYM(minratio), formal(tripleArray2(), SYM(p), false, false), formal(primPair(), SYM(b), false, false)); +#line 2101 "runarray.in" + addFunc(ve, run::gen_runarray78, primPair(), SYM(maxratio), formal(tripleArray2(), SYM(p), false, false), formal(primPair(), SYM(b), false, false)); +#line 2113 "runarray.in" + addFunc(ve, run::gen_runarray79, realArray(), SYM(_projection)); +} + +} // namespace trans -- cgit v1.2.3