diff options
Diffstat (limited to 'Build/source/utils/asymptote/runarray.cc')
-rw-r--r-- | Build/source/utils/asymptote/runarray.cc | 727 |
1 files changed, 365 insertions, 362 deletions
diff --git a/Build/source/utils/asymptote/runarray.cc b/Build/source/utils/asymptote/runarray.cc index 5478028e4e7..05c7f0a1a7c 100644 --- a/Build/source/utils/asymptote/runarray.cc +++ b/Build/source/utils/asymptote/runarray.cc @@ -71,7 +71,7 @@ function *realRealFunction(); #ifdef HAVE_LIBFFTW3 #include "fftw++.h" -static const char *rectangular="matrix must be rectangular"; + static const char *rectangular="matrix must be rectangular"; #else static const char *installFFTW= "Please install fftw3, run ./configure, and recompile"; @@ -113,7 +113,7 @@ void outOfBounds(const char *op, size_t len, Int n) error(buf); } -inline item& arrayRead(array *a, Int n) +inline item& arrayRead(array *a, Int n) { size_t len=checkArray(a); bool cyclic=a->cyclic(); @@ -126,7 +126,7 @@ inline item& arrayRead(array *a, Int n) static array* deepArray(Int depth, Int *dims) { assert(depth > 0); - + if (depth == 1) { return new array(dims[0]); } else { @@ -204,7 +204,7 @@ array *copyArray(array *a) { size_t size=checkArray(a); array *c=new array(size); - for(size_t i=0; i < size; i++) + for(size_t i=0; i < size; i++) (*c)[i]=(*a)[i]; return c; } @@ -218,7 +218,7 @@ array *copyArray2(array *a) size_t aisize=checkArray(ai); array *ci=new array(aisize); (*c)[i]=ci; - for(size_t j=0; j < aisize; j++) + for(size_t j=0; j < aisize; j++) (*ci)[j]=(*ai)[j]; } return c; @@ -229,12 +229,12 @@ double *copyTripleArray2Components(array *a, size_t &N, GCPlacement placement) size_t n=checkArray(a); N=0; for(size_t i=0; i < n; i++) - N += checkArray(read<array*>(a,i)); - + N += checkArray(read<array*>(a,i)); + double *A=(placement == NoGC) ? new double [3*N] : new(placement) double[3*N]; double *p=A; - + for(size_t i=0; i < n; i++) { array *ai=read<array*>(a,i); size_t m=checkArray(ai); @@ -254,12 +254,12 @@ triple *copyTripleArray2C(array *a, size_t &N, GCPlacement placement) size_t n=checkArray(a); N=0; for(size_t i=0; i < n; i++) - N += checkArray(read<array*>(a,i)); - + N += checkArray(read<array*>(a,i)); + triple *A=(placement == NoGC) ? new triple [N] : new(placement) triple[N]; triple *p=A; - + for(size_t i=0; i < n; i++) { array *ai=read<array*>(a,i); size_t m=checkArray(ai); @@ -277,20 +277,20 @@ triple operator *(const array& t, const triple& v) 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 || + + 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+ @@ -303,14 +303,14 @@ template<class T> array *mult(array *a, array *b) { size_t n=checkArray(a); - + size_t nb=checkArray(b); size_t na0=n == 0 ? 0 : checkArray(read<array*>(a,0)); - if(na0 != nb) + if(na0 != nb) error(incommensurate); - + size_t nb0=nb == 0 ? 0 : checkArray(read<array*>(b,0)); - + array *c=new array(n); T *A,*B; @@ -329,10 +329,10 @@ array *mult(array *a, array *b) (*ci)[j]=sum; } } - + delete[] B; delete[] A; - + return c; } @@ -342,7 +342,7 @@ array *AtA(array *a) { size_t n=checkArray(a); size_t m=n == 0 ? 0 : checkArray(read<array*>(a,0)); - + array *c=new array(m); T *A; @@ -360,12 +360,12 @@ array *AtA(array *a) (*ci)[j]=sum; } } - + delete[] A; return c; } -double norm(double *a, size_t n) +double norm(double *a, size_t n) { if(n == 0) return 0.0; double M=fabs(a[0]); @@ -374,7 +374,7 @@ double norm(double *a, size_t n) return M; } -double norm(triple *a, size_t n) +double norm(triple *a, size_t n) { if(n == 0) return 0.0; double M=a[0].abs2(); @@ -393,7 +393,7 @@ void transpose(double *a, size_t n) double temp=a[ij]; a[ij]=a[ji]; a[ji]=temp; - } + } } } @@ -420,16 +420,16 @@ void inverse(double *M, size_t n) 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; @@ -437,10 +437,10 @@ void inverse(double *M, size_t n) } 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++) { @@ -465,7 +465,7 @@ void inverse(double *M, size_t n) } } ++(pivot[col]); - + // Interchange rows, if needed, to put the pivot element on the diagonal. double *acol=M+n*col; if(row != col) { @@ -476,8 +476,8 @@ void inverse(double *M, size_t n) acol[k]=temp; } } - - Row[i]=row; + + Row[i]=row; Col[i]=col; // Divide the pivot row by the pivot element. @@ -488,12 +488,12 @@ void inverse(double *M, size_t n) } real pivinv=1.0/denom; acol[col]=1.0; - for(size_t k=0; k < n; k++) + 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) { + if(k != col) { double *ak=M+n*k; real akcol=ak[col]; ak[col]=0.0; @@ -502,7 +502,7 @@ void inverse(double *M, size_t n) } } } - + // Unscramble the inverse matrix in view of the column interchanges. for(size_t k=n; k > 0;) { k--; @@ -593,7 +593,7 @@ Int LUdecompose(double *a, size_t n, size_t* index, bool warn=true) swap *= -1; vv[imax]=vv[j]; } - if(index) + if(index) index[j]=imax; if(j != n) { double denom=aj[j]; @@ -619,7 +619,7 @@ void dividebyzero(size_t i) buf << "Divide by zero"; error(buf); } - + void integeroverflow(size_t i) { ostringstream buf; @@ -704,7 +704,7 @@ void newAppendedArray(stack *Stack) for (Int index = n-1; index >= 0; index--) (*a)[index] = pop(Stack); - + copy(tail->begin(), tail->end(), back_inserter(*a)); {Stack->push<array*>(a); return;} @@ -739,18 +739,19 @@ void copyArray(stack *Stack) Int depth=vm::pop<Int>(Stack,Int_MAX); array * a=vm::pop<array *>(Stack); #line 669 "runarray.in" + if(a == 0) vm::error(dereferenceNullArray); if(depth < 0) error("cannot copy to a negative depth"); if(depth > typeDepth) depth=typeDepth; {Stack->push<array*>(a->copyToDepth(depth)); return;} } // Read an element from an array. Checks for initialization & bounds. -#line 676 "runarray.in" +#line 677 "runarray.in" void arrayRead(stack *Stack) { Int n=vm::pop<Int>(Stack); array * a=vm::pop<array *>(Stack); -#line 677 "runarray.in" +#line 678 "runarray.in" item& i=arrayRead(a,n); if (i.empty()) { ostringstream buf; @@ -761,37 +762,37 @@ void arrayRead(stack *Stack) } // Slice a substring from an array. -#line 688 "runarray.in" +#line 689 "runarray.in" void arraySliceRead(stack *Stack) { Int right=vm::pop<Int>(Stack); Int left=vm::pop<Int>(Stack); array * a=vm::pop<array *>(Stack); -#line 689 "runarray.in" +#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 696 "runarray.in" +#line 697 "runarray.in" void arraySliceReadToEnd(stack *Stack) { Int left=vm::pop<Int>(Stack); array * a=vm::pop<array *>(Stack); -#line 697 "runarray.in" +#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 704 "runarray.in" +#line 705 "runarray.in" void arrayArrayRead(stack *Stack) { Int n=vm::pop<Int>(Stack); array * a=vm::pop<array *>(Stack); -#line 705 "runarray.in" +#line 706 "runarray.in" item& i=arrayRead(a,n); if (i.empty()) i=new array(0); {Stack->push(i); return;} @@ -799,13 +800,13 @@ void arrayArrayRead(stack *Stack) // Write an element to an array. Increase size if necessary. // TODO: Add arrayWriteAndPop -#line 713 "runarray.in" +#line 714 "runarray.in" void arrayWrite(stack *Stack) { item value=vm::pop(Stack); Int n=vm::pop<Int>(Stack); array * a=vm::pop<array *>(Stack); -#line 714 "runarray.in" +#line 715 "runarray.in" size_t len=checkArray(a); bool cyclic=a->cyclic(); if(cyclic && len > 0) n=imod(n,len); @@ -819,27 +820,27 @@ void arrayWrite(stack *Stack) {Stack->push(value); return;} } -#line 728 "runarray.in" +#line 729 "runarray.in" void arraySliceWrite(stack *Stack) { array * src=vm::pop<array *>(Stack); Int right=vm::pop<Int>(Stack); Int left=vm::pop<Int>(Stack); array * dest=vm::pop<array *>(Stack); -#line 729 "runarray.in" +#line 730 "runarray.in" checkArray(src); checkArray(dest); dest->setSlice(left, right, src); {Stack->push<array*>(src); return;} } -#line 736 "runarray.in" +#line 737 "runarray.in" void arraySliceWriteToEnd(stack *Stack) { array * src=vm::pop<array *>(Stack); Int left=vm::pop<Int>(Stack); array * dest=vm::pop<array *>(Stack); -#line 737 "runarray.in" +#line 738 "runarray.in" checkArray(src); size_t len=checkArray(dest); dest->setSlice(left, (Int) len, src); @@ -847,20 +848,20 @@ void arraySliceWriteToEnd(stack *Stack) } // Returns the length of an array. -#line 745 "runarray.in" +#line 746 "runarray.in" void arrayLength(stack *Stack) { array * a=vm::pop<array *>(Stack); -#line 746 "runarray.in" +#line 747 "runarray.in" {Stack->push<Int>((Int) checkArray(a)); return;} } // Returns an array of integers representing the keys of the array. -#line 751 "runarray.in" +#line 752 "runarray.in" void arrayKeys(stack *Stack) { array * a=vm::pop<array *>(Stack); -#line 752 "runarray.in" +#line 753 "runarray.in" size_t size=checkArray(a); array *keys=new array(); @@ -874,33 +875,33 @@ void arrayKeys(stack *Stack) } // Return the cyclic flag for an array. -#line 766 "runarray.in" +#line 767 "runarray.in" void arrayCyclicFlag(stack *Stack) { array * a=vm::pop<array *>(Stack); -#line 767 "runarray.in" +#line 768 "runarray.in" checkArray(a); {Stack->push<bool>(a->cyclic()); return;} } -#line 772 "runarray.in" +#line 773 "runarray.in" void arraySetCyclicFlag(stack *Stack) { array * a=vm::pop<array *>(Stack); bool b=vm::pop<bool>(Stack); -#line 773 "runarray.in" +#line 774 "runarray.in" checkArray(a); a->cyclic(b); {Stack->push<bool>(b); return;} } // Check to see if an array element is initialized. -#line 780 "runarray.in" +#line 781 "runarray.in" void arrayInitializedHelper(stack *Stack) { array * a=vm::pop<array *>(Stack); Int n=vm::pop<Int>(Stack); -#line 781 "runarray.in" +#line 782 "runarray.in" size_t len=checkArray(a); bool cyclic=a->cyclic(); if(cyclic && len > 0) n=imod(n,len); @@ -910,62 +911,62 @@ void arrayInitializedHelper(stack *Stack) } // Returns the initialize method for an array. -#line 791 "runarray.in" +#line 792 "runarray.in" void arrayInitialized(stack *Stack) { array * a=vm::pop<array *>(Stack); -#line 792 "runarray.in" +#line 793 "runarray.in" {Stack->push<callable*>(new thunk(new bfunc(arrayInitializedHelper),a)); return;} } // The helper function for the cyclic method that sets the cyclic flag. -#line 797 "runarray.in" +#line 798 "runarray.in" void arrayCyclicHelper(stack *Stack) { array * a=vm::pop<array *>(Stack); bool b=vm::pop<bool>(Stack); -#line 798 "runarray.in" +#line 799 "runarray.in" checkArray(a); a->cyclic(b); } // Set the cyclic flag for an array. -#line 804 "runarray.in" +#line 805 "runarray.in" void arrayCyclic(stack *Stack) { array * a=vm::pop<array *>(Stack); -#line 805 "runarray.in" +#line 806 "runarray.in" {Stack->push<callable*>(new thunk(new bfunc(arrayCyclicHelper),a)); return;} } // The helper function for the push method that does the actual operation. -#line 810 "runarray.in" +#line 811 "runarray.in" void arrayPushHelper(stack *Stack) { array * a=vm::pop<array *>(Stack); item x=vm::pop(Stack); -#line 811 "runarray.in" +#line 812 "runarray.in" checkArray(a); a->push(x); {Stack->push(x); return;} } // Returns the push method for an array. -#line 818 "runarray.in" +#line 819 "runarray.in" void arrayPush(stack *Stack) { array * a=vm::pop<array *>(Stack); -#line 819 "runarray.in" +#line 820 "runarray.in" {Stack->push<callable*>(new thunk(new bfunc(arrayPushHelper),a)); return;} } // The helper function for the append method that appends b to a. -#line 824 "runarray.in" +#line 825 "runarray.in" void arrayAppendHelper(stack *Stack) { array * a=vm::pop<array *>(Stack); array * b=vm::pop<array *>(Stack); -#line 825 "runarray.in" +#line 826 "runarray.in" checkArray(a); size_t size=checkArray(b); for(size_t i=0; i < size; i++) @@ -973,85 +974,85 @@ void arrayAppendHelper(stack *Stack) } // Returns the append method for an array. -#line 833 "runarray.in" +#line 834 "runarray.in" void arrayAppend(stack *Stack) { array * a=vm::pop<array *>(Stack); -#line 834 "runarray.in" +#line 835 "runarray.in" {Stack->push<callable*>(new thunk(new bfunc(arrayAppendHelper),a)); return;} } // The helper function for the pop method. -#line 839 "runarray.in" +#line 840 "runarray.in" void arrayPopHelper(stack *Stack) { array * a=vm::pop<array *>(Stack); -#line 840 "runarray.in" +#line 841 "runarray.in" size_t asize=checkArray(a); - if(asize == 0) + if(asize == 0) error("cannot pop element from empty array"); {Stack->push(a->pop()); return;} } // Returns the pop method for an array. -#line 848 "runarray.in" +#line 849 "runarray.in" void arrayPop(stack *Stack) { array * a=vm::pop<array *>(Stack); -#line 849 "runarray.in" +#line 850 "runarray.in" {Stack->push<callable*>(new thunk(new bfunc(arrayPopHelper),a)); return;} } // The helper function for the insert method. -#line 854 "runarray.in" +#line 855 "runarray.in" void arrayInsertHelper(stack *Stack) { array * a=vm::pop<array *>(Stack); array * x=vm::pop<array *>(Stack); Int i=vm::pop<Int>(Stack); -#line 855 "runarray.in" +#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) + 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 865 "runarray.in" +#line 866 "runarray.in" void arrayInsert(stack *Stack) { array * a=vm::pop<array *>(Stack); -#line 866 "runarray.in" +#line 867 "runarray.in" {Stack->push<callable*>(new thunk(new bfunc(arrayInsertHelper),a)); return;} } // Returns the delete method for an array. -#line 871 "runarray.in" +#line 872 "runarray.in" void arrayDelete(stack *Stack) { array * a=vm::pop<array *>(Stack); -#line 872 "runarray.in" +#line 873 "runarray.in" {Stack->push<callable*>(new thunk(new bfunc(arrayDeleteHelper),a)); return;} } -#line 876 "runarray.in" +#line 877 "runarray.in" void arrayAlias(stack *Stack) { array * b=vm::pop<array *>(Stack); array * a=vm::pop<array *>(Stack); -#line 877 "runarray.in" +#line 878 "runarray.in" {Stack->push<bool>(a==b); return;} } // Return array formed by indexing array a with elements of integer array b -#line 882 "runarray.in" +#line 883 "runarray.in" void arrayIntArray(stack *Stack) { array * b=vm::pop<array *>(Stack); array * a=vm::pop<array *>(Stack); -#line 883 "runarray.in" +#line 884 "runarray.in" size_t asize=checkArray(a); size_t bsize=checkArray(b); array *r=new array(bsize); @@ -1069,13 +1070,13 @@ void arrayIntArray(stack *Stack) // 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 901 "runarray.in" +#line 902 "runarray.in" // Intarray* complement(Intarray *a, Int n); void gen_runarray32(stack *Stack) { Int n=vm::pop<Int>(Stack); Intarray * a=vm::pop<Intarray *>(Stack); -#line 902 "runarray.in" +#line 903 "runarray.in" size_t asize=checkArray(a); array *r=new array(0); bool *keep=new bool[n]; @@ -1086,18 +1087,18 @@ void gen_runarray32(stack *Stack) } for(Int i=0; i < n; i++) if(keep[i]) r->push(i); - + delete[] keep; {Stack->push<Intarray*>(r); return;} } // Generate the sequence {f(i) : i=0,1,...n-1} given a function f and integer n -#line 919 "runarray.in" +#line 920 "runarray.in" void arraySequence(stack *Stack) { Int n=vm::pop<Int>(Stack); callable * f=vm::pop<callable *>(Stack); -#line 920 "runarray.in" +#line 921 "runarray.in" if(n < 0) n=0; array *a=new array(n); for(Int i=0; i < n; ++i) { @@ -1109,12 +1110,12 @@ void arraySequence(stack *Stack) } // Return the array {0,1,...n-1} -#line 932 "runarray.in" +#line 933 "runarray.in" // Intarray* sequence(Int n); void gen_runarray34(stack *Stack) { Int n=vm::pop<Int>(Stack); -#line 933 "runarray.in" +#line 934 "runarray.in" if(n < 0) n=0; array *a=new array(n); for(Int i=0; i < n; ++i) { @@ -1124,12 +1125,12 @@ void gen_runarray34(stack *Stack) } // Apply a function to each element of an array -#line 943 "runarray.in" +#line 944 "runarray.in" void arrayFunction(stack *Stack) { array * a=vm::pop<array *>(Stack); callable * f=vm::pop<callable *>(Stack); -#line 944 "runarray.in" +#line 945 "runarray.in" size_t size=checkArray(a); array *b=new array(size); for(size_t i=0; i < size; ++i) { @@ -1140,13 +1141,13 @@ void arrayFunction(stack *Stack) {Stack->push<array*>(b); return;} } -#line 955 "runarray.in" +#line 956 "runarray.in" void arraySort(stack *Stack) { bool stable=vm::pop<bool>(Stack,true); callable * less=vm::pop<callable *>(Stack); array * a=vm::pop<array *>(Stack); -#line 956 "runarray.in" +#line 957 "runarray.in" array *c=copyArray(a); compareFunc=less; FuncStack=Stack; @@ -1155,13 +1156,13 @@ void arraySort(stack *Stack) {Stack->push<array*>(c); return;} } -#line 965 "runarray.in" +#line 966 "runarray.in" void arraySearch(stack *Stack) { callable * less=vm::pop<callable *>(Stack); item key=vm::pop(Stack); array * a=vm::pop<array *>(Stack); -#line 966 "runarray.in" +#line 967 "runarray.in" size_t size=a->size(); compareFunc=less; FuncStack=Stack; @@ -1169,7 +1170,7 @@ void arraySearch(stack *Stack) size_t u=size-1; if(!compareFunction(key,(*a)[u])) {Stack->push<Int>(Intcast(u)); return;} size_t l=0; - + while (l < u) { size_t i=(l+u)/2; if(compareFunction(key,(*a)[i])) u=i; @@ -1179,12 +1180,12 @@ void arraySearch(stack *Stack) {Stack->push<Int>(0); return;} } -#line 984 "runarray.in" +#line 985 "runarray.in" // bool all(boolarray *a); void gen_runarray38(stack *Stack) { boolarray * a=vm::pop<boolarray *>(Stack); -#line 985 "runarray.in" +#line 986 "runarray.in" size_t size=checkArray(a); bool c=true; for(size_t i=0; i < size; i++) @@ -1192,12 +1193,12 @@ void gen_runarray38(stack *Stack) {Stack->push<bool>(c); return;} } -#line 993 "runarray.in" +#line 994 "runarray.in" // boolarray* !(boolarray* a); void gen_runarray39(stack *Stack) { boolarray* a=vm::pop<boolarray*>(Stack); -#line 994 "runarray.in" +#line 995 "runarray.in" size_t size=checkArray(a); array *c=new array(size); for(size_t i=0; i < size; i++) @@ -1205,12 +1206,12 @@ void gen_runarray39(stack *Stack) {Stack->push<boolarray*>(c); return;} } -#line 1002 "runarray.in" +#line 1003 "runarray.in" // Int sum(boolarray *a); void gen_runarray40(stack *Stack) { boolarray * a=vm::pop<boolarray *>(Stack); -#line 1003 "runarray.in" +#line 1004 "runarray.in" size_t size=checkArray(a); Int sum=0; for(size_t i=0; i < size; i++) @@ -1218,11 +1219,11 @@ void gen_runarray40(stack *Stack) {Stack->push<Int>(sum); return;} } -#line 1011 "runarray.in" +#line 1012 "runarray.in" void arrayConcat(stack *Stack) { array * a=vm::pop<array *>(Stack); -#line 1012 "runarray.in" +#line 1013 "runarray.in" // a is an array of arrays to be concatenated together. // The signature is // T[] concat(... T[][] a); @@ -1249,28 +1250,30 @@ void arrayConcat(stack *Stack) {Stack->push<array*>(result); return;} } -#line 1039 "runarray.in" +#line 1040 "runarray.in" void array2Transpose(stack *Stack) { array * a=vm::pop<array *>(Stack); -#line 1040 "runarray.in" +#line 1041 "runarray.in" size_t asize=checkArray(a); array *c=new array(0); + size_t csize=0; for(size_t i=0; i < asize; i++) { size_t ip=i+1; array *ai=read<array*>(a,i); size_t aisize=checkArray(ai); - size_t csize=checkArray(c); - if(csize < aisize) { + if(c->size() < aisize) { c->resize(aisize); - for(size_t j=csize; j < aisize; j++) { - (*c)[j]=new array(ip); - } + for(size_t j=csize; j < aisize; j++) + (*c)[j]=new array(0); + csize=aisize; } 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]; + if(!(*ai)[j].empty()) { + array *cj=read<array*>(c,j); + if(checkArray(cj) < ip) cj->resize(ip); + (*cj)[i]=(*ai)[j]; + } } } {Stack->push<array*>(c); return;} @@ -1279,13 +1282,13 @@ void array2Transpose(stack *Stack) // 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 1067 "runarray.in" +// return the array obtained by putting a[i][j][k] into position perm{ijk}. +#line 1070 "runarray.in" void array3Transpose(stack *Stack) { array * perm=vm::pop<array *>(Stack); array * a=vm::pop<array *>(Stack); -#line 1068 "runarray.in" +#line 1071 "runarray.in" const size_t DIM=3; if(checkArray(perm) != DIM) { @@ -1293,10 +1296,10 @@ void array3Transpose(stack *Stack) 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; @@ -1307,13 +1310,13 @@ void array3Transpose(stack *Stack) } 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); @@ -1328,15 +1331,15 @@ void array3Transpose(stack *Stack) 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); @@ -1346,9 +1349,9 @@ void array3Transpose(stack *Stack) (*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]) { @@ -1360,22 +1363,22 @@ void array3Transpose(stack *Stack) } } } - - delete[] i; - delete[] size; + + delete[] i; + delete[] size; {Stack->push<array*>(c); return;} } // Find the index of the nth true value in a boolean array or -1 if not found. // If n is negative, search backwards. -#line 1152 "runarray.in" +#line 1155 "runarray.in" // Int find(boolarray *a, Int n=1); void gen_runarray44(stack *Stack) { Int n=vm::pop<Int>(Stack,1); boolarray * a=vm::pop<boolarray *>(Stack); -#line 1153 "runarray.in" +#line 1156 "runarray.in" size_t size=checkArray(a); Int j=-1; if(n > 0) @@ -1391,13 +1394,13 @@ void gen_runarray44(stack *Stack) {Stack->push<Int>(j); return;} } -// Find all indices of true values in a boolean array. -#line 1170 "runarray.in" +// Find all indices of true values in a boolean array. +#line 1173 "runarray.in" // Intarray* findall(boolarray *a); void gen_runarray45(stack *Stack) { boolarray * a=vm::pop<boolarray *>(Stack); -#line 1171 "runarray.in" +#line 1174 "runarray.in" size_t size=checkArray(a); array *b=new array(0); for(size_t i=0; i < size; i++) { @@ -1410,13 +1413,13 @@ void gen_runarray45(stack *Stack) // 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 1184 "runarray.in" +#line 1187 "runarray.in" void arrayConditional(stack *Stack) { array * c=vm::pop<array *>(Stack); array * b=vm::pop<array *>(Stack); array * a=vm::pop<array *>(Stack); -#line 1185 "runarray.in" +#line 1188 "runarray.in" size_t size=checkArray(a); array *r=new array(size); if(b && c) { @@ -1440,22 +1443,22 @@ void arrayConditional(stack *Stack) } // Return an n x n identity matrix. -#line 1209 "runarray.in" +#line 1212 "runarray.in" // realarray2* identity(Int n); void gen_runarray47(stack *Stack) { Int n=vm::pop<Int>(Stack); -#line 1210 "runarray.in" +#line 1213 "runarray.in" {Stack->push<realarray2*>(Identity(n)); return;} } // Return the inverse of an n x n matrix a using Gauss-Jordan elimination. -#line 1215 "runarray.in" +#line 1218 "runarray.in" // realarray2* inverse(realarray2 *a); void gen_runarray48(stack *Stack) { realarray2 * a=vm::pop<realarray2 *>(Stack); -#line 1216 "runarray.in" +#line 1219 "runarray.in" size_t n=checkArray(a); double *A; copyArray2C(A,a,true,0,NoGC); @@ -1468,33 +1471,33 @@ void gen_runarray48(stack *Stack) // 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 1229 "runarray.in" +#line 1232 "runarray.in" // realarray* solve(realarray2 *a, realarray *b, bool warn=true); void gen_runarray49(stack *Stack) { bool warn=vm::pop<bool>(Stack,true); realarray * b=vm::pop<realarray *>(Stack); realarray2 * a=vm::pop<realarray2 *>(Stack); -#line 1230 "runarray.in" +#line 1233 "runarray.in" size_t n=checkArray(a); - + if(n == 0) {Stack->push<realarray*>(new array(0)); return;} - + size_t m=checkArray(b); if(m != n) error(incommensurate); - + real *A; copyArray2C(A,a); size_t *index=new size_t[n]; - + if(LUdecompose(A,n,index,warn) == 0) {Stack->push<realarray*>(new array(0)); return;} array *x=new array(n); - + real *B; copyArrayC(B,b); - + for(size_t i=0; i < n; ++i) { size_t ip=index[i]; real sum=B[ip]; @@ -1504,7 +1507,7 @@ void gen_runarray49(stack *Stack) sum -= Ai[j]*B[j]; B[i]=sum; } - + for(size_t i=n; i > 0;) { --i; real sum=B[i]; @@ -1513,46 +1516,46 @@ void gen_runarray49(stack *Stack) sum -= Ai[j]*B[j]; B[i]=sum/Ai[i]; } - + for(size_t i=0; i < n; ++i) (*x)[i]=B[i]; delete[] index; delete[] B; delete[] A; - + {Stack->push<realarray*>(x); return;} } // Solve the linear equation ax=b by LU decomposition, returning the // solution x, where a is an n x n matrix and b is an n x m matrix. // If no solution exists, return an empty array. -#line 1282 "runarray.in" +#line 1285 "runarray.in" // realarray2* solve(realarray2 *a, realarray2 *b, bool warn=true); void gen_runarray50(stack *Stack) { bool warn=vm::pop<bool>(Stack,true); realarray2 * b=vm::pop<realarray2 *>(Stack); realarray2 * a=vm::pop<realarray2 *>(Stack); -#line 1283 "runarray.in" +#line 1286 "runarray.in" size_t n=checkArray(a); - + if(n == 0) {Stack->push<realarray2*>(new array(0)); return;} - + if(checkArray(b) != n) error(incommensurate); size_t m=checkArray(read<array*>(b,0)); - + real *A,*B; copyArray2C(A,a); copyArray2C(B,b,false); - + size_t *index=new size_t[n]; - + if(LUdecompose(A,n,index,warn) == 0) {Stack->push<realarray2*>(new array(0)); return;} array *x=new array(n); - + for(size_t i=0; i < n; ++i) { real *Ai=A+i*n; real *Bi=B+i*m; @@ -1566,7 +1569,7 @@ void gen_runarray50(stack *Stack) Bi[k]=sum; } } - + for(size_t i=n; i > 0;) { --i; real *Ai=A+i*n; @@ -1579,7 +1582,7 @@ void gen_runarray50(stack *Stack) Bi[k]=sum/Ai[i]; } } - + for(size_t i=0; i < n; ++i) { real *Bi=B+i*m; array *xi=new array(m); @@ -1587,42 +1590,42 @@ void gen_runarray50(stack *Stack) for(size_t j=0; j < m; ++j) (*xi)[j]=Bi[j]; } - + delete[] index; delete[] B; delete[] A; - + {Stack->push<realarray2*>(x); return;} } // Compute the determinant of an n x n matrix. -#line 1345 "runarray.in" +#line 1348 "runarray.in" // real determinant(realarray2 *a); void gen_runarray51(stack *Stack) { realarray2 * a=vm::pop<realarray2 *>(Stack); -#line 1346 "runarray.in" +#line 1349 "runarray.in" real *A; copyArray2C(A,a); size_t n=checkArray(a); - + real det=LUdecompose(A,n,NULL,false); size_t n1=n+1; for(size_t i=0; i < n; ++i) det *= A[i*n1]; - + delete[] A; - + {Stack->push<real>(det); return;} } -#line 1361 "runarray.in" +#line 1364 "runarray.in" // realarray* *(realarray2 *a, realarray *b); void gen_runarray52(stack *Stack) { realarray * b=vm::pop<realarray *>(Stack); realarray2 * a=vm::pop<realarray2 *>(Stack); -#line 1362 "runarray.in" +#line 1365 "runarray.in" size_t n=checkArray(a); size_t m=checkArray(b); array *c=new array(n); @@ -1640,13 +1643,13 @@ void gen_runarray52(stack *Stack) {Stack->push<realarray*>(c); return;} } -#line 1380 "runarray.in" +#line 1383 "runarray.in" // realarray* *(realarray *a, realarray2 *b); void gen_runarray53(stack *Stack) { realarray2 * b=vm::pop<realarray2 *>(Stack); realarray * a=vm::pop<realarray *>(Stack); -#line 1381 "runarray.in" +#line 1384 "runarray.in" size_t n=checkArray(a); if(n != checkArray(b)) error(incommensurate); real *A; @@ -1674,62 +1677,62 @@ void gen_runarray53(stack *Stack) {Stack->push<realarray*>(c); return;} } -#line 1409 "runarray.in" +#line 1412 "runarray.in" // Intarray2* *(Intarray2 *a, Intarray2 *b); void gen_runarray54(stack *Stack) { Intarray2 * b=vm::pop<Intarray2 *>(Stack); Intarray2 * a=vm::pop<Intarray2 *>(Stack); -#line 1410 "runarray.in" +#line 1413 "runarray.in" {Stack->push<Intarray2*>(mult<Int>(a,b)); return;} } -#line 1414 "runarray.in" +#line 1417 "runarray.in" // realarray2* *(realarray2 *a, realarray2 *b); void gen_runarray55(stack *Stack) { realarray2 * b=vm::pop<realarray2 *>(Stack); realarray2 * a=vm::pop<realarray2 *>(Stack); -#line 1415 "runarray.in" +#line 1418 "runarray.in" {Stack->push<realarray2*>(mult<real>(a,b)); return;} } -#line 1419 "runarray.in" +#line 1422 "runarray.in" // pairarray2* *(pairarray2 *a, pairarray2 *b); void gen_runarray56(stack *Stack) { pairarray2 * b=vm::pop<pairarray2 *>(Stack); pairarray2 * a=vm::pop<pairarray2 *>(Stack); -#line 1420 "runarray.in" +#line 1423 "runarray.in" {Stack->push<pairarray2*>(mult<pair>(a,b)); return;} } -#line 1424 "runarray.in" +#line 1427 "runarray.in" // triple *(realarray2 *t, triple v); void gen_runarray57(stack *Stack) { triple v=vm::pop<triple>(Stack); realarray2 * t=vm::pop<realarray2 *>(Stack); -#line 1425 "runarray.in" +#line 1428 "runarray.in" {Stack->push<triple>(*t*v); return;} } -#line 1429 "runarray.in" +#line 1432 "runarray.in" // realarray2* AtA(realarray2 *a); void gen_runarray58(stack *Stack) { realarray2 * a=vm::pop<realarray2 *>(Stack); -#line 1430 "runarray.in" +#line 1433 "runarray.in" {Stack->push<realarray2*>(AtA<real>(a)); return;} } -#line 1434 "runarray.in" +#line 1437 "runarray.in" // pair project(triple v, realarray2 *t); void gen_runarray59(stack *Stack) { realarray2 * t=vm::pop<realarray2 *>(Stack); triple v=vm::pop<triple>(Stack); -#line 1435 "runarray.in" +#line 1438 "runarray.in" size_t n=checkArray(t); if(n != 4) error(incommensurate); array *t0=read<array*>(t,0); @@ -1737,16 +1740,16 @@ void gen_runarray59(stack *Stack) array *t3=read<array*>(t,3); if(checkArray(t0) != 4 || checkArray(t1) != 4 || checkArray(t3) != 4) error(incommensurate); - + real x=v.getx(); real y=v.gety(); real z=v.getz(); - + real f=read<real>(t3,0)*x+read<real>(t3,1)*y+read<real>(t3,2)*z+ read<real>(t3,3); if(f == 0.0) dividebyzero(); f=1.0/f; - + {Stack->push<pair>(pair((read<real>(t0,0)*x+read<real>(t0,1)*y+read<real>(t0,2)*z+ read<real>(t0,3))*f, (read<real>(t1,0)*x+read<real>(t1,1)*y+read<real>(t1,2)*z+ @@ -1754,13 +1757,13 @@ void gen_runarray59(stack *Stack) } // Compute the dot product of vectors a and b. -#line 1460 "runarray.in" +#line 1463 "runarray.in" // real dot(realarray *a, realarray *b); void gen_runarray60(stack *Stack) { realarray * b=vm::pop<realarray *>(Stack); realarray * a=vm::pop<realarray *>(Stack); -#line 1461 "runarray.in" +#line 1464 "runarray.in" size_t n=checkArrays(a,b); real sum=0.0; for(size_t i=0; i < n; ++i) @@ -1769,13 +1772,13 @@ void gen_runarray60(stack *Stack) } // Compute the complex dot product of vectors a and b. -#line 1470 "runarray.in" +#line 1473 "runarray.in" // pair dot(pairarray *a, pairarray *b); void gen_runarray61(stack *Stack) { pairarray * b=vm::pop<pairarray *>(Stack); pairarray * a=vm::pop<pairarray *>(Stack); -#line 1471 "runarray.in" +#line 1474 "runarray.in" size_t n=checkArrays(a,b); pair sum=zero; for(size_t i=0; i < n; ++i) @@ -1790,7 +1793,7 @@ void gen_runarray61(stack *Stack) // [ a[2] b[2] c[2] ] // [ ... ] // [ c[n-1] a[n-1] b[n-1] ] -#line 1486 "runarray.in" +#line 1489 "runarray.in" // realarray* tridiagonal(realarray *a, realarray *b, realarray *c, realarray *f); void gen_runarray62(stack *Stack) { @@ -1798,26 +1801,26 @@ void gen_runarray62(stack *Stack) realarray * c=vm::pop<realarray *>(Stack); realarray * b=vm::pop<realarray *>(Stack); realarray * a=vm::pop<realarray *>(Stack); -#line 1487 "runarray.in" +#line 1490 "runarray.in" size_t n=checkArrays(a,b); checkEqual(n,checkArray(c)); checkEqual(n,checkArray(f)); - + array *up=new array(n); array& u=*up; if(n == 0) {Stack->push<realarray*>(up); return;} - + // Special case: zero Dirichlet boundary conditions if(read<real>(a,0) == 0.0 && read<real>(c,n-1) == 0.0) { real temp=read<real>(b,0); if(temp == 0.0) dividebyzero(); temp=1.0/temp; - + real *work=new real[n]; u[0]=read<real>(f,0)*temp; work[0]=-read<real>(c,0)*temp; - + for(size_t i=1; i < n; i++) { real temp=(read<real>(b,i)+read<real>(a,i)*work[i-1]); if(temp == 0.0) {delete[] work; dividebyzero();} @@ -1828,15 +1831,15 @@ void gen_runarray62(stack *Stack) for(size_t i=n-1; i >= 1; i--) u[i-1]=read<real>(u,i-1)+work[i-1]*read<real>(u,i); - + delete[] work; {Stack->push<realarray*>(up); return;} } - + real binv=read<real>(b,0); if(binv == 0.0) dividebyzero(); binv=1.0/binv; - + if(n == 1) {u[0]=read<real>(f,0)*binv; {Stack->push<realarray*>(up); return;}} if(n == 2) { real factor=(read<real>(b,0)*read<real>(b,1)- @@ -1850,10 +1853,10 @@ void gen_runarray62(stack *Stack) u[1]=temp; {Stack->push<realarray*>(up); return;} } - + real *gamma=new real[n-2]; real *delta=new real[n-2]; - + gamma[0]=read<real>(c,0)*binv; delta[0]=read<real>(a,0)*binv; u[0]=read<real>(f,0)*binv; @@ -1872,7 +1875,7 @@ void gen_runarray62(stack *Stack) 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; @@ -1884,18 +1887,18 @@ void gen_runarray62(stack *Stack) if(temp == 0.0) {delete[] gamma; delete[] delta; dividebyzero();} u[n-1]=temp=(fn-beta*read<real>(u,n-2))/temp; u[n-2]=read<real>(u,n-2)-dnm1*temp; - + for(size_t i=n-2; i >= 1; i--) u[i-1]=read<real>(u,i-1)-gamma[i-1]*read<real>(u,i)-delta[i-1]*temp; - + delete[] delta; delete[] gamma; - + {Stack->push<realarray*>(up); return;} } // Root solve by Newton-Raphson -#line 1584 "runarray.in" +#line 1587 "runarray.in" // real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x, bool verbose=false); void gen_runarray63(stack *Stack) { @@ -1904,22 +1907,22 @@ void gen_runarray63(stack *Stack) callableReal * fprime=vm::pop<callableReal *>(Stack); callableReal * f=vm::pop<callableReal *>(Stack); Int iterations=vm::pop<Int>(Stack,100); -#line 1586 "runarray.in" +#line 1589 "runarray.in" static const real fuzz=1000.0*DBL_EPSILON; Int i=0; size_t oldPrec=0; - if(verbose) + 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; @@ -1928,14 +1931,14 @@ void gen_runarray63(stack *Stack) 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; @@ -1950,7 +1953,7 @@ void gen_runarray63(stack *Stack) // Root solve by Newton-Raphson bisection // cf. routine rtsafe (Press et al., Numerical Recipes, 1991). -#line 1632 "runarray.in" +#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) { @@ -1960,22 +1963,22 @@ void gen_runarray64(stack *Stack) callableReal * fprime=vm::pop<callableReal *>(Stack); callableReal * f=vm::pop<callableReal *>(Stack); Int iterations=vm::pop<Int>(Stack,100); -#line 1634 "runarray.in" +#line 1637 "runarray.in" static const real fuzz=1000.0*DBL_EPSILON; size_t oldPrec=0; - if(verbose) + if(verbose) oldPrec=cout.precision(DBL_DIG); Stack->push(x1); f->call(Stack); real f1=pop<real>(Stack); if(f1 == 0.0) {Stack->push<real>(x1); return;} - + Stack->push(x2); f->call(Stack); real f2=pop<real>(Stack); if(f2 == 0.0) {Stack->push<real>(x2); return;} - + if((f1 > 0.0 && f2 > 0.0) || (f1 < 0.0 && f2 < 0.0)) { ostringstream buf; buf << "root not bracketed, f(x1)=" << f1 << ", f(x2)=" << f2 << endl; @@ -1989,7 +1992,7 @@ void gen_runarray64(stack *Stack) x1=x2; x2=temp; } - + if(verbose) cout << "midpoint: " << x << endl; @@ -1997,7 +2000,7 @@ void gen_runarray64(stack *Stack) Stack->push(x); f->call(Stack); real y=pop<real>(Stack); - + Stack->push(x); fprime->call(Stack); real dy=pop<real>(Stack); @@ -2021,11 +2024,11 @@ void gen_runarray64(stack *Stack) if(temp == x) {Stack->push<real>(x); return;} } if(fabs(dx) < fuzz*fabs(x)) {Stack->push<real>(x); return;} - + Stack->push(x); f->call(Stack); y=pop<real>(Stack); - + Stack->push(x); fprime->call(Stack); dy=pop<real>(Stack); @@ -2046,7 +2049,7 @@ void gen_runarray64(stack *Stack) // 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 1720 "runarray.in" +#line 1723 "runarray.in" // real _findroot(callableReal *f, real a, real b, real tolerance, real fa, real fb); void gen_runarray65(stack *Stack) { @@ -2056,13 +2059,13 @@ void gen_runarray65(stack *Stack) real b=vm::pop<real>(Stack); real a=vm::pop<real>(Stack); callableReal * f=vm::pop<callableReal *>(Stack); -#line 1722 "runarray.in" +#line 1725 "runarray.in" if(fa == 0.0) {Stack->push<real>(a); return;} if(fb == 0.0) {Stack->push<real>(b); return;} - + const char* oppsign="fa and fb must have opposite signs"; int sign; - + if(fa < 0.0) { if(fb < 0.0) error(oppsign); sign=1; @@ -2079,7 +2082,7 @@ void gen_runarray65(stack *Stack) while(b-a > tolerance) { t=(a+b)*0.5; - + Stack->push(t); f->call(Stack); ft=sign*pop<double>(Stack); @@ -2093,11 +2096,11 @@ void gen_runarray65(stack *Stack) 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) { @@ -2110,7 +2113,7 @@ void gen_runarray65(stack *Stack) } } } - + if(found) { if(ft > 0.0) { b=t; @@ -2119,7 +2122,7 @@ void gen_runarray65(stack *Stack) a=t; fa=ft; } - + t=root; // If the interpolated value is close to one edge of @@ -2128,7 +2131,7 @@ void gen_runarray65(stack *Stack) real margin=(b-a)*1.0e-3; if(t-a < margin) t=a+2.0*(t-a); else if(b-t < margin) t=b-2.0*(b-t); - + Stack->push(t); f->call(Stack); ft=sign*pop<double>(Stack); @@ -2136,7 +2139,7 @@ void gen_runarray65(stack *Stack) if(ft == 0.0) {Stack->push<real>(t); return;} } } - + if(ft > 0.0) { b=t; fb=ft; @@ -2148,7 +2151,7 @@ void gen_runarray65(stack *Stack) {Stack->push<real>(a-(b-a)/(fb-fa)*fa); return;} } -#line 1814 "runarray.in" +#line 1817 "runarray.in" // real simpson(callableReal *f, real a, real b, real acc=DBL_EPSILON, real dxmax=0); void gen_runarray66(stack *Stack) { @@ -2157,7 +2160,7 @@ void gen_runarray66(stack *Stack) real b=vm::pop<real>(Stack); real a=vm::pop<real>(Stack); callableReal * f=vm::pop<callableReal *>(Stack); -#line 1816 "runarray.in" +#line 1819 "runarray.in" real integral; if(dxmax <= 0) dxmax=fabs(b-a); callable *oldFunc=Func; @@ -2170,26 +2173,26 @@ void gen_runarray66(stack *Stack) } // Compute the fast Fourier transform of a pair array -#line 1829 "runarray.in" +#line 1832 "runarray.in" // pairarray* fft(pairarray *a, Int sign=1); void gen_runarray67(stack *Stack) { Int sign=vm::pop<Int>(Stack,1); pairarray * a=vm::pop<pairarray *>(Stack); -#line 1830 "runarray.in" +#line 1833 "runarray.in" #ifdef HAVE_LIBFFTW3 unsigned n=(unsigned) checkArray(a); array *c=new array(n); if(n) { Complex *f=utils::ComplexAlign(n); fftwpp::fft1d Forward(n,intcast(sign),f); - + for(size_t i=0; i < n; i++) { pair z=read<pair>(a,i); f[i]=Complex(z.getx(),z.gety()); } Forward.fft(f); - + for(size_t i=0; i < n; i++) { Complex z=f[i]; (*c)[i]=pair(z.real(),z.imag()); @@ -2206,13 +2209,13 @@ void gen_runarray67(stack *Stack) } // Compute the fast Fourier transform of a 2D pair array -#line 1860 "runarray.in" +#line 1863 "runarray.in" // pairarray2* fft(pairarray2 *a, Int sign=1); void gen_runarray68(stack *Stack) { Int sign=vm::pop<Int>(Stack,1); pairarray2 * a=vm::pop<pairarray2 *>(Stack); -#line 1861 "runarray.in" +#line 1864 "runarray.in" #ifdef HAVE_LIBFFTW3 size_t n=checkArray(a); size_t m=n == 0 ? 0 : checkArray(read<array*>(a,0)); @@ -2257,13 +2260,13 @@ void gen_runarray68(stack *Stack) } // Compute the fast Fourier transform of a 3D pair array -#line 1906 "runarray.in" +#line 1909 "runarray.in" // pairarray3* fft(pairarray3 *a, Int sign=1); void gen_runarray69(stack *Stack) { Int sign=vm::pop<Int>(Stack,1); pairarray3 * a=vm::pop<pairarray3 *>(Stack); -#line 1907 "runarray.in" +#line 1910 "runarray.in" #ifdef HAVE_LIBFFTW3 size_t n=checkArray(a); array *a0=read<array*>(a,0); @@ -2320,26 +2323,26 @@ void gen_runarray69(stack *Stack) {Stack->push<pairarray3*>(c); return;} } -#line 1964 "runarray.in" +#line 1967 "runarray.in" // Intarray2* triangulate(pairarray *z); void gen_runarray70(stack *Stack) { pairarray * z=vm::pop<pairarray *>(Stack); -#line 1965 "runarray.in" +#line 1968 "runarray.in" size_t nv=checkArray(z); // Call robust version of Gilles Dumoulin's port of Paul Bourke's // triangulation code. XYZ *pxyz=new XYZ[nv+3]; ITRIANGLE *V=new ITRIANGLE[4*nv]; - + for(size_t i=0; i < nv; ++i) { pair w=read<pair>(z,i); pxyz[i].p[0]=w.getx(); pxyz[i].p[1]=w.gety(); pxyz[i].i=(Int) i; } - + Int ntri; Triangulate((Int) nv,pxyz,V,ntri,true,false); @@ -2353,18 +2356,18 @@ void gen_runarray70(stack *Stack) (*ti)[1]=pxyz[Vi->p2].i; (*ti)[2]=pxyz[Vi->p3].i; } - + delete[] V; delete[] pxyz; {Stack->push<Intarray2*>(t); return;} } -#line 1999 "runarray.in" +#line 2002 "runarray.in" // real norm(realarray *a); void gen_runarray71(stack *Stack) { realarray * a=vm::pop<realarray *>(Stack); -#line 2000 "runarray.in" +#line 2003 "runarray.in" size_t n=checkArray(a); real M=0.0; for(size_t i=0; i < n; ++i) { @@ -2374,12 +2377,12 @@ void gen_runarray71(stack *Stack) {Stack->push<real>(M); return;} } -#line 2010 "runarray.in" +#line 2013 "runarray.in" // real norm(realarray2 *a); void gen_runarray72(stack *Stack) { realarray2 * a=vm::pop<realarray2 *>(Stack); -#line 2011 "runarray.in" +#line 2014 "runarray.in" size_t n=checkArray(a); real M=0.0; for(size_t i=0; i < n; ++i) { @@ -2393,12 +2396,12 @@ void gen_runarray72(stack *Stack) {Stack->push<real>(M); return;} } -#line 2025 "runarray.in" +#line 2028 "runarray.in" // real norm(triplearray2 *a); void gen_runarray73(stack *Stack) { triplearray2 * a=vm::pop<triplearray2 *>(Stack); -#line 2026 "runarray.in" +#line 2029 "runarray.in" size_t n=checkArray(a); real M=0.0; for(size_t i=0; i < n; ++i) { @@ -2412,21 +2415,21 @@ void gen_runarray73(stack *Stack) {Stack->push<real>(sqrt(M)); return;} } -#line 2040 "runarray.in" +#line 2043 "runarray.in" // real change2(triplearray2 *a); void gen_runarray74(stack *Stack) { triplearray2 * a=vm::pop<triplearray2 *>(Stack); -#line 2041 "runarray.in" +#line 2044 "runarray.in" size_t n=checkArray(a); if(n == 0) {Stack->push<real>(0.0); return;} - + vm::array *a0=vm::read<vm::array*>(a,0); size_t m=checkArray(a0); if(m == 0) {Stack->push<real>(0.0); return;} triple a00=vm::read<triple>(a0,0); real M=0.0; - + for(size_t i=0; i < n; ++i) { vm::array *ai=vm::read<vm::array*>(a,i); size_t m=checkArray(ai); @@ -2438,13 +2441,13 @@ void gen_runarray74(stack *Stack) {Stack->push<real>(M); return;} } -#line 2062 "runarray.in" +#line 2065 "runarray.in" // triple minbezier(triplearray2 *P, triple b); void gen_runarray75(stack *Stack) { triple b=vm::pop<triple>(Stack); triplearray2 * P=vm::pop<triplearray2 *>(Stack); -#line 2063 "runarray.in" +#line 2066 "runarray.in" size_t N; real *A=copyTripleArray2Components(P,N); bound_double *B=bounddouble(N); @@ -2455,13 +2458,13 @@ void gen_runarray75(stack *Stack) {Stack->push<triple>(b); return;} } -#line 2074 "runarray.in" +#line 2077 "runarray.in" // triple maxbezier(triplearray2 *P, triple b); void gen_runarray76(stack *Stack) { triple b=vm::pop<triple>(Stack); triplearray2 * P=vm::pop<triplearray2 *>(Stack); -#line 2075 "runarray.in" +#line 2078 "runarray.in" size_t N; real *A=copyTripleArray2Components(P,N); bound_double *B=bounddouble(N); @@ -2472,13 +2475,13 @@ void gen_runarray76(stack *Stack) {Stack->push<triple>(b); return;} } -#line 2086 "runarray.in" +#line 2089 "runarray.in" // pair minratio(triplearray2 *P, pair b); void gen_runarray77(stack *Stack) { pair b=vm::pop<pair>(Stack); triplearray2 * P=vm::pop<triplearray2 *>(Stack); -#line 2087 "runarray.in" +#line 2090 "runarray.in" size_t N; triple *A=copyTripleArray2C(P,N); real fuzz=Fuzz*norm(A,N); @@ -2489,13 +2492,13 @@ void gen_runarray77(stack *Stack) {Stack->push<pair>(b); return;} } -#line 2098 "runarray.in" +#line 2101 "runarray.in" // pair maxratio(triplearray2 *P, pair b); void gen_runarray78(stack *Stack) { pair b=vm::pop<pair>(Stack); triplearray2 * P=vm::pop<triplearray2 *>(Stack); -#line 2099 "runarray.in" +#line 2102 "runarray.in" size_t N; triple *A=copyTripleArray2C(P,N); bound_triple *B=boundtriple(N); @@ -2506,35 +2509,35 @@ void gen_runarray78(stack *Stack) {Stack->push<pair>(b); return;} } -#line 2110 "runarray.in" +#line 2113 "runarray.in" // realarray* _projection(); void gen_runarray79(stack *Stack) { -#line 2111 "runarray.in" +#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 @@ -2559,153 +2562,153 @@ void gen_runarray_venv(venv &ve) REGISTER_BLTIN(run::copyArrayValue,"copyArrayValue"); #line 664 "runarray.in" REGISTER_BLTIN(run::copyArray,"copyArray"); -#line 675 "runarray.in" +#line 676 "runarray.in" REGISTER_BLTIN(run::arrayRead,"arrayRead"); -#line 687 "runarray.in" +#line 688 "runarray.in" REGISTER_BLTIN(run::arraySliceRead,"arraySliceRead"); -#line 694 "runarray.in" +#line 695 "runarray.in" REGISTER_BLTIN(run::arraySliceReadToEnd,"arraySliceReadToEnd"); -#line 702 "runarray.in" +#line 703 "runarray.in" REGISTER_BLTIN(run::arrayArrayRead,"arrayArrayRead"); -#line 711 "runarray.in" +#line 712 "runarray.in" REGISTER_BLTIN(run::arrayWrite,"arrayWrite"); -#line 728 "runarray.in" +#line 729 "runarray.in" REGISTER_BLTIN(run::arraySliceWrite,"arraySliceWrite"); -#line 736 "runarray.in" +#line 737 "runarray.in" REGISTER_BLTIN(run::arraySliceWriteToEnd,"arraySliceWriteToEnd"); -#line 744 "runarray.in" +#line 745 "runarray.in" REGISTER_BLTIN(run::arrayLength,"arrayLength"); -#line 750 "runarray.in" +#line 751 "runarray.in" REGISTER_BLTIN(run::arrayKeys,"arrayKeys"); -#line 765 "runarray.in" +#line 766 "runarray.in" REGISTER_BLTIN(run::arrayCyclicFlag,"arrayCyclicFlag"); -#line 772 "runarray.in" +#line 773 "runarray.in" REGISTER_BLTIN(run::arraySetCyclicFlag,"arraySetCyclicFlag"); -#line 779 "runarray.in" +#line 780 "runarray.in" REGISTER_BLTIN(run::arrayInitializedHelper,"arrayInitializedHelper"); -#line 790 "runarray.in" +#line 791 "runarray.in" REGISTER_BLTIN(run::arrayInitialized,"arrayInitialized"); -#line 796 "runarray.in" +#line 797 "runarray.in" REGISTER_BLTIN(run::arrayCyclicHelper,"arrayCyclicHelper"); -#line 803 "runarray.in" +#line 804 "runarray.in" REGISTER_BLTIN(run::arrayCyclic,"arrayCyclic"); -#line 809 "runarray.in" +#line 810 "runarray.in" REGISTER_BLTIN(run::arrayPushHelper,"arrayPushHelper"); -#line 817 "runarray.in" +#line 818 "runarray.in" REGISTER_BLTIN(run::arrayPush,"arrayPush"); -#line 823 "runarray.in" +#line 824 "runarray.in" REGISTER_BLTIN(run::arrayAppendHelper,"arrayAppendHelper"); -#line 832 "runarray.in" +#line 833 "runarray.in" REGISTER_BLTIN(run::arrayAppend,"arrayAppend"); -#line 838 "runarray.in" +#line 839 "runarray.in" REGISTER_BLTIN(run::arrayPopHelper,"arrayPopHelper"); -#line 847 "runarray.in" +#line 848 "runarray.in" REGISTER_BLTIN(run::arrayPop,"arrayPop"); -#line 853 "runarray.in" +#line 854 "runarray.in" REGISTER_BLTIN(run::arrayInsertHelper,"arrayInsertHelper"); -#line 864 "runarray.in" +#line 865 "runarray.in" REGISTER_BLTIN(run::arrayInsert,"arrayInsert"); -#line 870 "runarray.in" +#line 871 "runarray.in" REGISTER_BLTIN(run::arrayDelete,"arrayDelete"); -#line 876 "runarray.in" +#line 877 "runarray.in" REGISTER_BLTIN(run::arrayAlias,"arrayAlias"); -#line 881 "runarray.in" +#line 882 "runarray.in" REGISTER_BLTIN(run::arrayIntArray,"arrayIntArray"); -#line 899 "runarray.in" +#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 918 "runarray.in" +#line 919 "runarray.in" REGISTER_BLTIN(run::arraySequence,"arraySequence"); -#line 931 "runarray.in" +#line 932 "runarray.in" addFunc(ve, run::gen_runarray34, IntArray(), SYM(sequence), formal(primInt(), SYM(n), false, false)); -#line 942 "runarray.in" +#line 943 "runarray.in" REGISTER_BLTIN(run::arrayFunction,"arrayFunction"); -#line 955 "runarray.in" +#line 956 "runarray.in" REGISTER_BLTIN(run::arraySort,"arraySort"); -#line 965 "runarray.in" +#line 966 "runarray.in" REGISTER_BLTIN(run::arraySearch,"arraySearch"); -#line 984 "runarray.in" +#line 985 "runarray.in" addFunc(ve, run::gen_runarray38, primBoolean(), SYM(all), formal(booleanArray(), SYM(a), false, false)); -#line 993 "runarray.in" +#line 994 "runarray.in" addFunc(ve, run::gen_runarray39, booleanArray(), SYM_LOGNOT, formal(booleanArray(), SYM(a), false, false)); -#line 1002 "runarray.in" +#line 1003 "runarray.in" addFunc(ve, run::gen_runarray40, primInt(), SYM(sum), formal(booleanArray(), SYM(a), false, false)); -#line 1011 "runarray.in" +#line 1012 "runarray.in" REGISTER_BLTIN(run::arrayConcat,"arrayConcat"); -#line 1039 "runarray.in" +#line 1040 "runarray.in" REGISTER_BLTIN(run::array2Transpose,"array2Transpose"); -#line 1063 "runarray.in" +#line 1066 "runarray.in" REGISTER_BLTIN(run::array3Transpose,"array3Transpose"); -#line 1150 "runarray.in" +#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 1169 "runarray.in" +#line 1172 "runarray.in" addFunc(ve, run::gen_runarray45, IntArray(), SYM(findall), formal(booleanArray(), SYM(a), false, false)); -#line 1182 "runarray.in" +#line 1185 "runarray.in" REGISTER_BLTIN(run::arrayConditional,"arrayConditional"); -#line 1208 "runarray.in" +#line 1211 "runarray.in" addFunc(ve, run::gen_runarray47, realArray2(), SYM(identity), formal(primInt(), SYM(n), false, false)); -#line 1214 "runarray.in" +#line 1217 "runarray.in" addFunc(ve, run::gen_runarray48, realArray2(), SYM(inverse), formal(realArray2(), SYM(a), false, false)); -#line 1226 "runarray.in" +#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 1279 "runarray.in" +#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 1344 "runarray.in" +#line 1347 "runarray.in" addFunc(ve, run::gen_runarray51, primReal(), SYM(determinant), formal(realArray2(), SYM(a), false, false)); -#line 1361 "runarray.in" +#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 1380 "runarray.in" +#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 1409 "runarray.in" +#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 1414 "runarray.in" +#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 1419 "runarray.in" +#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 1424 "runarray.in" +#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 1429 "runarray.in" +#line 1432 "runarray.in" addFunc(ve, run::gen_runarray58, realArray2(), SYM(AtA), formal(realArray2(), SYM(a), false, false)); -#line 1434 "runarray.in" +#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 1459 "runarray.in" +#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 1469 "runarray.in" +#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 1479 "runarray.in" +#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 1583 "runarray.in" +#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 1630 "runarray.in" +#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 1712 "runarray.in" +#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 1814 "runarray.in" +#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 1828 "runarray.in" +#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 1859 "runarray.in" +#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 1905 "runarray.in" +#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 1964 "runarray.in" +#line 1967 "runarray.in" addFunc(ve, run::gen_runarray70, IntArray2(), SYM(triangulate), formal(pairArray(), SYM(z), false, false)); -#line 1999 "runarray.in" +#line 2002 "runarray.in" addFunc(ve, run::gen_runarray71, primReal(), SYM(norm), formal(realArray(), SYM(a), false, false)); -#line 2010 "runarray.in" +#line 2013 "runarray.in" addFunc(ve, run::gen_runarray72, primReal(), SYM(norm), formal(realArray2(), SYM(a), false, false)); -#line 2025 "runarray.in" +#line 2028 "runarray.in" addFunc(ve, run::gen_runarray73, primReal(), SYM(norm), formal(tripleArray2(), SYM(a), false, false)); -#line 2040 "runarray.in" +#line 2043 "runarray.in" addFunc(ve, run::gen_runarray74, primReal(), SYM(change2), formal(tripleArray2(), SYM(a), false, false)); -#line 2062 "runarray.in" +#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 2074 "runarray.in" +#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 2086 "runarray.in" +#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 2098 "runarray.in" +#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 2110 "runarray.in" +#line 2113 "runarray.in" addFunc(ve, run::gen_runarray79, realArray(), SYM(_projection)); } |