summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/runtime.in
diff options
context:
space:
mode:
authorKarl Berry <karl@freefriends.org>2009-05-16 00:19:13 +0000
committerKarl Berry <karl@freefriends.org>2009-05-16 00:19:13 +0000
commitbab45528d65eaafe68a705dbb2a57075c7b7cbd8 (patch)
tree10b4ae2b5195c8dede153ab89359ec00f55f325f /Build/source/utils/asymptote/runtime.in
parent8643d90372e9c31e0f461c15c596b60a545bd7d3 (diff)
asymptote 1.72 sources (not integrated into build yet)
git-svn-id: svn://tug.org/texlive/trunk@13110 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/utils/asymptote/runtime.in')
-rw-r--r--Build/source/utils/asymptote/runtime.in5818
1 files changed, 5818 insertions, 0 deletions
diff --git a/Build/source/utils/asymptote/runtime.in b/Build/source/utils/asymptote/runtime.in
new file mode 100644
index 00000000000..decee0727b1
--- /dev/null
+++ b/Build/source/utils/asymptote/runtime.in
@@ -0,0 +1,5818 @@
+/*****
+ * runtime.in
+ * Tom Prince 2005/4/15
+ *
+ * Generate the runtime functions used by the vm::stack machine.
+ *
+ *****/
+
+/* Autogenerated routines are specified like this (separated by a formfeed):
+ type asyname:cname(cparams)
+ {
+ C code
+ }
+
+*/
+
+// Use Void f() instead of void f() to force an explicit Stack argument.
+
+
+void => primVoid()
+Void => primVoid()
+Int => primInt()
+bool => primBoolean()
+double => primReal()
+real => primReal()
+string* => primString()
+string => primString()
+pen => primPen()
+pair => primPair()
+triple => primTriple()
+path => primPath()
+path3 => primPath3()
+guide* => primGuide()
+cycleToken => primCycleToken()
+tensionSpecifier => primTensionSpecifier()
+curlSpecifier => primCurlSpecifier()
+file* => primFile()
+picture* => primPicture()
+transform => primTransform()
+callable* => voidFunction()
+callableBp* => breakpointFunction()
+callableReal* => realRealFunction()
+callableTransform* => transformFunction()
+runnable* => primCode()
+boolarray* => boolArray()
+Intarray* => IntArray()
+Intarray2* => IntArray2()
+realarray* => realArray()
+realarray2* => realArray2()
+pairarray* => pairArray()
+pairarray2* => pairArray2()
+triplearray* => tripleArray()
+triplearray2* => tripleArray2()
+patharray* => pathArray()
+patharray2* => pathArray2()
+guidearray* => guideArray()
+transformarray* => transformArray()
+penarray* => penArray()
+penarray2* => penArray2()
+stringarray* => stringArray()
+stringarray2* => stringArray2()
+
+#include <cfloat>
+#include <time.h>
+#include <sys/times.h>
+#include <locale.h>
+
+#include "angle.h"
+#include "pair.h"
+#include "triple.h"
+#include "transform.h"
+#include "path.h"
+#include "path3.h"
+#include "pen.h"
+#include "guide.h"
+#include "picture.h"
+#include "drawpath.h"
+#include "drawpath3.h"
+#include "drawsurface.h"
+#include "drawfill.h"
+#include "drawclipbegin.h"
+#include "drawclipend.h"
+#include "drawlabel.h"
+#include "drawverbatim.h"
+#include "drawgsave.h"
+#include "drawgrestore.h"
+#include "drawlayer.h"
+#include "drawimage.h"
+#include "drawgroup.h"
+#include "fileio.h"
+#include "genv.h"
+#include "builtin.h"
+#include "texfile.h"
+#include "pipestream.h"
+#include "parser.h"
+#include "stack.h"
+#include "util.h"
+#include "locate.h"
+#include "mathop.h"
+#include "callable.h"
+#include "stm.h"
+#include "lexical.h"
+#include "process.h"
+#include "arrayop.h"
+#include "predicates.h"
+#include "Delaunay.h"
+
+#ifdef HAVE_LIBFFTW3
+#include "fftw++.h"
+#endif
+
+#if defined(HAVE_LIBREADLINE) && defined(HAVE_LIBCURSES)
+#include <readline/readline.h>
+#include <readline/history.h>
+#endif
+
+#if defined(USEGC) && defined(GC_DEBUG) && defined(GC_BACKTRACE)
+ extern "C" {
+ void *GC_generate_random_valid_address(void);
+ void GC_debug_print_heap_obj_proc(void *);
+ }
+#endif
+
+using namespace vm;
+using namespace camp;
+using namespace settings;
+
+namespace run {
+using camp::pair;
+using vm::array;
+using vm::frame;
+using vm::stack;
+using camp::transform;
+using absyntax::runnable;
+
+typedef double real;
+
+#define CURRENTPEN processData().currentpen
+
+typedef array boolarray;
+typedef array Intarray;
+typedef array Intarray2;
+typedef array realarray;
+typedef array realarray2;
+typedef array pairarray;
+typedef array pairarray2;
+typedef array triplearray;
+typedef array triplearray2;
+typedef array patharray;
+typedef array patharray2;
+typedef array guidearray;
+typedef array transformarray;
+typedef array penarray;
+typedef array penarray2;
+typedef array stringarray;
+typedef array stringarray2;
+
+typedef callable callableBp;
+typedef callable callableReal;
+typedef callable callableTransform;
+}
+
+using vm::array;
+using types::function;
+
+#define PRIMITIVE(name,Name,asyName) using types::prim##Name;
+#include <primitives.h>
+#undef PRIMITIVE
+
+using types::boolArray;
+using types::IntArray;
+using types::IntArray2;
+using types::realArray;
+using types::realArray2;
+using types::pairArray;
+using types::pairArray2;
+using types::tripleArray;
+using types::tripleArray2;
+using types::pathArray;
+using types::pathArray2;
+using types::guideArray;
+using types::transformArray;
+using types::penArray;
+using types::penArray2;
+using types::stringArray;
+using types::stringArray2;
+using types::formal;
+
+function *voidFunction()
+{
+ return new function(primVoid());
+}
+
+function *breakpointFunction()
+{
+ return new function(primString(),primString(),primInt(),primInt(),
+ primCode());
+}
+
+function *realRealFunction()
+{
+ return new function(primReal(),primReal());
+}
+
+function *transformFunction()
+{
+ return new function(primTransform());
+}
+
+function *realTripleFunction()
+{
+ return new function(primReal(),primTriple());
+}
+
+const size_t camp::ColorComponents[]={0,0,1,3,4,0};
+
+namespace vm {
+}
+
+namespace run {
+
+const char *invalidargument="invalid argument";
+const char *arrayempty="cannot take min or max of empty array";
+
+// Return the factorial of a non-negative integer using a lookup table.
+Int factorial(Int n)
+{
+ static Int *table;
+ static Int size=0;
+ if(size == 0) {
+ Int f=1;
+ size=2;
+ while(f <= Int_MAX/size)
+ f *= (size++);
+ table=new Int[size];
+ table[0]=f=1;
+ for(Int i=1; i < size; ++i) {
+ f *= i;
+ table[i]=f;
+ }
+ }
+ if(n >= size) integeroverflow(0);
+ return table[n];
+}
+
+static inline Int Round(double x)
+{
+ return Int(x+((x >= 0) ? 0.5 : -0.5));
+}
+
+inline Int sgn(double x)
+{
+ return (x > 0.0 ? 1 : (x < 0.0 ? -1 : 0));
+}
+
+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;
+ }
+}
+
+array *nop(array *a)
+{
+ return a;
+}
+
+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;
+}
+
+array *copyArray(array *a)
+{
+ size_t size=checkArray(a);
+ array *c=new array(size);
+ for(size_t i=0; i < size; i++)
+ (*c)[i]=(*a)[i];
+ return c;
+}
+
+inline size_t checkdimension(array *a, size_t dim)
+{
+ size_t size=checkArray(a);
+ if(dim && size != dim) {
+ ostringstream buf;
+ buf << "array of length " << dim << " expected" << endl;
+ error(buf);
+ }
+ return size;
+}
+
+double *copyArrayC(array *a, size_t dim=0)
+{
+ size_t size=checkdimension(a,dim);
+ double *c=new double[size];
+ for(size_t i=0; i < size; i++)
+ c[i]=read<double>(a,i);
+ return c;
+}
+
+triple *copyTripleArrayC(array *a, size_t dim=0)
+{
+ size_t size=checkdimension(a,dim);
+ triple *c=new triple[size];
+ for(size_t i=0; i < size; i++)
+ c[i]=read<triple>(a,i);
+ return c;
+}
+
+array *copyArray2(array *a)
+{
+ size_t size=checkArray(a);
+ array *c=new array(size);
+ for(size_t i=0; i < size; i++) {
+ array *ai=read<array*>(a,i);
+ size_t aisize=checkArray(ai);
+ array *ci=new array(aisize);
+ (*c)[i]=ci;
+ for(size_t j=0; j < aisize; j++)
+ (*ci)[j]=(*ai)[j];
+ }
+ return c;
+}
+
+array *copyArray3(array *a)
+{
+ size_t size=checkArray(a);
+ array *c=new array(size);
+ for(size_t i=0; i < size; i++) {
+ array *ai=read<array*>(a,i);
+ size_t aisize=checkArray(ai);
+ array *ci=new array(aisize);
+ (*c)[i]=ci;
+ for(size_t j=0; j < aisize; j++) {
+ array *aij=read<array*>(ai,j);
+ size_t aijsize=checkArray(aij);
+ array *cij=new array(aijsize);
+ (*ci)[j]=cij;
+ for(size_t k=0; k < aijsize; k++)
+ (*cij)[k]=(*aij)[k];
+ }
+ }
+ return c;
+}
+
+double *copyArray2C(array *a, bool square=true, size_t dim2=0)
+{
+ size_t n=checkArray(a);
+ size_t m=(square || n == 0) ? n : checkArray(read<array*>(a,0));
+ if(n > 0 && dim2 && m != dim2) {
+ ostringstream buf;
+ buf << "second matrix dimension must be " << dim2 << endl;
+ error(buf);
+ }
+
+ double *c=new double[n*m];
+ for(size_t i=0; i < n; i++) {
+ array *ai=read<array*>(a,i);
+ size_t aisize=checkArray(ai);
+ if(aisize == m) {
+ double *ci=c+i*m;
+ for(size_t j=0; j < m; j++)
+ ci[j]=read<double>(ai,j);
+ } else
+ error(square ? "matrix must be square" : "matrix must be rectangular");
+ }
+ return c;
+}
+
+static const char *incommensurate="Incommensurate matrices";
+static const char *singular="Singular matrix";
+static size_t *pivot,*Row,*Col;
+
+triple operator *(const array& t, const triple& v)
+{
+ size_t n=checkArray(&t);
+ if(n != 4) error(incommensurate);
+ array *t0=read<array*>(t,0);
+ array *t1=read<array*>(t,1);
+ array *t2=read<array*>(t,2);
+ array *t3=read<array*>(t,3);
+
+ if(checkArray(t0) != 4 || checkArray(t1) != 4 ||
+ checkArray(t2) != 4 || checkArray(t3) != 4)
+ error(incommensurate);
+
+ double x=v.getx();
+ double y=v.gety();
+ double z=v.getz();
+
+ double f=read<real>(t3,0)*x+read<real>(t3,1)*y+read<real>(t3,2)*z+
+ read<real>(t3,3);
+ if(f == 0.0) run::dividebyzero();
+ f=1.0/f;
+
+ return triple((read<real>(t0,0)*x+read<real>(t0,1)*y+read<real>(t0,2)*z+
+ read<real>(t0,3))*f,
+ (read<real>(t1,0)*x+read<real>(t1,1)*y+read<real>(t1,2)*z+
+ read<real>(t1,3))*f,
+ (read<real>(t2,0)*x+read<real>(t2,1)*y+read<real>(t2,2)*z+
+ read<real>(t2,3))*f);
+}
+
+triple multshiftless(const array& t, const triple& v)
+{
+ size_t n=checkArray(&t);
+ if(n != 4) error(incommensurate);
+ array *t0=read<array*>(t,0);
+ array *t1=read<array*>(t,1);
+ array *t2=read<array*>(t,2);
+ array *t3=read<array*>(t,3);
+
+ if(checkArray(t0) != 4 || checkArray(t1) != 4 ||
+ checkArray(t2) != 4 || checkArray(t3) != 4)
+ error(incommensurate);
+
+ double x=v.getx();
+ double y=v.gety();
+ double z=v.getz();
+
+ double f=read<real>(t3,0)*x+read<real>(t3,1)*y+read<real>(t3,2)*z+
+ read<real>(t3,3);
+ if(f == 0.0) run::dividebyzero();
+ f=1.0/f;
+
+ return triple((read<real>(t0,0)*x+read<real>(t0,1)*y+read<real>(t0,2)*z)*f,
+ (read<real>(t1,0)*x+read<real>(t1,1)*y+read<real>(t1,2)*z)*f,
+ (read<real>(t2,0)*x+read<real>(t2,1)*y+read<real>(t2,2)*z)*f);
+}
+
+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;
+}
+
+void writestring(stack *s)
+{
+ callable *suffix=pop<callable *>(s,NULL);
+ string S=pop<string>(s);
+ vm::item it=pop(s);
+ bool defaultfile=isdefault(it);
+ camp::file *f=defaultfile ? &camp::Stdout : vm::get<camp::file*>(it);
+ if(!f->isOpen()) return;
+ if(S != "") f->write(S);
+ if(f->text()) {
+ if(suffix) {
+ s->push(f);
+ suffix->call(s);
+ } else if(defaultfile) f->writeline();
+ }
+}
+
+void checkSquare(array *a)
+{
+ size_t n=checkArray(a);
+ for(size_t i=0; i < n; i++)
+ if(checkArray(read<array*>(a,i)) != n)
+ error("matrix a must be square");
+}
+
+// Crout's algorithm for computing the LU decomposition of a square matrix.
+// cf. routine ludcmp (Press et al., Numerical Recipes, 1991).
+Int LUdecompose(double *a, size_t n, size_t* index, bool warn=true)
+{
+ double *vv=new double[n];
+ Int swap=1;
+ for(size_t i=0; i < n; ++i) {
+ double big=0.0;
+ double *ai=a+i*n;
+ for(size_t j=0; j < n; ++j) {
+ double temp=fabs(ai[j]);
+ if(temp > big) big=temp;
+ }
+ if(big == 0.0) {
+ delete[] vv;
+ if(warn) error(singular);
+ else return 0;
+ }
+ vv[i]=1.0/big;
+ }
+ for(size_t j=0; j < n; ++j) {
+ for(size_t i=0; i < j; ++i) {
+ double *ai=a+i*n;
+ double sum=ai[j];
+ for(size_t k=0; k < i; ++k) {
+ sum -= ai[k]*a[k*n+j];
+ }
+ ai[j]=sum;
+ }
+ double big=0.0;
+ size_t imax=j;
+ for(size_t i=j; i < n; ++i) {
+ double *ai=a+i*n;
+ double sum=ai[j];
+ for(size_t k=0; k < j; ++k)
+ sum -= ai[k]*a[k*n+j];
+ ai[j]=sum;
+ double temp=vv[i]*fabs(sum);
+ if(temp >= big) {
+ big=temp;
+ imax=i;
+ }
+ }
+ double *aj=a+j*n;
+ double *aimax=a+imax*n;
+ if(j != imax) {
+ for(size_t k=0; k < n; ++k) {
+ double temp=aimax[k];
+ aimax[k]=aj[k];
+ aj[k]=temp;
+ }
+ swap *= -1;
+ vv[imax]=vv[j];
+ }
+ if(index)
+ index[j]=imax;
+ if(j != n) {
+ double denom=aj[j];
+ if(denom == 0.0) {
+ delete[] vv;
+ if(warn) error(singular);
+ else return 0;
+ }
+ for(size_t i=j+1; i < n; ++i)
+ a[i*n+j] /= denom;
+ }
+ }
+ delete[] vv;
+ return swap;
+}
+
+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);
+}
+
+#if defined(HAVE_LIBREADLINE) && defined(HAVE_LIBCURSES)
+struct historyState {
+ bool store;
+ HISTORY_STATE state;
+};
+
+typedef mem::map<CONST string, historyState> historyMap_t;
+historyMap_t historyMap;
+static HISTORY_STATE history_save;
+
+// Store a deep copy of the current readline history in dest.
+void store_history(HISTORY_STATE *dest)
+{
+ HISTORY_STATE *src=history_get_history_state();
+ if(src) {
+ *dest=*src;
+ for(Int i=0; i < src->length; ++i)
+ dest->entries[i]=src->entries[i];
+ free(src);
+ }
+}
+
+stringarray* get_history(Int n)
+{
+ int N=intcast(n);
+ if(N <= 0) N=history_length;
+ else N=Min(N,history_length);
+ array *a=new array((size_t) N);
+ int offset=history_length-N+1;
+ for(int i=0; i < N; ++i) {
+ HIST_ENTRY *last=history_get(offset+i);
+ string s=last ? last->line : "";
+ (*a)[i]=s;
+ }
+ return a;
+}
+
+string historyfilename(const string &name)
+{
+ return historyname+"_"+name;
+}
+#endif
+
+#if defined(HAVE_LIBREADLINE) && defined(HAVE_LIBCURSES)
+int readline_startup_hook()
+{
+#ifdef __CYGWIN__
+ rl_set_key("\\M-[3~",rl_delete,rl_get_keymap());
+ rl_set_key("\\M-[2~",rl_overwrite_mode,rl_get_keymap());
+#endif
+ return 0;
+}
+
+void init_readline(bool tabcompletion=true)
+{
+ static bool first=true;
+ if(first) {
+ first=false;
+#ifdef __CYGWIN__
+ rl_startup_hook=readline_startup_hook;
+#endif
+ }
+ rl_bind_key('\t',tabcompletion ? rl_complete : rl_insert);
+}
+#endif
+
+void cleanup()
+{
+ processDataStruct &pd=processData();
+ pd.atExitFunction=NULL;
+ pd.atUpdateFunction=NULL;
+ pd.atBreakpointFunction=NULL;
+
+#if defined(HAVE_LIBREADLINE) && defined(HAVE_LIBCURSES)
+ store_history(&history_save);
+ int nlines=intcast(getSetting<Int>("historylines"));
+ for(historyMap_t::iterator h=historyMap.begin(); h != historyMap.end();
+ ++h) {
+ history_set_history_state(&h->second.state);
+ stifle_history(nlines);
+ if(h->second.store) write_history(historyfilename(h->first).c_str());
+ }
+ history_set_history_state(&history_save);
+#endif
+}
+
+void purge(Int divisor=0)
+{
+#ifdef USEGC
+ if(divisor > 0) GC_set_free_space_divisor((GC_word) divisor);
+ GC_gcollect();
+#endif
+}
+
+void updateFunction(stack *Stack)
+{
+ callable *atUpdateFunction=processData().atUpdateFunction;
+ if(atUpdateFunction && !nullfunc::instance()->compare(atUpdateFunction))
+ atUpdateFunction->call(Stack);
+}
+
+void exitFunction(stack *Stack)
+{
+ callable *atExitFunction=processData().atExitFunction;
+ if(atExitFunction && !nullfunc::instance()->compare(atExitFunction))
+ atExitFunction->call(Stack);
+ cleanup();
+}
+
+default_t def;
+string emptystring;
+array *emptyarray=new array(0);
+string commentchar="#";
+pair zero;
+
+void breakpoint(stack *Stack, runnable *r)
+{
+ callable *atBreakpointFunction=processData().atBreakpointFunction;
+ if(atBreakpointFunction &&
+ !nullfunc::instance()->compare(atBreakpointFunction)) {
+ position curPos=getPos();
+ Stack->push<string>(curPos.filename());
+ Stack->push<Int>((Int) curPos.Line());
+ Stack->push<Int>((Int) curPos.Column());
+ Stack->push(r ? r : item(def));
+ atBreakpointFunction->call(Stack); // returns a string
+ } else Stack->push<string>("");
+}
+
+}
+
+namespace types {
+extern const char *names[];
+}
+
+void checkformat(const char *ptr, bool intformat)
+{
+ while(*ptr != '\0') {
+ if(*ptr != '%') /* While we have regular characters, print them. */
+ ptr++;
+ else { /* We've got a format specifier. */
+ ptr++;
+
+ while(*ptr && strchr ("-+ #0'I", *ptr)) /* Move past flags. */
+ *ptr++;
+
+ if(*ptr == '*')
+ ptr++;
+ else while(isdigit(*ptr)) /* Handle explicit numeric value. */
+ ptr++;
+
+ if(*ptr == '.') {
+ *ptr++; /* Go past the period. */
+ if(*ptr == '*') {
+ ptr++;
+ } else
+ while(isdigit(*ptr)) /* Handle explicit numeric value. */
+ *ptr++;
+ }
+ while(*ptr && strchr ("hlL", *ptr))
+ *ptr++;
+
+ if(*ptr == '%') ++ptr;
+ else if(*ptr != '\0') {
+ if(intformat) {
+ switch(*ptr) {
+ case 'd':
+ case 'i':
+ case 'o':
+ case 'u':
+ case 'x':
+ case 'X':
+ case 'c':
+ break;
+ default:
+ ostringstream buf;
+ buf << "Invalid format '" << *ptr << "' for type "
+ << types::names[types::ty_Int];
+ error(buf);
+ break;
+ }
+ } else {
+ switch(*ptr) {
+ case 'f':
+ case 'F':
+ case 'e':
+ case 'E':
+ case 'g':
+ case 'G':
+ break;
+ default:
+ ostringstream buf;
+ buf << "Invalid format '" << *ptr << "' for type "
+ << types::names[types::ty_real];
+ error(buf);
+ break;
+ }
+ }
+ }
+ } /* End of else statement */
+ }
+}
+
+// Return an angle in the interval [0,360).
+inline double principalBranch(double deg)
+{
+ if(deg < 0) deg += 360;
+ return deg;
+}
+
+static string defaulttransparency=string("Compatible");
+
+static const string defaulttimeformat=string("%a %b %d %T %Z %Y");
+#ifdef HAVE_STRFTIME
+static const size_t nTime=256;
+static char Time[nTime];
+#endif
+
+void clear(string file, Int line, bool warn=false)
+{
+ bpinfo bp(file,line);
+ for(mem::list<bpinfo>::iterator p=bplist.begin(); p != bplist.end(); ++p) {
+ if(*p == bp) {
+ cout << "cleared breakpoint at " << file << ": " << line << endl;
+ bplist.remove(bp);
+ return;
+ }
+ }
+ if(warn)
+ cout << "No such breakpoint at " << file << ": " << line << endl;
+}
+
+Int windingnumber(array *p, camp::pair z)
+{
+ size_t size=checkArray(p);
+ Int count=0;
+ for(size_t i=0; i < size; i++)
+ count += read<path *>(p,i)->windingnumber(z);
+ return count;
+}
+
+string convertname(string name, const string& format, bool check=true)
+{
+ if(name.empty())
+ return buildname(outname(),format,"");
+ else
+ if(check) checkLocal(name);
+ return format.empty() ? name : format+":"+name;
+}
+
+callable *Func;
+stack *FuncStack;
+double wrapFunction(double x)
+{
+ FuncStack->push(x);
+ Func->call(FuncStack);
+ return pop<double>(FuncStack);
+}
+
+callable *compareFunc;
+bool compareFunction(const vm::item& i, const vm::item& j)
+{
+ FuncStack->push(i);
+ FuncStack->push(j);
+ compareFunc->call(FuncStack);
+ return pop<bool>(FuncStack);
+}
+
+void unused(void *)
+{
+}
+
+pair readpair(stringstream& s, bool tex)
+{
+ static const double f=12.0/100.0;
+ double x,y;
+ s >> y;
+ s >> x;
+ return tex ? pair(f*x,-f*y) : pair(x,y);
+}
+
+// Ignore unclosed begingroups but not spurious endgroups.
+const char *nobegin="endgroup without matching begingroup";
+
+// Return the component of vector v perpendicular to a unit vector u.
+inline triple perp(triple v, triple u)
+{
+ return v-dot(v,u)*u;
+}
+
+string ASYo="/ASYo {( ) print 12 string cvs print} bind def";
+string pathforall="{(M) print ASYo ASYo} {(L) print ASYo ASYo} {(C) print ASYo ASYo ASYo ASYo ASYo ASYo} {(c) print} pathforall";
+string currentpoint="print currentpoint ASYo ASYo ";
+
+array *readpath(const string& psname, bool tex, bool keep)
+{
+ array *P=new array(0);
+ ostringstream pipe;
+ pipe << "'" << getSetting<string>("gs")
+ << "' -q -dNOPAUSE -dBATCH";
+ if(safe) pipe << " -dSAFER";
+#ifdef __CYGWIN__
+ string null="NUL";
+#else
+ string null="/dev/null";
+#endif
+ pipe << " -sDEVICE=epswrite -sOutputFile="+null+" '" << psname << "'";
+ iopipestream gs(pipe.str().c_str(),"gs","Ghostscript");
+ stringstream buf;
+ while(true) {
+ string out;
+ gs >> out;
+ buf << out;
+ string s=buf.str();
+ if(gs.tailequals(s.c_str(),s.size(),"E",1)) break;
+ }
+ if(verbose > 2) cout << endl;
+
+ mem::vector<solvedKnot> nodes;
+ solvedKnot node;
+ bool cyclic=false;
+ bool active=false;
+
+ pair offset;
+ while(!buf.eof()) {
+ char c;
+ buf >> c;
+ switch(c) {
+ case 'M':
+ {
+ if(active) {
+ if(cyclic) {
+ if(node.point == nodes[0].point)
+ nodes[0].pre=node.pre;
+ else {
+ pair delta=(nodes[0].point-node.point)*third;
+ node.post=node.point+delta;
+ nodes[0].pre=nodes[0].point-delta;
+ node.straight=true;
+ nodes.push_back(node);
+ }
+ } else {
+ node.post=node.point;
+ node.straight=false;
+ nodes.push_back(node);
+ }
+ P->push(path(nodes,nodes.size(),cyclic));
+ nodes.clear();
+ }
+ active=false;
+ cyclic=false;
+ node.pre=node.point=readpair(buf,tex)-offset;
+ node.straight=false;
+ break;
+ }
+ case 'L':
+ {
+ pair point=readpair(buf,tex)-offset;
+ pair delta=(point-node.point)*third;
+ node.post=node.point+delta;
+ node.straight=true;
+ nodes.push_back(node);
+ active=true;
+ node.pre=point-delta;
+ node.point=point;
+ break;
+ }
+ case 'C':
+ {
+ pair point=readpair(buf,tex)-offset;
+ pair pre=readpair(buf,tex)-offset;
+ node.post=readpair(buf,tex)-offset;
+ node.straight=false;
+ nodes.push_back(node);
+ active=true;
+ node.pre=pre;
+ node.point=point;
+ break;
+ }
+ case 'c':
+ {
+ cyclic=true;
+ break;
+ }
+ case 'Z':
+ {
+ offset=readpair(buf,tex);
+ break;
+ }
+ }
+ }
+
+ if(!keep)
+ unlink(psname.c_str());
+ return P;
+}
+
+pair sin(pair z)
+{
+ return pair(sin(z.getx())*cosh(z.gety()),cos(z.getx())*sinh(z.gety()));
+}
+
+pair exp(pair z)
+{
+ return exp(z.getx())*expi(z.gety());
+}
+
+pair gamma(pair z)
+{
+ static double p[]={0.99999999999980993,676.5203681218851,-1259.1392167224028,
+ 771.32342877765313,-176.61502916214059,12.507343278686905,
+ -0.13857109526572012,9.9843695780195716e-6,
+ 1.5056327351493116e-7};
+ static int n=sizeof(p)/sizeof(double);
+ static double root2pi=sqrt(2*PI);
+ if(z.getx() < 0.5)
+ return PI/(sin(PI*z)*gamma(1.0-z));
+ z -= 1.0;
+ pair x=p[0];
+ for(int i=1; i < n; ++i)
+ x += p[i]/(z+i);
+ pair t=n-1.5+z;
+ return root2pi*pow(t,z+0.5)*exp(-t)*x;
+}
+
+// Autogenerated routines:
+
+
+
+// Initializers
+
+Int :IntZero()
+{
+ return 0;
+}
+
+real :realZero()
+{
+ return 0.0;
+}
+
+bool :boolFalse()
+{
+ return false;
+}
+
+array* :pushNullArray()
+{
+ return 0;
+}
+
+frame* :pushNullRecord()
+{
+ return 0;
+}
+
+item :pushNullFunction()
+{
+ return nullfunc::instance();
+}
+
+
+// Default operations
+
+// Put the default value token on the stack (in place of an argument when
+// making a function call).
+item :pushDefault()
+{
+ return def;
+}
+
+
+// Test if the value on the stack is the default value token.
+bool :isDefault(item i)
+{
+ return isdefault(i);
+}
+
+
+// Logical operations
+
+bool !(bool b)
+{
+ return !b;
+}
+
+bool :boolMemEq(frame *a, frame *b)
+{
+ return a == b;
+}
+
+bool :boolMemNeq(frame *a, frame *b)
+{
+ return a != b;
+}
+
+bool :boolFuncEq(callable *a, callable *b)
+{
+ return a->compare(b);
+}
+
+bool :boolFuncNeq(callable *a, callable *b)
+{
+ return !(a->compare(b));
+}
+
+
+// Bit operations
+
+Int AND(Int a, Int b)
+{
+ return a & b;
+}
+
+Int OR(Int a, Int b)
+{
+ return a | b;
+}
+
+Int XOR(Int a, Int b)
+{
+ return a ^ b;
+}
+
+Int NOT(Int a)
+{
+ return ~a;
+}
+
+
+// Casts
+
+guide* :pairToGuide(pair z)
+{
+ return new pairguide(z);
+}
+
+guide* :pathToGuide(path p)
+{
+ return new pathguide(p);
+}
+
+path :guideToPath(guide *g)
+{
+ return g->solve();
+}
+
+
+// Picture operations
+
+picture* :newPicture()
+{
+ return new picture();
+}
+
+bool empty(picture *f)
+{
+ return f->null();
+}
+
+void erase(picture *f)
+{
+ f->nodes.clear();
+}
+
+pair min(picture *f)
+{
+ return f->bounds().Min();
+}
+
+pair max(picture *f)
+{
+ return f->bounds().Max();
+}
+
+void label(picture *f, string *s, string *size, transform t, pair position,
+ pair align, pen p)
+{
+ f->append(new drawLabel(*s,*size,t,position,align,p));
+}
+
+bool labels(picture *f)
+{
+ return f->havelabels();
+}
+
+patharray *_texpath(string *s, pen p=CURRENTPEN)
+{
+ array *P=new array(0);
+ if(s->empty()) return P;
+
+ string prefix=outname();
+ spaceToUnderscore(prefix);
+ string psname=auxname(prefix,"ps");
+
+ string texname=auxname(prefix,"tex");
+ string dviname=auxname(prefix,"dvi");
+ bbox b;
+ texfile tex(texname,b,true);
+ tex.miniprologue();
+
+ tex.setfont(p);
+ tex.verbatimline("\\special{ps:");
+ tex.verbatimline(ASYo);
+ tex.verbatimline("/ASY1 true def");
+ tex.verbatimline("/v {neg exch 4 copy 4 2 roll 2 copy 6 2 roll 2 copy (M) print ASYo ASYo (L) print ASYo add ASYo (L) print add ASYo add ASYo (L) print add ASYo ASYo (c) print} bind def");
+ tex.verbatimline("/show {ASY1 {(Z) "+currentpoint+
+ "/ASY1 false def} if currentpoint newpath moveto false charpath "+
+ pathforall+"} bind def}");
+ tex.verbatim(*s);
+ tex.verbatimline("\\special{ps: (E) print}");
+ tex.epilogue(true);
+ tex.close();
+
+ ostringstream cmd;
+ string aux=auxname(prefix,"aux");
+ unlink(aux.c_str());
+ cmd << texprogram(true) << " \\nonstopmode\\input '" << texname << "'";
+ bool quiet=verbose <= 2;
+ int status=System(cmd,quiet ? 1 : 0,"texpath",texpathmessage());
+ if(!status && getSetting<bool>("twice"))
+ status=System(cmd,quiet ? 1 : 0,"texpath",texpathmessage());
+ if(status) {
+ if(quiet) {
+ ostringstream cmd;
+ cmd << texprogram(true) << " \\scrollmode\\input '" << texname << "'";
+ System(cmd,0);
+ }
+ }
+
+ if(!status) {
+ ostringstream cmd;
+ cmd << "'" << getSetting<string>("dvips") << "' -R -Pdownload35 -D600 "
+ << getSetting<string>("dvipsOptions");
+ if(verbose <= 2) cmd << " -q";
+ cmd << " -o '" << psname << "' '" << dviname << "'";
+ status=System(cmd,0,true,"dvips");
+ }
+
+ if(status != 0)
+ error("texpath failed");
+
+ bool keep=getSetting<bool>("keep");
+ if(!keep) { // Delete temporary files.
+ unlink(texname.c_str());
+ if(!getSetting<bool>("keepaux"))
+ unlink(aux.c_str());
+ unlink(auxname(prefix,"log").c_str());
+ unlink(dviname.c_str());
+ }
+ return readpath(psname,true,keep);
+}
+
+patharray *_strokepath(path g, pen p=CURRENTPEN)
+{
+ array *P=new array(0);
+ if(g.size() == 0) return P;
+
+ string prefix=outname();
+ spaceToUnderscore(prefix);
+ string psname=auxname(prefix,"ps");
+ bbox b;
+ psfile ps(psname, false);
+ ps.prologue(b);
+ ps.verbatimline(ASYo);
+ ps.verbatimline("/stroke {(Z) "+currentpoint+pathforall+"} bind def");
+ ps.resetpen();
+ ps.setpen(p);
+ ps.write(g);
+ ps.strokepath();
+ ps.stroke();
+ ps.verbatimline("(M) "+currentpoint+"(E) print");
+ ps.epilogue();
+ ps.close();
+ return readpath(psname,false,getSetting<bool>("keep"));
+}
+
+void _draw(picture *f, path g, pen p)
+{
+ f->append(new drawPath(g,p));
+}
+
+void fill(picture *f, patharray *g, pen p=CURRENTPEN, bool copy=true)
+{
+ array *(*copyarray)(array *a)=copy ? copyArray: nop;
+ f->append(new drawFill(*copyarray(g),false,p));
+}
+
+void latticeshade(picture *f, patharray *g, bool stroke=false,
+ pen fillrule=CURRENTPEN, penarray2 *p, bool copy=true)
+{
+ array *(*copyarray)(array *a)=copy ? copyArray: nop;
+ f->append(new drawLatticeShade(*copyarray(g),stroke,fillrule,*copyarray(p)));
+}
+
+void axialshade(picture *f, patharray *g, bool stroke=false, pen pena, pair a,
+ pen penb, pair b, bool copy=true)
+{
+ array *(*copyarray)(array *a)=copy ? copyArray: nop;
+ f->append(new drawAxialShade(*copyarray(g),stroke,pena,a,penb,b));
+}
+
+void radialshade(picture *f, patharray *g, bool stroke=false, pen pena,
+ pair a, real ra, pen penb, pair b, real rb, bool copy=true)
+{
+ array *(*copyarray)(array *a)=copy ? copyArray: nop;
+ f->append(new drawRadialShade(*copyarray(g),stroke,pena,a,ra,penb,b,rb));
+}
+
+void gouraudshade(picture *f, patharray *g, bool stroke=false,
+ pen fillrule=CURRENTPEN, penarray *p, pairarray *z,
+ Intarray *edges, bool copy=true)
+{
+ array *(*copyarray)(array *a)=copy ? copyArray: nop;
+ checkArrays(p,z);
+ checkArrays(z,edges);
+ f->append(new drawGouraudShade(*copyarray(g),stroke,fillrule,*copyarray(p),
+ *copyarray(z),*copyarray(edges)));
+}
+
+void gouraudshade(picture *f, patharray *g, bool stroke=false,
+ pen fillrule=CURRENTPEN, penarray *p, Intarray *edges,
+ bool copy=true)
+{
+ array *(*copyarray)(array *a)=copy ? copyArray: nop;
+ size_t n=checkArrays(p,edges);
+ size_t m=checkArray(g);
+ array *z=new array(n);
+ Int k=0;
+ Int in=(Int) n;
+ for(size_t j=0; j < m; ++j) {
+ path *P=read<path *>(g,j);
+ assert(P);
+ Int stop=Min(P->size(),in-k);
+ mem::vector<solvedKnot>& nodes=P->Nodes();
+ for(Int i=0; i < stop; ++i)
+ (*z)[k++]=nodes[i].point;
+ }
+ checkArrays(p,z);
+
+ f->append(new drawGouraudShade(*copyarray(g),stroke,fillrule,*copyarray(p),
+ *z,*copyarray(edges)));
+}
+
+void tensorshade(picture *f, patharray *g, bool stroke=false,
+ pen fillrule=CURRENTPEN, penarray2 *p, patharray *b=NULL,
+ pairarray2 *z=emptyarray, bool copy=true)
+{
+ array *(*copyarray)(array *a)=copy ? copyArray: nop;
+ array *(*copyarray2)(array *a)=copy ? copyArray2: nop;
+ if(b == NULL) b=g;
+ size_t n=checkArrays(p,b);
+ size_t nz=checkArray(z);
+ if(nz != 0)
+ checkEqual(nz,n);
+ f->append(new drawTensorShade(*copyarray(g),stroke,fillrule,*copyarray2(p),
+ *copyarray(b),*copyarray2(z)));
+}
+
+void functionshade(picture *f, patharray *g, bool stroke=false,
+ pen fillrule=CURRENTPEN, string shader=emptystring,
+ bool copy=true)
+{
+ array *(*copyarray)(array *a)=copy ? copyArray: nop;
+ f->append(new drawFunctionShade(*copyarray(g),stroke,fillrule,shader));
+}
+
+// Clip a picture to a superpath using the given fill rule.
+// Subsequent additions to the picture will not be affected by the clipping.
+void clip(picture *f, patharray *g, bool stroke=false,
+ pen fillrule=CURRENTPEN, bool copy=true)
+{
+ array *(*copyarray)(array *a)=copy ? copyArray: nop;
+ drawClipBegin *begin=new drawClipBegin(*copyarray(g),stroke,fillrule,true);
+ f->enclose(begin,new drawClipEnd(true,begin));
+}
+
+void beginclip(picture *f, patharray *g, bool stroke=false,
+ pen fillrule=CURRENTPEN, bool copy=true)
+{
+ array *(*copyarray)(array *a)=copy ? copyArray: nop;
+ f->append(new drawClipBegin(*copyarray(g),stroke,fillrule,false));
+}
+
+void endclip(picture *f)
+{
+ f->append(new drawClipEnd(false));
+}
+
+void gsave(picture *f)
+{
+ f->append(new drawGsave());
+}
+
+void grestore(picture *f)
+{
+ f->append(new drawGrestore());
+}
+
+void begingroup(picture *f)
+{
+ f->append(new drawBegin());
+}
+
+void endgroup(picture *f)
+{
+ f->append(new drawEnd());
+}
+
+void add(picture *dest, picture *src)
+{
+ dest->add(*src);
+}
+
+void prepend(picture *dest, picture *src)
+{
+ dest->prepend(*src);
+}
+
+void postscript(picture *f, string s)
+{
+ f->append(new drawVerbatim(PostScript,s));
+}
+
+void tex(picture *f, string s)
+{
+ f->append(new drawVerbatim(TeX,s));
+}
+
+void postscript(picture *f, string s, pair min, pair max)
+{
+ f->append(new drawVerbatim(PostScript,s,min,max));
+}
+
+void tex(picture *f, string s, pair min, pair max)
+{
+ f->append(new drawVerbatim(TeX,s,min,max));
+}
+
+void texpreamble(string s)
+{
+ string t=s+"\n";
+ processDataStruct &pd=processData();
+ pd.TeXpipepreamble.push_back(t);
+ pd.TeXpreamble.push_back(t);
+}
+
+void deletepreamble()
+{
+ if(getSetting<bool>("inlinetex")) {
+ unlink(auxname(outname(),"pre").c_str());
+ }
+}
+
+void _labelpath(picture *f, string s, string size, path g, string justify,
+ pair offset, pen p)
+{
+ f->append(new drawLabelPath(s,size,g,justify,offset,p));
+}
+
+void texreset()
+{
+ processDataStruct &pd=processData();
+ pd.TeXpipepreamble.clear();
+ pd.TeXpreamble.clear();
+ pd.tex.pipeclose();
+}
+
+void layer(picture *f)
+{
+ f->append(new drawLayer());
+}
+
+void newpage(picture *f)
+{
+ f->append(new drawNewPage());
+}
+
+void _image(picture *f, realarray2 *data, pair initial, pair final,
+ penarray *palette=NULL, transform t=identity, bool copy=true,
+ bool antialias=false)
+{
+ array *(*copyarray)(array *a)=copy ? copyArray: nop;
+ array *(*copyarray2)(array *a)=copy ? copyArray2: nop;
+ f->append(new drawImage(*copyarray2(data),*copyarray(palette),
+ t*matrix(initial,final),antialias));
+}
+
+void _image(picture *f, penarray2 *data, pair initial, pair final,
+ transform t=identity, bool copy=true, bool antialias=false)
+{
+ array *(*copyarray2)(array *a)=copy ? copyArray2: nop;
+ f->append(new drawImage(*copyarray2(data),t*matrix(initial,final),antialias));
+}
+
+string nativeformat()
+{
+ return nativeformat();
+}
+
+bool latex()
+{
+ return latex(getSetting<string>("tex"));
+}
+
+bool pdf()
+{
+ return pdf(getSetting<string>("tex"));
+}
+
+void shipout(string prefix=emptystring, picture *f, picture *preamble=NULL,
+ string format=emptystring, bool wait=false, bool view=true,
+ callableTransform *xform)
+{
+ if(prefix.empty()) prefix=outname();
+
+ picture *result=new picture;
+ unsigned level=0;
+ picture::nodelist::iterator p;
+ for(p = f->nodes.begin(); p != f->nodes.end(); ++p) {
+ xform->call(Stack);
+ transform t=pop<transform>(Stack);
+ static transform Zero=transform(0.0,0.0,0.0,0.0,0.0,0.0);
+ bool Delete=(t == Zero);
+ picture *group=new picture;
+ assert(*p);
+ if((*p)->endgroup()) error(nobegin);
+ if((*p)->begingroup()) {
+ ++level;
+ while(p != f->nodes.end() && level) {
+ if(!Delete) {
+ drawElement *e=t.isIdentity() ? *p : (*p)->transformed(t);
+ group->append(e);
+ }
+ ++p;
+ if(p == f->nodes.end()) break;
+ assert(*p);
+ if((*p)->begingroup()) ++level;
+ if((*p)->endgroup()) {
+ if(level) --level;
+ else error(nobegin);
+ }
+ }
+ }
+ if(p == f->nodes.end()) break;
+ assert(*p);
+ if(!Delete) {
+ drawElement *e=t.isIdentity() ? *p : (*p)->transformed(t);
+ group->append(e);
+ result->add(*group);
+ }
+ }
+
+ result->shipout(preamble,prefix,format,0.0,wait,view);
+}
+
+void shipout3(string prefix, picture *f, string format=emptystring,
+ real width, real height, real angle, triple m, triple M,
+ realarray2 *t, triplearray *lights, realarray2 *diffuse,
+ realarray2 *ambient, realarray2 *specular, bool viewportlighting,
+ bool view=true)
+{
+ size_t n=checkArrays(lights,diffuse);
+ checkEqual(n,checkArray(ambient));
+ checkEqual(n,checkArray(specular));
+
+ f->shipout3(prefix,format,width,height,angle,m,M,copyArray2C(t,true,4),n,
+ copyTripleArrayC(lights),copyArray2C(diffuse,false,4),
+ copyArray2C(ambient,false,4),copyArray2C(specular,false,4),
+ viewportlighting,view);
+}
+
+void shipout3(string prefix, picture *f)
+{
+ f->shipout3(prefix);
+}
+
+void deconstruct(picture *f, picture *preamble=NULL, real magnification=1,
+ callableTransform *xform)
+{
+ unsigned level=0;
+ unsigned n=0;
+
+ string prefix=outname();
+ string xformat=getSetting<string>("xformat");
+
+ static long arg_max=sysconf(_SC_ARG_MAX);
+ const unsigned maxargs=::min(arg_max/(prefix.size()+xformat.size()+25ul),
+ 256ul);
+
+ cout << xformat << newl;
+ cout << maxargs << newl;
+
+ string preformat=nativeformat();
+ const string Done="Done";
+ const string Error="Error";
+
+ ostringstream cmd;
+
+ // Enforce ghostscript limitations.
+ magnification=::max(magnification,0.0001);
+ real res=::min(::max(magnification*72.0,2.0),8192.0);
+
+ const char *converter=NULL, *hint=NULL;
+
+ bool png=xformat == "png";
+
+ if(magnification > 0.0) {
+ mem::list<string> nameStack;
+ string outname;
+ unsigned arg=0;
+ unsigned batch=0;
+ for(picture::nodelist::iterator p=f->nodes.begin();;) {
+ if(p == f->nodes.end()) break;
+ if(arg == 0) {
+ cmd.str("");
+ ostringstream buf;
+ buf << batch << "_";
+ outname=buildname(prefix+buf.str()+"%d",xformat,"");
+ if(png) {
+ cmd << "'" << getSetting<string>("gs")
+ << "' -q -dNOPAUSE -dBATCH -sDEVICE=pngalpha -dEPSCrop";
+ if(safe)
+ cmd << " -dSAFER";
+ cmd << " -r" << res << "x" << res
+ << " -sOutputFile='" << outname << "'";
+ converter="gs";
+ hint="Ghostscript";
+ } else {
+ cmd << "'" << getSetting<string>("convert")
+ << "' -density " << res << "x" << res
+ << " -transparent white";
+ hint=converter="convert";
+ }
+ }
+
+ picture *group=new picture;
+ xform->call(Stack);
+ transform t=pop<transform>(Stack);
+ assert(*p);
+ if((*p)->endgroup()) {
+ cout << Error << endl;
+ error(nobegin);
+ }
+ if((*p)->begingroup()) {
+ ++level;
+ while(p != f->nodes.end() && level) {
+ drawElement *e=t.isIdentity() ? *p : (*p)->transformed(t);
+ group->append(e);
+ ++p;
+ if(p == f->nodes.end()) break;
+ assert(*p);
+ if((*p)->begingroup()) ++level;
+ if((*p)->endgroup()) {
+ if(level) --level;
+ else {
+ cout << Error << endl;
+ error(nobegin);
+ }
+ }
+ }
+ }
+ if(p != f->nodes.end()) {
+ assert(*p);
+ drawElement *e=t.isIdentity() ? *p : (*p)->transformed(t);
+ group->append(e);
+ bbox b;
+ ostringstream buf;
+ buf << prefix << "_" << n;
+ group->shipout(preamble,buf.str(),preformat,magnification,false,false);
+ string Preformat=png && group->Transparency() ? "pdf" : preformat;
+ string name=buildname(buf.str(),Preformat);
+ nameStack.push_back(name);
+ cmd << " '";
+ if(!png) cmd << preformat+":";
+ cmd << name << "'";
+ b=group->bounds();
+ b *= magnification;
+
+ cout << b << newl;
+ ++n;
+ ++p;
+ ++arg;
+ }
+
+ if(p == f->nodes.end() || arg >= maxargs) {
+ arg=0;
+ ++batch;
+ cout.flush();
+ if(!png)
+ cmd << " -scene 1 '" << xformat << ":" << outname << "'";
+
+ int status=System(cmd,0,true,converter,hint);
+ if(status) {
+ cout << Error << endl;
+ error("deconstruct failed");
+ }
+ }
+ }
+
+ if(!getSetting<bool>("keep")) {
+ for(mem::list<string>::iterator p=nameStack.begin();
+ p != nameStack.end(); ++p)
+ unlink(p->c_str());
+ }
+
+ cout << Done << endl;
+ }
+}
+
+void purge(Int divisor=0)
+{
+ purge(divisor);
+}
+
+
+// Pen operations
+
+pen :newPen()
+{
+ return pen();
+}
+
+bool ==(pen a, pen b)
+{
+ return a == b;
+}
+
+bool !=(pen a, pen b)
+{
+ return a != b;
+}
+
+pen +(pen a, pen b)
+{
+ return a+b;
+}
+
+pen Operator *(real a, pen b)
+{
+ return a*b;
+}
+
+pen Operator *(pen a, real b)
+{
+ return b*a;
+}
+
+pair max(pen p)
+{
+ return p.bounds().Max();
+}
+
+pair min(pen p)
+{
+ return p.bounds().Min();
+}
+
+// Reset the meaning of pen default attributes.
+void resetdefaultpen()
+{
+ processData().defaultpen=camp::pen::initialpen();
+}
+
+void defaultpen(pen p)
+{
+ processData().defaultpen=pen(resolvepen,p);
+}
+
+pen defaultpen()
+{
+ return processData().defaultpen;
+}
+
+bool invisible(pen p)
+{
+ return p.invisible();
+}
+
+pen invisible()
+{
+ return pen(invisiblepen);
+}
+
+pen gray(pen p)
+{
+ p.togrey();
+ return p;
+}
+
+pen rgb(pen p)
+{
+ p.torgb();
+ return p;
+}
+
+pen cmyk(pen p)
+{
+ p.tocmyk();
+ return p;
+}
+
+pen interp(pen a, pen b, real t)
+{
+ return interpolate(a,b,t);
+}
+
+pen rgb(real r, real g, real b)
+{
+ return pen(r,g,b);
+}
+
+pen cmyk(real c, real m, real y, real k)
+{
+ return pen(c,m,y,k);
+}
+
+pen gray(real gray)
+{
+ return pen(gray);
+}
+
+realarray *colors(pen p)
+{
+ size_t n=ColorComponents[p.colorspace()];
+ array *a=new array(n);
+
+ switch(n) {
+ case 0:
+ break;
+ case 1:
+ (*a)[0]=p.gray();
+ break;
+ case 3:
+ (*a)[0]=p.red();
+ (*a)[1]=p.green();
+ (*a)[2]=p.blue();
+ break;
+ case 4:
+ (*a)[0]=p.cyan();
+ (*a)[1]=p.magenta();
+ (*a)[2]=p.yellow();
+ (*a)[3]=p.black();
+ break;
+ default:
+ break;
+ }
+ return a;
+}
+
+string colorspace(pen p)
+{
+ string s=ColorDeviceSuffix[p.colorspace()];
+ std::transform(s.begin(),s.end(),s.begin(),tolower);
+ return s;
+}
+
+pen pattern(string *s)
+{
+ return pen(setpattern,*s);
+}
+
+string pattern(pen p)
+{
+ return p.fillpattern();
+}
+
+pen fillrule(Int n)
+{
+ return pen(n >= 0 && n < nFill ? (FillRule) n : DEFFILL);
+}
+
+Int fillrule(pen p)
+{
+ return p.Fillrule();
+}
+
+pen opacity(real opacity=1.0, string blend=defaulttransparency)
+{
+ for(Int i=0; i < nBlendMode; ++i)
+ if(blend == BlendMode[i]) return pen(Transparency(blend,opacity));
+
+ ostringstream buf;
+ buf << "Unknown blend mode: " << "'" << blend << "'";
+ error(buf);
+}
+
+real opacity(pen p)
+{
+ return p.opacity();
+}
+
+string blend(pen p)
+{
+ return p.blend();
+}
+
+pen linetype(string *s, real offset=0, bool scale=true, bool adjust=true)
+{
+ return pen(LineType(*s,offset,scale,adjust));
+}
+
+string linetype(pen p=CURRENTPEN)
+{
+ return p.stroke();
+}
+
+real offset(pen p)
+{
+ return p.linetype().offset;
+}
+
+bool scale(pen p)
+{
+ return p.linetype().scale;
+}
+
+bool adjust(pen p)
+{
+ return p.linetype().adjust;
+}
+
+pen adjust(pen p, real arclength, bool cyclic)
+{
+ return adjustdash(p,arclength,cyclic);
+}
+
+pen linecap(Int n)
+{
+ return pen(setlinecap,n >= 0 && n < nCap ? n : DEFCAP);
+}
+
+Int linecap(pen p=CURRENTPEN)
+{
+ return p.cap();
+}
+
+pen linejoin(Int n)
+{
+ return pen(setlinejoin,n >= 0 && n < nJoin ? n : DEFJOIN);
+}
+
+Int linejoin(pen p=CURRENTPEN)
+{
+ return p.join();
+}
+
+pen miterlimit(real x)
+{
+ return pen(setmiterlimit,x >= 1.0 ? x : DEFJOIN);
+}
+
+real miterlimit(pen p=CURRENTPEN)
+{
+ return p.miter();
+}
+
+pen linewidth(real x)
+{
+ return pen(setlinewidth,x >= 0.0 ? x : DEFWIDTH);
+}
+
+real linewidth(pen p=CURRENTPEN)
+{
+ return p.width();
+}
+
+pen fontcommand(string *s)
+{
+ return pen(setfont,*s);
+}
+
+string font(pen p=CURRENTPEN)
+{
+ return p.Font();
+}
+
+pen fontsize(real size, real lineskip)
+{
+ return pen(setfontsize,size > 0.0 ? size : 0.0,
+ lineskip > 0.0 ? lineskip : 0.0);
+}
+
+real fontsize(pen p=CURRENTPEN)
+{
+ return p.size();
+}
+
+real lineskip(pen p=CURRENTPEN)
+{
+ return p.Lineskip();
+}
+
+pen overwrite(Int n)
+{
+ return pen(setoverwrite,n >= 0 && n < nOverwrite ? (overwrite_t) n :
+ DEFWRITE);
+}
+
+Int overwrite(pen p=CURRENTPEN)
+{
+ return p.Overwrite();
+}
+
+pen basealign(Int n)
+{
+ return pen(n >= 0 && n < nBaseLine ? (BaseLine) n : DEFBASE);
+}
+
+Int basealign(pen p=CURRENTPEN)
+{
+ return p.Baseline();
+}
+
+transform transform(pen p)
+{
+ return p.getTransform();
+}
+
+path nib(pen p)
+{
+ return p.Path();
+}
+
+pen makepen(path p)
+{
+ return pen(p);
+}
+
+pen colorless(pen p)
+{
+ p.colorless();
+ return p;
+}
+
+// Interactive mode
+
+bool interactive()
+{
+ return interact::interactive;
+}
+
+bool uptodate()
+{
+ return interact::uptodate;
+}
+
+
+// System commands
+
+Int system(string s)
+{
+ if(safe) error("system() call disabled; override with option -nosafe");
+ if(s.empty()) return 0;
+ else return System(s.c_str());
+}
+
+bool view()
+{
+ return view();
+}
+
+string asydir()
+{
+ return systemDir;
+}
+
+string locale(string s=emptystring)
+{
+ char *L=setlocale(LC_ALL,s.empty() ? NULL : s.c_str());
+ return L != NULL ? string(L) : "";
+}
+
+void abort(string s=emptystring)
+{
+ if(s.empty()) throw handled_error();
+ error(s.c_str());
+}
+
+void exit()
+{
+ throw quit();
+}
+
+void assert(bool b, string s=emptystring)
+{
+ flush(cout);
+ if(!b) {
+ ostringstream buf;
+ buf << "assert FAILED";
+ if(s != "") buf << ": " << s << endl;
+ error(buf);
+ }
+}
+
+void sleep(Int seconds)
+{
+ if(seconds <= 0) return;
+ sleep(seconds);
+}
+
+void usleep(Int microseconds)
+{
+ if(microseconds <= 0) return;
+ usleep((unsigned long) microseconds);
+}
+
+void _eval(string *s, bool embedded, bool interactiveWrite=false)
+{
+ if (embedded) {
+ trans::coenv *e=Stack->getEnvironment();
+ vm::interactiveStack *is=dynamic_cast<vm::interactiveStack *>(Stack);
+ if (e && is) {
+ runStringEmbedded(*s, *e, *is);
+ } else {
+ cerr << "no runtime environment for embedded eval" << endl;
+ }
+ } else {
+ runString(*s,interactiveWrite);
+ }
+}
+
+void _eval(runnable *s, bool embedded)
+{
+ absyntax::block *ast=new absyntax::block(s->getPos(), false);
+ ast->add(s);
+
+ if (embedded) {
+ trans::coenv *e=Stack->getEnvironment();
+ vm::interactiveStack *is=dynamic_cast<vm::interactiveStack *>(Stack);
+ if (e && is) {
+ runCodeEmbedded(ast, *e, *is);
+ } else {
+ cerr << "no runtime environment for embedded eval" << endl;
+ }
+ } else {
+ runCode(ast);
+ }
+}
+
+string location() {
+ ostringstream buf;
+ buf << getPos();
+ return buf.str();
+}
+
+// Wrapper for the stack::load() method.
+void :loadModule(string *index)
+{
+ Stack->load(*index);
+}
+
+string cd(string s=emptystring)
+{
+ if(!s.empty() && !globalwrite()) writeDisabled();
+ return setPath(s.c_str());
+}
+
+void list(string *s, bool imports=false)
+{
+ if(*s == "-") return;
+ trans::genv ge;
+ symbol *name=symbol::trans(*s);
+ record *r=ge.getModule(name,*s);
+ r->e.list(imports ? 0 : r);
+}
+
+
+// Path operations
+
+path :nullPath()
+{
+ return nullpath;
+}
+
+bool ==(path a, path b)
+{
+ return a == b;
+}
+
+bool !=(path a, path b)
+{
+ return !(a == b);
+}
+
+pair point(path p, Int t)
+{
+ return p.point((Int) t);
+}
+
+pair point(path p, real t)
+{
+ return p.point(t);
+}
+
+pair precontrol(path p, Int t)
+{
+ return p.precontrol((Int) t);
+}
+
+pair precontrol(path p, real t)
+{
+ return p.precontrol(t);
+}
+
+pair postcontrol(path p, Int t)
+{
+ return p.postcontrol((Int) t);
+}
+
+pair postcontrol(path p, real t)
+{
+ return p.postcontrol(t);
+}
+
+pair dir(path p, Int t, Int sign=0, bool normalize=true)
+{
+ return p.dir(t,sign,normalize);
+}
+
+pair dir(path p, real t, bool normalize=true)
+{
+ return p.dir(t,normalize);
+}
+
+pair accel(path p, Int t, Int sign=0)
+{
+ return p.accel(t,sign);
+}
+
+pair accel(path p, real t)
+{
+ return p.accel(t);
+}
+
+real radius(path p, real t)
+{
+ pair v=p.dir(t,false);
+ pair a=p.accel(t);
+ real d=dot(a,v);
+ real v2=v.abs2();
+ real a2=a.abs2();
+ real denom=v2*a2-d*d;
+ real r=v2*sqrt(v2);
+ return denom > 0 ? r/sqrt(denom) : 0.0;
+}
+
+path reverse(path p)
+{
+ return p.reverse();
+}
+
+path subpath(path p, Int a, Int b)
+{
+ return p.subpath((Int) a, (Int) b);
+}
+
+path subpath(path p, real a, real b)
+{
+ return p.subpath(a,b);
+}
+
+path nurb(pair z0, pair z1, pair z2, pair z3,
+ real w0, real w1, real w2, real w3, Int m)
+{
+ return nurb(z0,z1,z2,z3,w0,w1,w2,w3,m);
+}
+
+Int length(path p)
+{
+ return p.length();
+}
+
+bool cyclic(path p)
+{
+ return p.cyclic();
+}
+
+bool straight(path p, Int t)
+{
+ return p.straight(t);
+}
+
+path unstraighten(path p)
+{
+ return p.unstraighten();
+}
+
+bool piecewisestraight(path p)
+{
+ return p.piecewisestraight();
+}
+
+real arclength(path p)
+{
+ return p.arclength();
+}
+
+real arctime(path p, real dval)
+{
+ return p.arctime(dval);
+}
+
+real dirtime(path p, pair z)
+{
+ return p.directiontime(z);
+}
+
+realarray* intersect(path p, path q, real fuzz=-1)
+{
+ bool exact=fuzz <= 0.0;
+ if(fuzz < 0)
+ fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())),
+ ::max(length(q.max()),length(q.min())));
+ std::vector<real> S,T;
+ real s,t;
+ if(intersections(s,t,S,T,p,q,fuzz,true,exact)) {
+ array *V=new array(2);
+ (*V)[0]=s;
+ (*V)[1]=t;
+ return V;
+ } else
+ return new array(0);
+}
+
+realarray2* intersections(path p, path q, real fuzz=-1)
+{
+ bool exact=fuzz <= 0.0;
+ if(fuzz < 0.0)
+ fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())),
+ ::max(length(q.max()),length(q.min())));
+ real s,t;
+ std::vector<real> S,T;
+ intersections(s,t,S,T,p,q,fuzz,false,true);
+ size_t n=S.size();
+ if(n == 0 && !exact) {
+ if(intersections(s,t,S,T,p,q,fuzz,true,false)) {
+ array *V=new array(1);
+ array *Vi=new array(2);
+ (*V)[0]=Vi;
+ (*Vi)[0]=s;
+ (*Vi)[1]=t;
+ return V;
+ }
+ }
+ array *V=new array(n);
+ for(size_t i=0; i < n; ++i) {
+ array *Vi=new array(2);
+ (*V)[i]=Vi;
+ (*Vi)[0]=S[i];
+ (*Vi)[1]=T[i];
+ }
+ stable_sort(V->begin(),V->end(),run::compare2<real>());
+ return V;
+}
+
+realarray* intersections(path p, explicit pair a, explicit pair b, real fuzz=-1)
+{
+ if(fuzz < 0)
+ fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())),
+ ::max(length(a),length(b)));
+ std::vector<real> S;
+ intersections(S,p,a,b,fuzz);
+ sort(S.begin(),S.end());
+ size_t n=S.size();
+ array *V=new array(n);
+ for(size_t i=0; i < n; ++i)
+ (*V)[i]=S[i];
+ return V;
+}
+
+// Return the intersection point of the extensions of the line segments
+// PQ and pq.
+pair extension(pair P, pair Q, pair p, pair q)
+{
+ pair ac=P-Q;
+ pair bd=q-p;
+ real det=ac.getx()*bd.gety()-ac.gety()*bd.getx();
+ if(det == 0) return pair(infinity,infinity);
+ return P+((p.getx()-P.getx())*bd.gety()-(p.gety()-P.gety())*bd.getx())*ac/det;
+}
+
+Int size(path p)
+{
+ return p.size();
+}
+
+path &(path p, path q)
+{
+ return camp::concat(p,q);
+}
+
+pair min(path p)
+{
+ return p.min();
+}
+
+pair max(path p)
+{
+ return p.max();
+}
+
+realarray *mintimes(path p)
+{
+ array *V=new array(2);
+ pair z=p.mintimes();
+ (*V)[0]=z.getx();
+ (*V)[1]=z.gety();
+ return V;
+}
+
+realarray *maxtimes(path p)
+{
+ array *V=new array(2);
+ pair z=p.maxtimes();
+ (*V)[0]=z.getx();
+ (*V)[1]=z.gety();
+ return V;
+}
+
+real relativedistance(real theta, real phi, real t, bool atleast)
+{
+ return camp::velocity(theta,phi,tension(t,atleast));
+}
+
+Int windingnumber(patharray *p, pair z)
+{
+ return windingnumber(p,z);
+}
+
+bool inside(explicit patharray *g, pair z, pen fillrule=CURRENTPEN)
+{
+ return fillrule.inside(windingnumber(g,z));
+}
+
+bool inside(path g, pair z, pen fillrule=CURRENTPEN)
+{
+ return fillrule.inside(g.windingnumber(z));
+}
+
+// Determine the side of a--b that c lies on
+// (negative=left, zero=on line, positive=right).
+real side(pair a, pair b, pair c)
+{
+ return orient2d(a,b,c);
+}
+
+// Determine the side of the counterclockwise circle through a,b,c that d
+// lies on (negative=inside, 0=on circle, positive=right).
+real incircle(pair a, pair b, pair c, pair d)
+{
+ return incircle(a.getx(),a.gety(),b.getx(),b.gety(),c.getx(),c.gety(),
+ d.getx(),d.gety());
+}
+
+
+// Path3 operations
+
+path3 path3(triplearray *pre, triplearray *point, triplearray *post,
+ boolarray *straight, bool cyclic)
+{
+ size_t n=checkArrays(pre,point);
+ checkEqual(n,checkArray(post));
+ checkEqual(n,checkArray(straight));
+ mem::vector<solvedKnot3> nodes(n);
+ for(size_t i=0; i < n; ++i) {
+ nodes[i].pre=read<triple>(pre,i);
+ nodes[i].point=read<triple>(point,i);
+ nodes[i].post=read<triple>(post,i);
+ nodes[i].straight=read<bool>(straight,i);
+ }
+
+ return path3(nodes,(Int) n,cyclic);
+}
+
+path3 :nullPath3()
+{
+ return nullpath3;
+}
+
+bool ==(path3 a, path3 b)
+{
+ return a == b;
+}
+
+bool !=(path3 a, path3 b)
+{
+ return !(a == b);
+}
+
+triple point(path3 p, Int t)
+{
+ return p.point((Int) t);
+}
+
+triple point(path3 p, real t)
+{
+ return p.point(t);
+}
+
+triple precontrol(path3 p, Int t)
+{
+ return p.precontrol((Int) t);
+}
+
+triple precontrol(path3 p, real t)
+{
+ return p.precontrol(t);
+}
+
+triple postcontrol(path3 p, Int t)
+{
+ return p.postcontrol((Int) t);
+}
+
+triple postcontrol(path3 p, real t)
+{
+ return p.postcontrol(t);
+}
+
+triple dir(path3 p, Int t, Int sign=0, bool normalize=true)
+{
+ return p.dir(t,sign,normalize);
+}
+
+triple dir(path3 p, real t, bool normalize=true)
+{
+ return p.dir(t,normalize);
+}
+
+triple accel(path3 p, Int t, Int sign=0)
+{
+ return p.accel(t,sign);
+}
+
+triple accel(path3 p, real t)
+{
+ return p.accel(t);
+}
+
+real radius(path3 p, real t)
+{
+ triple v=p.dir(t,false);
+ triple a=p.accel(t);
+ real d=dot(a,v);
+ real v2=v.abs2();
+ real a2=a.abs2();
+ real denom=v2*a2-d*d;
+ real r=v2*sqrt(v2);
+ return denom > 0 ? r/sqrt(denom) : 0.0;
+}
+
+real radius(triple z0, triple c0, triple c1, triple z1, real t)
+{
+ triple v=(3.0*(z1-z0)+9.0*(c0-c1))*t*t+(6.0*(z0+c1)-12.0*c0)*t+3.0*(c0-z0);
+ triple a=6.0*(z1-z0+3.0*(c0-c1))*t+6.0*(z0+c1)-12.0*c0;
+ real d=dot(a,v);
+ real v2=v.abs2();
+ real a2=a.abs2();
+ real denom=v2*a2-d*d;
+ real r=v2*sqrt(v2);
+ return denom > 0 ? r/sqrt(denom) : 0.0;
+}
+
+path3 reverse(path3 p)
+{
+ return p.reverse();
+}
+
+path3 subpath(path3 p, Int a, Int b)
+{
+ return p.subpath((Int) a, (Int) b);
+}
+
+path3 subpath(path3 p, real a, real b)
+{
+ return p.subpath(a,b);
+}
+
+Int length(path3 p)
+{
+ return p.length();
+}
+
+bool cyclic(path3 p)
+{
+ return p.cyclic();
+}
+
+bool straight(path3 p, Int t)
+{
+ return p.straight(t);
+}
+
+// Return the component of vector v perpendicular to a unit vector u.
+triple perp(triple v, triple u)
+{
+ return v-dot(v,u)*u;
+}
+
+// Return the maximum perpendicular deviation of segment i of path3 g
+// from a straight line.
+real straightness(path3 p, Int t)
+{
+ if(p.straight(t)) return 0;
+ triple z0=p.point(t);
+ triple u=unit(p.point(t+1)-z0);
+ return ::max(length(perp(p.postcontrol(t)-z0,u)),
+ length(perp(p.precontrol(t+1)-z0,u)));
+}
+
+// Return the maximum perpendicular deviation of z0..controls c0 and c1..z1
+// from a straight line.
+real straightness(triple z0, triple c0, triple c1, triple z1)
+{
+ triple u=unit(z1-z0);
+ return ::max(length(perp(c0-z0,u)),length(perp(c1-z0,u)));
+}
+
+bool piecewisestraight(path3 p)
+{
+ return p.piecewisestraight();
+}
+
+real arclength(path3 p)
+{
+ return p.arclength();
+}
+
+real arctime(path3 p, real dval)
+{
+ return p.arctime(dval);
+}
+
+realarray* intersect(path3 p, path3 q, real fuzz=-1)
+{
+ bool exact=fuzz <= 0.0;
+ if(fuzz < 0)
+ fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())),
+ ::max(length(q.max()),length(q.min())));
+
+ std::vector<real> S,T;
+ real s,t;
+ if(intersections(s,t,S,T,p,q,fuzz,true,exact)) {
+ array *V=new array(2);
+ (*V)[0]=s;
+ (*V)[1]=t;
+ return V;
+ } else
+ return new array(0);
+}
+
+realarray2* intersections(path3 p, path3 q, real fuzz=-1)
+{
+ bool exact=fuzz <= 0.0;
+ if(fuzz < 0)
+ fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())),
+ ::max(length(q.max()),length(q.min())));
+ bool single=!exact;
+
+ real s,t;
+ std::vector<real> S,T;
+ bool found=intersections(s,t,S,T,p,q,fuzz,single,exact);
+ if(!found) return new array(0);
+ array *V;
+ if(single) {
+ V=new array(1);
+ array *Vi=new array(2);
+ (*V)[0]=Vi;
+ (*Vi)[0]=s;
+ (*Vi)[1]=t;
+ } else {
+ size_t n=S.size();
+ V=new array(n);
+ for(size_t i=0; i < n; ++i) {
+ array *Vi=new array(2);
+ (*V)[i]=Vi;
+ (*Vi)[0]=S[i];
+ (*Vi)[1]=T[i];
+ }
+ }
+ stable_sort(V->begin(),V->end(),run::compare2<real>());
+ return V;
+}
+
+Int size(path3 p)
+{
+ return p.size();
+}
+
+path3 &(path3 p, path3 q)
+{
+ return camp::concat(p,q);
+}
+
+triple min(path3 p)
+{
+ return p.min();
+}
+
+triple max(path3 p)
+{
+ return p.max();
+}
+
+realarray *mintimes(path3 p)
+{
+ array *V=new array(3);
+ triple v=p.mintimes();
+ (*V)[0]=v.getx();
+ (*V)[1]=v.gety();
+ (*V)[2]=v.getz();
+ return V;
+}
+
+realarray *maxtimes(path3 p)
+{
+ array *V=new array(3);
+ triple v=p.maxtimes();
+ (*V)[0]=v.getx();
+ (*V)[1]=v.gety();
+ (*V)[2]=v.getz();
+ return V;
+}
+
+path3 Operator *(realarray2 *t, path3 g)
+{
+ return transformed(*t,g);
+}
+
+
+// Guide operations
+
+guide* :nullGuide()
+{
+ return new pathguide(path());
+}
+
+guide* :dotsGuide(guidearray *a)
+{
+ guidevector v;
+ size_t size=checkArray(a);
+ for (size_t i=0; i < size; ++i)
+ v.push_back(a->read<guide*>(i));
+
+ return new multiguide(v);
+}
+
+guide* :dashesGuide(guidearray *a)
+{
+ static camp::curlSpec curly;
+ static specguide curlout(&curly, camp::OUT);
+ static specguide curlin(&curly, camp::IN);
+
+ size_t n=checkArray(a);
+
+ // a--b is equivalent to a{curl 1}..{curl 1}b
+ guidevector v;
+ if (n > 0)
+ v.push_back(a->read<guide*>(0));
+
+ if (n==1) {
+ v.push_back(&curlout);
+ v.push_back(&curlin);
+ }
+ else
+ for (size_t i=1; i<n; ++i) {
+ v.push_back(&curlout);
+ v.push_back(&curlin);
+ v.push_back(a->read<guide*>(i));
+ }
+
+ return new multiguide(v);
+}
+
+cycleToken :newCycleToken()
+{
+ return cycleToken();
+}
+
+guide *operator cast(cycleToken tok)
+{
+// Avoid unused variable warning messages.
+ unused(&tok);
+ return new cycletokguide();
+}
+
+guide* operator spec(pair z, Int p)
+{
+ camp::side d=(camp::side) p;
+ camp::dirSpec *sp=new camp::dirSpec(z);
+
+ return new specguide(sp,d);
+}
+
+curlSpecifier operator curl(real gamma, Int p)
+{
+ camp::side s=(camp::side) p;
+ return curlSpecifier(gamma,s);
+}
+
+real :curlSpecifierValuePart(curlSpecifier spec)
+{
+ return spec.getValue();
+}
+
+Int :curlSpecifierSidePart(curlSpecifier spec)
+{
+ return spec.getSide();
+}
+
+guide *operator cast(curlSpecifier spec)
+{
+ return new specguide(spec);
+}
+
+tensionSpecifier operator tension(real tout, real tin, bool atleast)
+{
+ return tensionSpecifier(tout, tin, atleast);
+}
+
+real :tensionSpecifierOutPart(tensionSpecifier t)
+{
+ return t.getOut();
+}
+
+real :tensionSpecifierInPart(tensionSpecifier t)
+{
+ return t.getIn();
+}
+
+bool :tensionSpecifierAtleastPart(tensionSpecifier t)
+{
+ return t.getAtleast();
+}
+
+guide *operator cast(tensionSpecifier t)
+{
+ return new tensionguide(t);
+}
+
+guide* operator controls(pair zout, pair zin)
+{
+ return new controlguide(zout, zin);
+}
+
+Int size(guide *g)
+{
+ flatguide f;
+ g->flatten(f,false);
+ return f.size();
+}
+
+Int length(guide *g)
+{
+ flatguide f;
+ g->flatten(f,false);
+ return g->cyclic() ? f.size() : f.size()-1;
+}
+
+bool cyclic(guide *g)
+{
+ flatguide f;
+ g->flatten(f,false);
+ return g->cyclic();
+}
+
+pair point(guide *g, Int t)
+{
+ flatguide f;
+ g->flatten(f,false);
+ return f.Nodes(adjustedIndex(t,f.size(),g->cyclic())).z;
+}
+
+pairarray *dirSpecifier(guide *g, Int t)
+{
+ flatguide f;
+ g->flatten(f,false);
+ Int n=f.size();
+ if(!g->cyclic() && (t < 0 || t >= n-1)) return new array(0);
+ array *c=new array(2);
+ (*c)[0]=f.Nodes(t).out->dir();
+ (*c)[1]=f.Nodes(t+1).in->dir();
+ return c;
+}
+
+pairarray *controlSpecifier(guide *g, Int t)
+{
+ flatguide f;
+ g->flatten(f,false);
+ Int n=f.size();
+ if(!g->cyclic() && (t < 0 || t >= n-1)) return new array(0);
+ knot curr=f.Nodes(t);
+ knot next=f.Nodes(t+1);
+ if(curr.out->controlled()) {
+ assert(next.in->controlled());
+ array *c=new array(2);
+ (*c)[0]=curr.out->control();
+ (*c)[1]=next.in->control();
+ return c;
+ } else return new array(0);
+}
+
+tensionSpecifier tensionSpecifier(guide *g, Int t)
+{
+ flatguide f;
+ g->flatten(f,false);
+ Int n=f.size();
+ if(!g->cyclic() && (t < 0 || t >= n-1)) return tensionSpecifier(1.0,1.0,false);
+ knot curr=f.Nodes(t);
+ return tensionSpecifier(curr.tout.val,f.Nodes(t+1).tin.val,curr.tout.atleast);
+}
+
+realarray *curlSpecifier(guide *g, Int t)
+{
+ flatguide f;
+ g->flatten(f,false);
+ Int n=f.size();
+ if(!g->cyclic() && (t < 0 || t >= n-1)) return new array(0);
+ array *c=new array(2);
+ real c0=f.Nodes(t).out->curl();
+ real c1=f.Nodes(t+1).in->curl();
+ (*c)[0]=c0 >= 0.0 ? c0 : 1.0;
+ (*c)[1]=c1 >= 0.0 ? c1 : 1.0;
+ return c;
+}
+
+guide *reverse(guide *g)
+{
+ flatguide f;
+ g->flatten(f,false);
+ if(f.precyclic())
+ return new pathguide(g->solve().reverse());
+
+ size_t n=f.size();
+ bool cyclic=g->cyclic();
+ guidevector v;
+ if(n >= 0) {
+ size_t start=cyclic ? n : n-1;
+ knot curr=f.Nodes(start);
+ knot next;
+ for(size_t i=start; i > 0; --i) {
+ next=f.Nodes(i-1);
+ v.push_back(new pairguide(curr.z));
+ if(next.out->controlled()) {
+ assert(curr.in->controlled());
+ v.push_back(new controlguide(curr.in->control(),next.out->control()));
+ } else {
+ pair d=curr.in->dir();
+ if(d != zero)
+ v.push_back(new specguide(new dirSpec(-d),camp::OUT));
+ else {
+ real C=curr.in->curl();
+ if(C >= 0.0)
+ v.push_back(new specguide(new curlSpec(C),camp::OUT));
+ }
+ real tout=curr.tin.val;
+ real tin=next.tout.val;
+ bool atleast=next.tout.atleast;
+ if(tout != 1.0 || tin != 1.0 || next.tout.atleast)
+ v.push_back(new tensionguide(tensionSpecifier(tout,tin,atleast)));
+ d=next.out->dir();
+ if(d != zero)
+ v.push_back(new specguide(new dirSpec(-d),camp::IN));
+ else {
+ real C=next.out->curl();
+ if(C >= 0.0)
+ v.push_back(new specguide(new curlSpec(C),camp::IN));
+ }
+ }
+ curr=next;
+ }
+ if(cyclic)
+ v.push_back(new cycletokguide());
+ else
+ v.push_back(new pairguide(next.z));
+
+ }
+ return new multiguide(v);
+}
+
+
+// Three-dimensional picture and surface operations
+
+void _draw(picture *f, path3 g, pen p)
+{
+ f->append(new drawPath3(g,p));
+}
+
+void draw(picture *f, triplearray2 *g, bool straight, penarray *p, real opacity,
+ real shininess, real granularity, triple normal, bool lighton,
+ penarray *colors)
+{
+ f->append(new drawSurface(*g,straight,*p,opacity,shininess,granularity,
+ normal,lighton,*colors));
+}
+
+triple min3(picture *f)
+{
+ return f->bounds3().Min();
+}
+
+triple max3(picture *f)
+{
+ return f->bounds3().Max();
+}
+
+pair min(picture *f, realarray2 *t)
+{
+ real *T=copyArray2C(t,4);
+ pair m=f->bounds(::min,xproject,yproject,T);
+ delete[] T;
+ return m;
+}
+
+pair max(picture *f, realarray2 *t)
+{
+ real *T=copyArray2C(t,4);
+ pair M=f->bounds(::max,xproject,yproject,T);
+ delete[] T;
+ return M;
+}
+
+pair minratio(picture *f)
+{
+ return f->bounds(::min,xratio,yratio);
+}
+
+pair maxratio(picture *f)
+{
+ return f->bounds(::max,xratio,yratio);
+}
+
+real minbound(realarray *a, real b)
+{
+ real *A=copyArrayC(a,16);
+ b=bound(A,::min,b);
+ delete[] A;
+ return b;
+}
+
+real maxbound(realarray *a, real b)
+{
+ real *A=copyArrayC(a,16);
+ b=bound(A,::max,b);
+ delete[] A;
+ return b;
+}
+
+pair minbound(triplearray *a, realarray2 *t, pair b)
+{
+ triple *A=copyTripleArrayC(a,16);
+ real *T=copyArray2C(t,4);
+ b=pair(bound(A,::min,xproject,T,b.getx()),
+ bound(A,::min,yproject,T,b.gety()));
+ delete[] T;
+ delete[] A;
+ return b;
+}
+
+pair maxbound(triplearray *a, realarray2 *t, pair b)
+{
+ triple *A=copyTripleArrayC(a,16);
+ real *T=copyArray2C(t,4);
+ b=pair(bound(A,::max,xproject,T,b.getx()),
+ bound(A,::max,yproject,T,b.gety()));
+ delete[] T;
+ delete[] A;
+ return b;
+}
+
+pair max(path3 g, realarray2 *t)
+{
+ real *T=copyArray2C(t,4);
+ pair b=g.bounds(::max,xproject,yproject,T);
+ delete[] T;
+ return b;
+}
+
+pair min(path3 g, realarray2 *t)
+{
+ real *T=copyArray2C(t,4);
+ pair b=g.bounds(::min,xproject,yproject,T);
+ delete[] T;
+ return b;
+}
+
+real change2(triplearray2 *a)
+{
+ size_t n=checkArray(a);
+ if(n == 0) return 0.0;
+
+ vm::array *a0=vm::read<vm::array*>(a,0);
+ size_t m=checkArray(a0);
+ if(m == 0) return 0.0;
+ triple a00=vm::read<triple>(a0,0);
+ real M=0.0;
+
+ for(size_t i=0; i < n; ++i) {
+ vm::array *ai=vm::read<vm::array*>(a,i);
+ size_t m=checkArray(ai);
+ for(size_t j=0; j < m; ++j) {
+ real a=(vm::read<triple>(ai,j)-a00).abs2();
+ if(a > M) M=a;
+ }
+ }
+ return M;
+}
+
+bool is3D(picture *f)
+{
+ return f->have3D();
+}
+
+pair bezier(pair a, pair b, pair c, pair d, real t)
+{
+ real onemt=1-t;
+ real onemt2=onemt*onemt;
+ return onemt2*onemt*a+t*(3.0*(onemt2*b+t*onemt*c)+t*t*d);
+}
+
+pair bezierP(pair a, pair b, pair c, pair d, real t)
+{
+ return 3.0*(t*t*(d-a+3.0*(b-c))+t*(2.0*(a+c)-4.0*b)+b-a);
+}
+
+pair bezierPP(pair a, pair b, pair c, pair d, real t)
+{
+ return 6.0*(t*(d-a+3.0*(b-c))+a+c-2.0*b);
+}
+
+pair bezierPPP(pair a, pair b, pair c, pair d)
+{
+ return 6.0*(d-a+3.0*(b-c));
+}
+
+triple bezier(triple a, triple b, triple c, triple d, real t)
+{
+ real onemt=1-t;
+ real onemt2=onemt*onemt;
+ return onemt2*onemt*a+t*(3.0*(onemt2*b+t*onemt*c)+t*t*d);
+}
+
+triple bezierP(triple a, triple b, triple c, triple d, real t)
+{
+ return 3.0*(t*t*(d-a+3.0*(b-c))+t*(2.0*(a+c)-4.0*b)+b-a);
+}
+
+triple bezierPP(triple a, triple b, triple c, triple d, real t)
+{
+ return 6.0*(t*(d-a+3.0*(b-c))+a+c-2.0*b);
+}
+
+triple bezierPPP(triple a, triple b, triple c, triple d)
+{
+ return 6.0*(d-a+3.0*(b-c));
+}
+
+
+// String operations
+
+string :emptyString()
+{
+ return emptystring;
+}
+
+Int length(string *s)
+{
+ return (Int) s->length();
+}
+
+Int find(string *s, string t, Int pos=0)
+{
+ size_t n=s->find(t,pos);
+ return n == string::npos ? (Int) -1 : (Int) n;
+}
+
+Int rfind(string *s, string t, Int pos=-1)
+{
+ size_t n=s->rfind(t,pos);
+ return n == string::npos ? (Int) -1 : (Int) n;
+}
+
+string reverse(string s)
+{
+ reverse(s.begin(),s.end());
+ return s;
+}
+
+string insert(string s, Int pos, string t)
+{
+ if ((size_t) pos < s.length())
+ return s.insert(pos,t);
+ return s;
+}
+
+string substr(string* s, Int pos, Int n=-1)
+{
+ if ((size_t) pos < s->length())
+ return s->substr(pos,n);
+ return emptystring;
+}
+
+string erase(string s, Int pos, Int n)
+{
+ if ((size_t) pos < s.length())
+ return s.erase(pos,n);
+ return s;
+}
+
+string downcase(string s)
+{
+ std::transform(s.begin(),s.end(),s.begin(),tolower);
+ return s;
+}
+
+string upcase(string s)
+{
+ std::transform(s.begin(),s.end(),s.begin(),toupper);
+ return s;
+}
+
+// returns a string constructed by translating all occurrences of the string
+// from in an array of string pairs {from,to} to the string to in string s.
+string replace(string *S, stringarray2 *translate)
+{
+ size_t size=checkArray(translate);
+ for(size_t i=0; i < size; i++) {
+ array *a=read<array*>(translate,i);
+ checkArray(a);
+ }
+ const char *p=S->c_str();
+ ostringstream buf;
+ while(*p) {
+ for(size_t i=0; i < size;) {
+ array *a=read<array*>(translate,i);
+ string* from=read<string*>(a,0);
+ size_t len=from->length();
+ if(strncmp(p,from->c_str(),len) != 0) {i++; continue;}
+ buf << read<string>(a,1);
+ p += len;
+ if(*p == 0) return buf.str();
+ i=0;
+ }
+ buf << *(p++);
+ }
+ return buf.str();
+}
+
+string format(string *format, Int x)
+{
+ const char *f=format->c_str();
+
+ checkformat(f,true);
+
+ Int size=snprintf(NULL,0,f,x)+1;
+ if(size < 1) size=255; // Workaround for non-C99 compliant systems.
+ char *buf=new char[size];
+ snprintf(buf,size,f,x);
+ string s=string(buf);
+ delete[] buf;
+ return s;
+}
+
+string format(string *format, real x, string locale=emptystring)
+{
+ ostringstream out;
+
+ checkformat(format->c_str(),false);
+
+ const char *phantom="\\phantom{+}";
+ const char *p0=format->c_str();
+
+ const char *p=p0;
+ const char *start=NULL;
+ while (*p != 0) {
+ if(*p == '%') {
+ p++;
+ if(*p != '%') {start=p-1; break;}
+ }
+ out << *(p++);
+ }
+
+ if(!start) return out.str();
+
+ // Allow at most 1 argument
+ while (*p != 0) {
+ if(*p == '*' || *p == '$') return out.str();
+ if(isupper(*p) || islower(*p)) {p++; break;}
+ p++;
+ }
+
+ const char *tail=p;
+ string f=format->substr(start-p0,tail-start);
+
+ const char *oldlocale=NULL;
+ if(!locale.empty()) {
+ oldlocale=setlocale(LC_ALL,NULL);
+ if(oldlocale) oldlocale=StrdupNoGC(oldlocale);
+ setlocale(LC_ALL,locale.c_str());
+ }
+
+ Int size=snprintf(NULL,0,f.c_str(),x)+1;
+ if(size < 1) size=255; // Workaround for non-C99 compliant systems.
+ char *buf=new char[size];
+ snprintf(buf,size,f.c_str(),x);
+
+ if(oldlocale) {
+ setlocale(LC_ALL,oldlocale);
+ delete[] oldlocale;
+ }
+
+ bool trailingzero=f.find("#") < string::npos;
+ bool plus=f.find("+") < string::npos;
+ bool space=f.find(" ") < string::npos;
+
+ char *q=buf; // beginning of formatted number
+
+ if(*q == ' ') {
+ out << phantom;
+ q++;
+ }
+
+ const char decimal=*(localeconv()->decimal_point);
+
+ // Remove any spurious sign
+ if(*q == '-' || *q == '+') {
+ p=q+1;
+ bool zero=true;
+ while(*p != 0) {
+ if(!isdigit(*p) && *p != decimal) break;
+ if(isdigit(*p) && *p != '0') {zero=false; break;}
+ p++;
+ }
+ if(zero) {
+ q++;
+ if(plus || space) out << phantom;
+ }
+ }
+
+ const char *r=p=q;
+ bool dp=false;
+ while(*r != 0 && (isdigit(*r) || *r == decimal || *r == '+' || *r == '-')) {
+ if(*r == decimal) dp=true;
+ r++;
+ }
+ if(dp) { // Remove trailing zeros and/or decimal point
+ r--;
+ unsigned n=0;
+ while(r > q && *r == '0') {r--; n++;}
+ if(*r == decimal) {r--; n++;}
+ while(q <= r) out << *(q++);
+ if(!trailingzero) q += n;
+ }
+
+ bool zero=(r == p && *r == '0') && !trailingzero;
+
+ // Translate "E+/E-/e+/e-" exponential notation to TeX
+ while(*q != 0) {
+ if((*q == 'E' || *q == 'e') && (*(q+1) == '+' || *(q+1) == '-')) {
+ if(!zero) out << "\\!\\times\\!10^{";
+ bool plus=(*(q+1) == '+');
+ q++;
+ if(plus) q++;
+ if(*q == '-') out << *(q++);
+ while(*q == '0' && (zero || isdigit(*(q+1)))) q++;
+ while(isdigit(*q)) out << *(q++);
+ if(!zero)
+ out << "}";
+ break;
+ }
+ out << *(q++);
+ }
+
+ while(*tail != 0)
+ out << *(tail++);
+
+ delete[] buf;
+
+ return out.str();
+}
+
+Int hex(string s)
+{
+ istringstream is(s);
+ is.setf(std::ios::hex,std::ios::basefield);
+ Int value;
+ if(is && is >> value && ((is >> std::ws).eof())) return value;
+ ostringstream buf;
+ buf << "invalid hexidecimal cast from string \"" << s << "\"";
+ error(buf);
+}
+
+string string(Int x)
+{
+ ostringstream buf;
+ buf << x;
+ return buf.str();
+}
+
+string string(real x, Int digits=DBL_DIG)
+{
+ ostringstream buf;
+ buf.precision(digits);
+ buf << x;
+ return buf.str();
+}
+
+string time(string format=defaulttimeformat)
+{
+#ifdef HAVE_STRFTIME
+ const time_t bintime=time(NULL);
+ if(!strftime(Time,nTime,format.c_str(),localtime(&bintime))) return "";
+ return Time;
+#else
+ return format;
+#endif
+}
+
+string time(Int seconds, string format=defaulttimeformat)
+{
+#ifdef HAVE_STRFTIME
+ const time_t bintime=seconds;
+ if(!strftime(Time,nTime,format.c_str(),localtime(&bintime))) return "";
+ return Time;
+#else
+// Avoid unused variable warning messages
+ unused(&seconds);
+ return format;
+#endif
+}
+
+Int seconds(string t=emptystring, string format=emptystring)
+{
+#if defined(HAVE_STRPTIME)
+ const time_t bintime=time(NULL);
+ tm tm=*localtime(&bintime);
+ if(t != "" && !strptime(t.c_str(),format.c_str(),&tm)) return -1;
+ return (Int) mktime(&tm);
+#else
+ return -1;
+#endif
+}
+
+realarray *_cputime()
+{
+ static const real ticktime=1.0/sysconf(_SC_CLK_TCK);
+ struct tms buf;
+
+ ::times(&buf);
+ array *t=new array(4);
+ (*t)[0] = ((real) buf.tms_utime)*ticktime;
+ (*t)[1] = ((real) buf.tms_stime)*ticktime;
+ (*t)[2] = ((real) buf.tms_cutime)*ticktime;
+ (*t)[3] = ((real) buf.tms_cstime)*ticktime;
+ return t;
+}
+
+
+// Math
+
+real ^(real x, Int y)
+{
+ return pow(x,y);
+}
+
+pair ^(pair z, Int y)
+{
+ return pow(z,y);
+}
+
+Int quotient(Int x, Int y)
+{
+ if(y == 0) dividebyzero();
+ if(y == -1) return Negate(x);
+// Implementation-independent definition of integer division: round down
+ return (x-portableMod(x,y))/y;
+}
+
+Int abs(Int x)
+{
+ return Abs(x);
+}
+
+Int sgn(real x)
+{
+ return sgn(x);
+}
+
+Int rand()
+{
+ return rand();
+}
+
+void srand(Int seed)
+{
+ srand(intcast(seed));
+}
+
+// a random number uniformly distributed in the interval [0,1]
+real unitrand()
+{
+ return ((real) rand())/RAND_MAX;
+}
+
+Int ceil(real x)
+{
+ return Intcast(ceil(x));
+}
+
+Int floor(real x)
+{
+ return Intcast(floor(x));
+}
+
+Int round(real x)
+{
+ if(validInt(x)) return Round(x);
+ integeroverflow(0);
+}
+
+Int Ceil(real x)
+{
+ return Ceil(x);
+}
+
+Int Floor(real x)
+{
+ return Floor(x);
+}
+
+Int Round(real x)
+{
+ return Round(Intcap(x));
+}
+
+real fmod(real x, real y)
+{
+ if (y == 0.0) dividebyzero();
+ return fmod(x,y);
+}
+
+real atan2(real y, real x)
+{
+ return atan2(y,x);
+}
+
+real hypot(real x, real y)
+{
+ return hypot(x,y);
+}
+
+real remainder(real x, real y)
+{
+ return remainder(x,y);
+}
+
+real J(Int n, real x)
+{
+ return jn(n,x);
+}
+
+real Y(Int n, real x)
+{
+ return yn(n,x);
+}
+
+real erf(real x)
+{
+ return erf(x);
+}
+
+real erfc(real x)
+{
+ return erfc(x);
+}
+
+Int factorial(Int n) {
+ if(n < 0) error(invalidargument);
+ return factorial(n);
+}
+
+Int choose(Int n, Int k) {
+ if(n < 0 || k < 0 || k > n) error(invalidargument);
+ Int f=1;
+ Int r=n-k;
+ for(Int i=n; i > r; --i) {
+ if(f > Int_MAX/i) integeroverflow(0);
+ f=(f*i)/(n-i+1);
+ }
+ return f;
+}
+
+real gamma(real x)
+{
+#ifdef HAVE_TGAMMA
+ return tgamma(x);
+#else
+ real lg = lgamma(x);
+ return signgam*exp(lg);
+#endif
+}
+
+// Complex Gamma function
+pair gamma(explicit pair z)
+{
+ return gamma(z);
+}
+
+realarray *quadraticroots(real a, real b, real c)
+{
+ quadraticroots q(a,b,c);
+ array *roots=new array(q.roots);
+ if(q.roots >= 1) (*roots)[0]=q.t1;
+ if(q.roots == 2) (*roots)[1]=q.t2;
+ return roots;
+}
+
+pairarray *quadraticroots(explicit pair a, explicit pair b, explicit pair c)
+{
+ Quadraticroots q(a,b,c);
+ array *roots=new array(q.roots);
+ if(q.roots >= 1) (*roots)[0]=q.z1;
+ if(q.roots == 2) (*roots)[1]=q.z2;
+ return roots;
+}
+
+realarray *cubicroots(real a, real b, real c, real d)
+{
+ cubicroots q(a,b,c,d);
+ array *roots=new array(q.roots);
+ if(q.roots >= 1) (*roots)[0]=q.t1;
+ if(q.roots >= 2) (*roots)[1]=q.t2;
+ if(q.roots == 3) (*roots)[2]=q.t3;
+ return roots;
+}
+
+
+// Transforms
+
+bool ==(transform a, transform b)
+{
+ return a == b;
+}
+
+bool !=(transform a, transform b)
+{
+ return a != b;
+}
+
+transform +(transform a, transform b)
+{
+ return a+b;
+}
+
+transform Operator *(transform a, transform b)
+{
+ return a*b;
+}
+
+pair Operator *(transform t, pair z)
+{
+ return t*z;
+}
+
+path Operator *(transform t, path g)
+{
+ return transformed(t,g);
+}
+
+pen Operator *(transform t, pen p)
+{
+ return transformed(t,p);
+}
+
+picture * Operator *(transform t, picture *f)
+{
+ return transformed(t,f);
+}
+
+picture * Operator *(realarray2 *t, picture *f)
+{
+ return transformed(*t,f);
+}
+
+transform ^(transform t, Int n)
+{
+ transform T;
+ if(n < 0) {
+ n=-n;
+ t=inverse(t);
+ }
+ for(Int i=0; i < n; i++) T=T*t;
+ return T;
+}
+
+real :transformXPart(transform t)
+{
+ return t.getx();
+}
+
+real :transformYPart(transform t)
+{
+ return t.gety();
+}
+
+real :transformXXPart(transform t)
+{
+ return t.getxx();
+}
+
+real :transformXYPart(transform t)
+{
+ return t.getxy();
+}
+
+real :transformYXPart(transform t)
+{
+ return t.getyx();
+}
+
+real :transformYYPart(transform t)
+{
+ return t.getyy();
+}
+
+transform :real6ToTransform(real x, real y, real xx, real xy,
+ real yx, real yy)
+{
+ return transform(x,y,xx,xy,yx,yy);
+}
+
+transform shift(transform t)
+{
+ return transform(t.getx(),t.gety(),0,0,0,0);
+}
+
+transform shiftless(transform t)
+{
+ return transform(0,0,t.getxx(),t.getxy(),t.getyx(),t.getyy());
+}
+
+transform identity:transformIdentity()
+{
+ return identity;
+}
+
+transform inverse(transform t)
+{
+ return inverse(t);
+}
+
+transform shift(pair z)
+{
+ return shift(z);
+}
+
+transform shift(real x, real y)
+{
+ return shift(pair(x,y));
+}
+
+transform xscale(real x)
+{
+ return xscale(x);
+}
+
+transform yscale(real y)
+{
+ return yscale(y);
+}
+
+transform scale(real x)
+{
+ return scale(x);
+}
+
+transform scale(real x, real y)
+{
+ return xscale(x)*yscale(y);
+}
+
+transform slant(real s)
+{
+ return slant(s);
+}
+
+transform rotate(real angle, pair z=0)
+{
+ return rotatearound(z,radians(angle));
+}
+
+transform reflect(pair a, pair b)
+{
+ return reflectabout(a,b);
+}
+
+
+// Pair operations
+
+pair :pairZero()
+{
+ return zero;
+}
+
+pair :realRealToPair(real x, real y)
+{
+ return pair(x,y);
+}
+
+pair :pairNegate(pair z)
+{
+ return -z;
+}
+
+real xpart:pairXPart(pair z)
+{
+ return z.getx();
+}
+
+real ypart:pairYPart(pair z)
+{
+ return z.gety();
+}
+
+real length(pair z)
+{
+ return z.length();
+}
+
+real abs(pair z)
+{
+ return z.length();
+}
+
+pair sqrt(explicit pair z)
+{
+ return Sqrt(z);
+}
+
+// Return the angle of z in radians.
+real angle(pair z, bool warn=true)
+{
+ if(!warn && z.getx() == 0.0 && z.gety() == 0.0) return 0.0;
+ return z.angle();
+}
+
+// Return the angle of z in degrees in the interval [0,360).
+real degrees(pair z, bool warn=true)
+{
+ if(!warn && z.getx() == 0.0 && z.gety() == 0.0) return 0.0;
+ return principalBranch(degrees(z.angle()));
+}
+
+// Convert degrees to radians.
+real radians(real degrees)
+{
+ return radians(degrees);
+}
+
+// Convert radians to degrees.
+real degrees(real radians)
+{
+ return degrees(radians);
+}
+
+// Convert radians to degrees in [0,360).
+real Degrees(real radians)
+{
+ return principalBranch(degrees(radians));
+}
+
+real Sin(real deg)
+{
+ return sin(radians(deg));
+}
+
+real Cos(real deg)
+{
+ return cos(radians(deg));
+}
+
+real Tan(real deg)
+{
+ return tan(radians(deg));
+}
+
+real aSin(real x)
+{
+ return degrees(asin(x));
+}
+
+real aCos(real x)
+{
+ return degrees(acos(x));
+}
+
+real aTan(real x)
+{
+ return degrees(atan(x));
+}
+
+pair unit(pair z)
+{
+ return unit(z);
+}
+
+pair dir(real degrees)
+{
+ return expi(radians(degrees));
+}
+
+pair dir(explicit pair z)
+{
+ return unit(z);
+}
+
+pair expi(real angle)
+{
+ return expi(angle);
+}
+
+pair exp(explicit pair z)
+{
+ return exp(z);
+}
+
+pair log(explicit pair z)
+{
+ return pair(log(z.length()),z.angle());
+}
+
+pair sin(explicit pair z)
+{
+ return sin(z);
+}
+
+pair cos(explicit pair z)
+{
+ return pair(cos(z.getx())*cosh(z.gety()),-sin(z.getx())*sinh(z.gety()));
+}
+
+pair conj(pair z)
+{
+ return conj(z);
+}
+
+pair realmult(pair z, pair w)
+{
+ return pair (z.getx()*w.getx(),z.gety()*w.gety());
+}
+
+triple realmult(triple u, triple v)
+{
+ return triple (u.getx()*v.getx(),u.gety()*v.gety(),u.getz()*v.getz());
+}
+
+// To avoid confusion, a dot product requires explicit pair arguments.
+real dot(explicit pair z, explicit pair w)
+{
+ return dot(z,w);
+}
+
+
+// Triple operations
+
+triple :tripleZero()
+{
+ static triple zero;
+ return zero;
+}
+
+triple :realRealRealToTriple(real x, real y, real z)
+{
+ return triple(x,y,z);
+}
+
+real xpart:tripleXPart(triple v)
+{
+ return v.getx();
+}
+
+real ypart:tripleYPart(triple v)
+{
+ return v.gety();
+}
+
+real zpart:tripleZPart(triple v)
+{
+ return v.getz();
+}
+
+triple Operator *(real x, triple v)
+{
+ return x*v;
+}
+
+triple Operator *(triple v, real x)
+{
+ return v*x;
+}
+
+triple /(triple v, real x)
+{
+ return v/x;
+}
+
+real length(triple v)
+{
+ return v.length();
+}
+
+real abs(triple v)
+{
+ return v.length();
+}
+
+real polar(triple v)
+{
+ return v.polar();
+}
+
+real azimuth(triple v)
+{
+ return v.azimuth();
+}
+
+real colatitude(triple v)
+{
+ return degrees(v.polar());
+}
+
+real latitude(triple v)
+{
+ return 90.0-degrees(v.polar());
+}
+
+// Return the longitude of v in [0,360).
+real longitude(triple v, bool warn=true)
+{
+ if(!warn && v.getx() == 0.0 && v.gety() == 0.0) return 0.0;
+ return principalBranch(degrees(v.azimuth()));
+}
+
+triple unit(triple v)
+{
+ return unit(v);
+}
+
+real dot(triple u, triple v)
+{
+ return dot(u,v);
+}
+
+triple cross(triple u, triple v)
+{
+ return cross(u,v);
+}
+
+triple expi(real polar, real azimuth)
+{
+ return expi(polar,azimuth);
+}
+
+triple dir(real colatitude, real longitude)
+{
+ return expi(radians(colatitude),radians(longitude));
+}
+
+
+// System routines
+
+void atupdate(callable *f)
+{
+ processData().atUpdateFunction=f;
+}
+
+callable *atupdate()
+{
+ return processData().atUpdateFunction;
+}
+
+void atexit(callable *f)
+{
+ processData().atExitFunction=f;
+}
+
+callable *atexit()
+{
+ return processData().atExitFunction;
+}
+
+void atbreakpoint(callableBp *f)
+{
+ processData().atBreakpointFunction=f;
+}
+
+void breakpoint(runnable *s=NULL)
+{
+ breakpoint(Stack,s);
+}
+
+string locatefile(string file)
+{
+ return locateFile(file);
+}
+
+void stop(string file, Int line, runnable *s=NULL)
+{
+ file=locateFile(file);
+ clear(file,line);
+ cout << "setting breakpoint at " << file << ": " << line << endl;
+ bplist.push_back(bpinfo(file,line,s));
+}
+
+void breakpoints()
+{
+ for(mem::list<bpinfo>::iterator p=bplist.begin(); p != bplist.end(); ++p)
+ cout << p->f.name() << ": " << p->f.line() << endl;
+}
+
+void clear(string file, Int line)
+{
+ file=locateFile(file);
+ clear(file,line,true);
+}
+
+void clear()
+{
+ bplist.clear();
+}
+
+// Strip directory from string
+string stripdirectory(string *s)
+{
+ return stripDir(*s);
+}
+
+// Strip directory from string
+string stripfile(string *s)
+{
+ return stripFile(*s);
+}
+
+// Strip file extension from string
+string stripextension(string *s)
+{
+ return stripExt(*s);
+}
+
+// Call ImageMagick convert.
+Int convert(string args=emptystring, string file=emptystring,
+ string format=emptystring)
+{
+ ostringstream cmd;
+
+ string name=convertname(file,format);
+ cmd << "'" << getSetting<string>("convert") << "' " << args
+ << " '" << name << "'";
+ bool quiet=verbose <= 1;
+ Int ret=System(cmd,quiet ? 1 : 0,true,"convert","your ImageMagick convert utility");
+
+ if(ret == 0 && verbose > 0)
+ cout << "Wrote " << ((file.empty()) ? name : file) << endl;
+
+ return ret;
+}
+
+// Call ImageMagick animate.
+Int animate(string args=emptystring, string file=emptystring,
+ string format=emptystring)
+{
+#ifndef __CYGWIN__
+ string name=convertname(file,format,false);
+ if(view()) {
+ ostringstream cmd;
+ cmd << "'" << getSetting<string>("animate") << "' " << args
+ << " '" << name << "'";
+ return System(cmd,0,false,"animate","your animated GIF viewer");
+ }
+#endif
+ return 0;
+}
+
+// Delete file named s.
+Int delete(string *s)
+{
+ checkLocal(*s);
+ Int rc=unlink(s->c_str());
+ if(rc == 0 && verbose > 0)
+ cout << "Deleted " << *s << endl;
+ return rc;
+}
+
+// Rename file "from" to file "to".
+Int rename(string *from, string *to)
+{
+ checkLocal(*from);
+ checkLocal(*to);
+ Int rc=rename(from->c_str(),to->c_str());
+ if(rc == 0 && verbose > 0)
+ cout << "Renamed " << *from << " to " << *to << endl;
+ return rc;
+}
+
+
+// Array operations
+
+// Create an empty array.
+array* :emptyArray()
+{
+ return new array(0);
+}
+
+// Create a new array (technically a vector).
+// This array will be multidimensional. First the number of dimensions
+// is popped off the stack, followed by each dimension in reverse order.
+// The array itself is technically a one dimensional array of one
+// dimension arrays and so on.
+array* :newDeepArray(Int depth)
+{
+ assert(depth > 0);
+
+ Int *dims = new Int[depth];
+
+ for (Int index = depth-1; index >= 0; index--) {
+ Int i=pop<Int>(Stack);
+ if(i < 0) error("cannot create a negative length array");
+ dims[index]=i;
+ }
+
+ array *a=deepArray(depth, dims);
+ delete[] dims;
+ return a;
+}
+
+// Creates an array with elements already specified. First, the number
+// of elements is popped off the stack, followed by each element in
+// reverse order.
+array* :newInitializedArray(Int n)
+{
+ assert(n >= 0);
+
+ array *a = new array(n);
+
+ for (Int index = n-1; index >= 0; index--)
+ (*a)[index] = pop(Stack);
+
+ return a;
+}
+
+// Similar to newInitializedArray, but after the n elements, append another
+// array to it.
+array* :newAppendedArray(array* tail, Int n)
+{
+ assert(n >= 0);
+
+ array *a = new array(n);
+
+ for (Int index = n-1; index >= 0; index--)
+ (*a)[index] = pop(Stack);
+
+ copy(tail->begin(), tail->end(), back_inserter(*a));
+
+ return a;
+}
+
+// The function T[] array(int n, T value, int depth=0) produces a array of n
+// copies of x, where each copy is copied up to depth.
+array* :newDuplicateArray(Int n, item value, Int depth=Int_MAX)
+{
+ if(n < 0) error("cannot create a negative length array");
+ if(depth < 0) error("cannot copy to a negative depth");
+
+ return new array(n, value, depth);
+}
+
+// Read an element from an array. Checks for initialization & bounds.
+item :arrayRead(array *a, Int n)
+{
+ item& i=arrayRead(a,n);
+ if (i.empty()) {
+ ostringstream buf;
+ buf << "read uninitialized value from array at index " << n;
+ error(buf);
+ }
+ return i;
+}
+
+// Slice a substring from an array.
+item :arraySliceRead(array *a, Int left, Int right)
+{
+ checkArray(a);
+ return a->slice(left, right);
+}
+
+// Slice a substring from an array. This implements the cases a[i:] and a[:]
+// where the endpoint is not given, and assumed to be the length of the array.
+item :arraySliceReadToEnd(array *a, Int left)
+{
+ size_t len=checkArray(a);
+ return a->slice(left, (Int)len);
+}
+
+// Read an element from an array of arrays. Check bounds and initialize
+// as necessary.
+item :arrayArrayRead(array *a, Int n)
+{
+ item& i=arrayRead(a,n);
+ if (i.empty()) i=new array(0);
+ return i;
+}
+
+// Write an element to an array. Increase size if necessary.
+item :arrayWrite(item value, array *a, Int n)
+{
+ size_t len=checkArray(a);
+ bool cyclic=a->cyclic();
+ if(cyclic && len > 0) n=imod(n,len);
+ else {
+ if(cyclic) outOfBounds("writing cyclic",len,n);
+ if(n < 0) outOfBounds("writing",len,n);
+ if(len <= (size_t) n)
+ a->resize(n+1);
+ }
+ (*a)[n] = value;
+ return value;
+}
+
+array * :arraySliceWrite(array *src, array *dest, Int left, Int right)
+{
+ checkArray(src);
+ checkArray(dest);
+ dest->setSlice(left, right, src);
+ return src;
+}
+
+array * :arraySliceWriteToEnd(array *src, array *dest, Int left)
+{
+ checkArray(src);
+ size_t len=checkArray(dest);
+ dest->setSlice(left, (Int) len, src);
+ return src;
+}
+
+// Returns the length of an array.
+Int :arrayLength(array *a)
+{
+ return (Int) checkArray(a);
+}
+
+// Returns an array of integers representing the keys of the array.
+array * :arrayKeys(array *a)
+{
+ size_t size=checkArray(a);
+
+ array *keys=new array();
+ for (size_t i=0; i<size; ++i) {
+ item& cell = (*a)[i];
+ if (!cell.empty())
+ keys->push((Int)i);
+ }
+
+ return keys;
+}
+
+// Return the cyclic flag for an array.
+bool :arrayCyclicFlag(array *a)
+{
+ checkArray(a);
+ return a->cyclic();
+}
+
+// Check to see if an array element is initialized.
+bool :arrayInitializedHelper(Int n, array *a)
+{
+ size_t len=checkArray(a);
+ bool cyclic=a->cyclic();
+ if(cyclic && len > 0) n=imod(n,len);
+ else if(n < 0 || n >= (Int) len) return false;
+ item&i=(*a)[(unsigned) n];
+ return !i.empty();
+}
+
+// Returns the initialize method for an array.
+callable* :arrayInitialized(array *a)
+{
+ return new thunk(new bfunc(arrayInitializedHelper),a);
+}
+
+// The helper function for the cyclic method that sets the cyclic flag.
+void :arrayCyclicHelper(bool b, array *a)
+{
+ checkArray(a);
+ a->cyclic(b);
+}
+
+// Set the cyclic flag for an array.
+callable* :arrayCyclic(array *a)
+{
+ return new thunk(new bfunc(arrayCyclicHelper),a);
+}
+
+// The helper function for the push method that does the actual operation.
+item :arrayPushHelper(item x, array *a)
+{
+ checkArray(a);
+ a->push(x);
+ return x;
+}
+
+// Returns the push method for an array.
+callable* :arrayPush(array *a)
+{
+ return new thunk(new bfunc(arrayPushHelper),a);
+}
+
+// The helper function for the append method that appends b to a.
+void :arrayAppendHelper(array *b, array *a)
+{
+ checkArray(a);
+ size_t size=checkArray(b);
+ for(size_t i=0; i < size; i++)
+ a->push((*b)[i]);
+}
+
+// Returns the append method for an array.
+callable* :arrayAppend(array *a)
+{
+ return new thunk(new bfunc(arrayAppendHelper),a);
+}
+
+// The helper function for the pop method.
+item :arrayPopHelper(array *a)
+{
+ size_t asize=checkArray(a);
+ if(asize == 0)
+ error("cannot pop element from empty array");
+ return a->pop();
+}
+
+// Returns the pop method for an array.
+callable* :arrayPop(array *a)
+{
+ return new thunk(new bfunc(arrayPopHelper),a);
+}
+
+// The helper function for the insert method.
+item :arrayInsertHelper(Int i, array *x, array *a)
+{
+ size_t asize=checkArray(a);
+ checkArray(x);
+ if(a->cyclic() && asize > 0) i=imod(i,asize);
+ if(i < 0 || i > (Int) asize)
+ outOfBounds("inserting",asize,i);
+ (*a).insert((*a).begin()+i,(*x).begin(),(*x).end());
+}
+
+// Returns the insert method for an array.
+callable* :arrayInsert(array *a)
+{
+ return new thunk(new bfunc(arrayInsertHelper),a);
+}
+
+// Returns the delete method for an array.
+callable* :arrayDelete(array *a)
+{
+ return new thunk(new bfunc(arrayDeleteHelper),a);
+}
+
+bool :arrayAlias(array *a, array *b)
+{
+ return a==b;
+}
+
+// Return array formed by indexing array a with elements of integer array b
+array* :arrayIntArray(array *a, array *b)
+{
+ size_t asize=checkArray(a);
+ size_t bsize=checkArray(b);
+ array *r=new array(bsize);
+ bool cyclic=a->cyclic();
+ for(size_t i=0; i < bsize; i++) {
+ Int index=read<Int>(b,i);
+ if(cyclic && asize > 0) index=imod(index,asize);
+ else
+ if(index < 0 || index >= (Int) asize)
+ outOfBounds("reading",asize,index);
+ (*r)[i]=(*a)[index];
+ }
+ return r;
+}
+
+// returns the complement of the integer array a in {0,2,...,n-1},
+// so that b[complement(a,b.length)] yields the complement of b[a].
+Intarray* complement(Intarray *a, Int n)
+{
+ size_t asize=checkArray(a);
+ array *r=new array(0);
+ bool *keep=new bool[n];
+ for(Int i=0; i < n; ++i) keep[i]=true;
+ for(size_t i=0; i < asize; ++i) {
+ Int j=read<Int>(a,i);
+ if(j >= 0 && j < n) keep[j]=false;
+ }
+ for(Int i=0; i < n; i++)
+ if(keep[i]) r->push(i);
+
+ delete[] keep;
+ return r;
+}
+
+// Generate the sequence {f(i) : i=0,1,...n-1} given a function f and integer n
+Intarray* :arraySequence(callable *f, Int n)
+{
+ if(n < 0) n=0;
+ array *a=new array(n);
+ for(Int i=0; i < n; ++i) {
+ Stack->push(i);
+ f->call(Stack);
+ (*a)[i]=pop(Stack);
+ }
+ return a;
+}
+
+// Return the array {0,1,...n-1}
+Intarray *sequence(Int n)
+{
+ if(n < 0) n=0;
+ array *a=new array(n);
+ for(Int i=0; i < n; ++i) {
+ (*a)[i]=i;
+ }
+ return a;
+}
+
+// Apply a function to each element of an array
+array* :arrayFunction(callable *f, array *a)
+{
+ size_t size=checkArray(a);
+ array *b=new array(size);
+ for(size_t i=0; i < size; ++i) {
+ Stack->push((*a)[i]);
+ f->call(Stack);
+ (*b)[i]=pop(Stack);
+ }
+ return b;
+}
+
+array* :arraySort(array *a, callable *f)
+{
+ array *c=copyArray(a);
+ compareFunc=f;
+ FuncStack=Stack;
+ stable_sort(c->begin(),c->end(),compareFunction);
+ return c;
+}
+
+bool all(boolarray *a)
+{
+ size_t size=checkArray(a);
+ bool c=true;
+ for(size_t i=0; i < size; i++)
+ if(!get<bool>((*a)[i])) {c=false; break;}
+ return c;
+}
+
+boolarray* !(boolarray* a)
+{
+ size_t size=checkArray(a);
+ array *c=new array(size);
+ for(size_t i=0; i < size; i++)
+ (*c)[i]=!read<bool>(a,i);
+ return c;
+}
+
+Int sum(boolarray *a)
+{
+ size_t size=checkArray(a);
+ Int sum=0;
+ for(size_t i=0; i < size; i++)
+ sum += read<bool>(a,i) ? 1 : 0;
+ return sum;
+}
+
+array* :arrayCopy(array *a)
+{
+ return copyArray(a);
+}
+
+array* :arrayConcat(array *a)
+{
+ // a is an array of arrays to be concatenated together.
+ // The signature is
+ // T[] concat(... T[][] a);
+
+ size_t numArgs=checkArray(a);
+ size_t resultSize=0;
+ for (size_t i=0; i < numArgs; ++i) {
+ resultSize += checkArray(a->read<array *>(i));
+ }
+
+ array *result=new array(resultSize);
+
+ size_t ri=0;
+ for (size_t i=0; i < numArgs; ++i) {
+ array *arg=a->read<array *>(i);
+ size_t size=checkArray(arg);
+
+ for (size_t j=0; j < size; ++j) {
+ (*result)[ri]=(*arg)[j];
+ ++ri;
+ }
+ }
+
+ return result;
+}
+
+array* :array2Copy(array *a)
+{
+ return copyArray2(a);
+}
+
+array* :array3Copy(array *a)
+{
+ return copyArray3(a);
+}
+
+array* :array2Transpose(array *a)
+{
+ size_t asize=checkArray(a);
+ array *c=new array(0);
+ for(size_t i=0; i < asize; i++) {
+ size_t ip=i+1;
+ array *ai=read<array*>(a,i);
+ size_t aisize=checkArray(ai);
+ size_t csize=checkArray(c);
+ if(csize < aisize) {
+ c->resize(aisize);
+ for(size_t j=csize; j < aisize; j++) {
+ (*c)[j]=new array(ip);
+ }
+ }
+ for(size_t j=0; j < aisize; j++) {
+ array *cj=read<array*>(c,j);
+ if(checkArray(cj) < ip) cj->resize(ip);
+ (*cj)[i]=(*ai)[j];
+ }
+ }
+ return c;
+}
+
+// a is a rectangular 3D array; perm is an Int array indicating the type of
+// permutation (021 or 120, etc; original is 012).
+// Transpose by sending respective members to the permutated locations:
+// return the array obtained by putting a[i][j][k] into position perm{ijk}.
+array* :array3Transpose(array *a, array *perm)
+{
+ const size_t DIM=3;
+
+ if(checkArray(perm) != DIM) {
+ ostringstream buf;
+ buf << "permutation array must have length " << DIM;
+ error(buf);
+ }
+
+ size_t* size=new size_t[DIM];
+ for(size_t i=0; i < DIM; ++i) size[i]=DIM;
+
+ for(size_t i=0; i < DIM; ++i) {
+ Int p=read<Int>(perm,i);
+ size_t P=(size_t) p;
+ if(p < 0 || P >= DIM) {
+ ostringstream buf;
+ buf << "permutation index out of range: " << p;
+ error(buf);
+ }
+ size[P]=P;
+ }
+
+ for(size_t i=0; i < DIM; ++i)
+ if(size[i] == DIM) error("permutation indices must be distinct");
+
+ static const char *rectangular=
+ "3D transpose implemented for rectangular matrices only";
+
+ size_t isize=size[0]=checkArray(a);
+ array *a0=read<array*>(a,0);
+ size[1]=checkArray(a0);
+ array *a00=read<array*>(a0,0);
+ size[2]=checkArray(a00);
+ for(size_t i=0; i < isize; i++) {
+ array *ai=read<array*>(a,i);
+ size_t jsize=checkArray(ai);
+ if(jsize != size[1]) error(rectangular);
+ for(size_t j=0; j < jsize; j++) {
+ array *aij=read<array*>(ai,j);
+ if(checkArray(aij) != size[2]) error(rectangular);
+ }
+ }
+
+ size_t perm0=(size_t) read<Int>(perm,0);
+ size_t perm1=(size_t) read<Int>(perm,1);
+ size_t perm2=(size_t) read<Int>(perm,2);
+
+ size_t sizep0=size[perm0];
+ size_t sizep1=size[perm1];
+ size_t sizep2=size[perm2];
+
+ array *c=new array(sizep0);
+ for(size_t i=0; i < sizep0; ++i) {
+ array *ci=new array(sizep1);
+ (*c)[i]=ci;
+ for(size_t j=0; j < sizep1; ++j) {
+ array *cij=new array(sizep2);
+ (*ci)[j]=cij;
+ }
+ }
+
+ size_t* i=new size_t[DIM];
+
+ for(i[0]=0; i[0] < size[0]; ++i[0]) {
+ array *a0=read<array*>(a,i[0]);
+ for(i[1]=0; i[1] < size[1]; ++i[1]) {
+ array *a1=read<array*>(a0,i[1]);
+ for(i[2]=0; i[2] < size[2]; ++i[2]) {
+ array *c0=read<array*>(c,i[perm0]);
+ array *c1=read<array*>(c0,i[perm1]);
+ (*c1)[i[perm2]]=read<real>(a1,i[2]);
+ }
+ }
+ }
+
+ delete [] i;
+ delete [] size;
+
+ return c;
+}
+
+// In a boolean array, find the index of the nth true value or -1 if not found
+// If n is negative, search backwards.
+Int find(boolarray *a, Int n=1)
+{
+
+ size_t size=checkArray(a);
+ Int j=-1;
+ if(n > 0)
+ for(size_t i=0; i < size; i++)
+ if(read<bool>(a,i)) {
+ n--; if(n == 0) {j=(Int) i; break;}
+ }
+ if(n < 0)
+ for(size_t i=size; i > 0;)
+ if(read<bool>(a,--i)) {
+ n++; if(n == 0) {j=(Int) i; break;}
+ }
+ return j;
+}
+
+bool Operator ==(realarray2 *a, realarray2 *b)
+{
+ size_t n=checkArray(a);
+ if(n != checkArray(b)) return false;
+
+ size_t n0=n == 0 ? 0 : checkArray(read<array*>(a,0));
+ if(n0 != checkArray(read<array*>(b,0))) return false;
+
+ for(size_t i=0; i < n; ++i) {
+ array *ai=read<array*>(a,i);
+ array *bi=read<array*>(b,i);
+ for(size_t j=0; j < n0; ++j) {
+ if(read<real>(ai,j) != read<real>(bi,j))
+ return false;
+ }
+ }
+ return true;
+}
+
+// construct vector obtained by replacing those elements of b for which the
+// corresponding elements of a are false by the corresponding element of c.
+array* :arrayConditional(array *a, array *b, array *c)
+{
+ size_t size=checkArray(a);
+ array *r=new array(size);
+ if(b && c) {
+ checkArrays(a,b);
+ checkArrays(b,c);
+ for(size_t i=0; i < size; i++)
+ (*r)[i]=read<bool>(a,i) ? (*b)[i] : (*c)[i];
+ } else {
+ r->clear();
+ if(b) {
+ checkArrays(a,b);
+ for(size_t i=0; i < size; i++)
+ if(read<bool>(a,i)) r->push((*b)[i]);
+ } else if(c) {
+ checkArrays(a,c);
+ for(size_t i=0; i < size; i++)
+ if(!read<bool>(a,i)) r->push((*c)[i]);
+ }
+ }
+ return r;
+}
+
+// Return an n x n identity matrix.
+realarray2 *identity(Int n)
+{
+ return Identity(n);
+}
+
+// Return the diagonal matrix with diagonal entries given by a.
+realarray2* :diagonal(realarray *a)
+{
+ size_t n=checkArray(a);
+ array *c=new array(n);
+ for(size_t i=0; i < n; ++i) {
+ array *ci=new array(n);
+ (*c)[i]=ci;
+ for(size_t j=0; j < i; ++j)
+ (*ci)[j]=0.0;
+ (*ci)[i]=read<real>(a,i);
+ for(size_t j=i+1; j < n; ++j)
+ (*ci)[j]=0.0;
+ }
+ return c;
+}
+
+// Return the inverse of an n x n matrix a using Gauss-Jordan elimination.
+realarray2 *inverse(realarray2 *a)
+{
+ a=copyArray2(a);
+ size_t n=checkArray(a);
+ checkSquare(a);
+
+ inverseAllocate(n);
+
+ for(size_t i=0; i < n; i++)
+ pivot[i]=0;
+
+ size_t col=0, row=0;
+ // This is the main loop over the columns to be reduced.
+ for(size_t i=0; i < n; i++) {
+ real big=0.0;
+ // This is the outer loop of the search for a pivot element.
+ for(size_t j=0; j < n; j++) {
+ array *aj=read<array*>(a,j);
+ if(pivot[j] != 1) {
+ for(size_t k=0; k < n; k++) {
+ if(pivot[k] == 0) {
+ real temp=fabs(read<real>(aj,k));
+ if(temp >= big) {
+ big=temp;
+ row=j;
+ col=k;
+ }
+ } else if(pivot[k] > 1) {
+ inverseDeallocate();
+ error(singular);
+ }
+ }
+ }
+ }
+ ++(pivot[col]);
+
+ // Interchange rows, if needed, to put the pivot element on the diagonal.
+ array *acol=read<array*>(a,col);
+ if(row != col) {
+ array *arow=read<array*>(a,row);
+ for(size_t l=0; l < n; l++) {
+ real temp=read<real>(arow,l);
+ (*arow)[l]=read<real>(acol,l);
+ (*acol)[l]=temp;
+ }
+ }
+
+ Row[i]=row;
+ Col[i]=col;
+
+ // Divide the pivot row by the pivot element.
+ real denom=read<real>(acol,col);
+ if(denom == 0.0) {
+ inverseDeallocate();
+ error(singular);
+ }
+ real pivinv=1.0/denom;
+ (*acol)[col]=1.0;
+ for(size_t l=0; l < n; l++)
+ (*acol)[l]=read<real>(acol,l)*pivinv;
+
+ // Reduce all rows except for the pivoted one.
+ for(size_t k=0; k < n; k++) {
+ if(k != col) {
+ array *ak=read<array*>(a,k);
+ real akcol=read<real>(ak,col);
+ (*ak)[col]=0.0;
+ for(size_t l=0; l < n; l++)
+ (*ak)[l]=read<real>(ak,l)-read<real>(acol,l)*akcol;
+ }
+ }
+ }
+
+ // Unscramble the inverse matrix in view of the column interchanges.
+ for(size_t l=n; l > 0;) {
+ l--;
+ size_t r=Row[l];
+ size_t c=Col[l];
+ if(r != c) {
+ for(size_t k=0; k < n; k++) {
+ array *ak=read<array*>(a,k);
+ real temp=read<real>(ak,r);
+ (*ak)[r]=read<real>(ak,c);
+ (*ak)[c]=temp;
+ }
+ }
+ }
+ inverseDeallocate();
+ return a;
+}
+
+// Solve the linear equation ax=b by LU decomposition, returning the
+// solution x, where a is an n x n matrix and b is an array of length n.
+// If no solution exists, return an empty array.
+realarray *solve(realarray2 *a, realarray *b, bool warn=true)
+{
+ size_t n=checkArray(a);
+
+ if(n == 0) return new array(0);
+
+ size_t m=checkArray(b);
+ if(m != n) error(incommensurate);
+
+ real *A=copyArray2C(a);
+ size_t *index=new size_t[n];
+
+ if(LUdecompose(A,n,index,warn) == 0)
+ return new array(0);
+
+ array *x=new array(n);
+
+ real *B=copyArrayC(b);
+
+ for(size_t i=0; i < n; ++i) {
+ size_t ip=index[i];
+ real sum=B[ip];
+ B[ip]=B[i];
+ real *Ai=A+i*n;
+ for(size_t j=0; j < i; ++j)
+ sum -= Ai[j]*B[j];
+ B[i]=sum;
+ }
+
+ for(size_t i=n; i > 0;) {
+ --i;
+ real sum=B[i];
+ real *Ai=A+i*n;
+ for(size_t j=i+1; j < n; ++j)
+ sum -= Ai[j]*B[j];
+ B[i]=sum/Ai[i];
+ }
+
+ for(size_t i=0; i < n; ++i)
+ (*x)[i]=B[i];
+
+ delete[] index;
+ delete[] B;
+ delete[] A;
+
+ return x;
+}
+
+// Solve the linear equation ax=b by LU decomposition, returning the
+// solution x, where a is an n x n matrix and b is an n x m matrix.
+// If no solution exists, return an empty array.
+realarray2 *solve(realarray2 *a, realarray2 *b, bool warn=true)
+{
+ size_t n=checkArray(a);
+
+ if(n == 0) return new array(0);
+
+ if(checkArray(b) != n) error(incommensurate);
+ size_t m=checkArray(read<array*>(b,0));
+
+ real *A=copyArray2C(a);
+ real *B=copyArray2C(b,false);
+
+ size_t *index=new size_t[n];
+
+ if(LUdecompose(A,n,index,warn) == 0)
+ return new array(0);
+
+ array *x=new array(n);
+
+ for(size_t i=0; i < n; ++i) {
+ real *Ai=A+i*n;
+ real *Bi=B+i*m;
+ real *Bip=B+index[i]*m;
+ for(size_t k=0; k < m; ++k) {
+ real sum=Bip[k];
+ Bip[k]=Bi[k];
+ size_t jk=k;
+ for(size_t j=0; j < i; ++j, jk += m)
+ sum -= Ai[j]*B[jk];
+ Bi[k]=sum;
+ }
+ }
+
+ for(size_t i=n; i > 0;) {
+ --i;
+ real *Ai=A+i*n;
+ real *Bi=B+i*m;
+ for(size_t k=0; k < m; ++k) {
+ real sum=Bi[k];
+ size_t jk=(i+1)*m+k;
+ for(size_t j=i+1; j < n; ++j, jk += m)
+ sum -= Ai[j]*B[jk];
+ Bi[k]=sum/Ai[i];
+ }
+ }
+
+ for(size_t i=0; i < n; ++i) {
+ real *Bi=B+i*m;
+ array *xi=new array(m);
+ (*x)[i]=xi;
+ for(size_t j=0; j < m; ++j)
+ (*xi)[j]=Bi[j];
+ }
+
+ delete[] index;
+ delete[] B;
+ delete[] A;
+
+ return x;
+}
+
+// Compute the determinant of an n x n matrix.
+real determinant(realarray2 *a)
+{
+ real *A=copyArray2C(a);
+ size_t n=checkArray(a);
+
+ real det=LUdecompose(A,n,NULL,false);
+ size_t n1=n+1;
+ for(size_t i=0; i < n; ++i)
+ det *= A[i*n1];
+
+ delete[] A;
+
+ return det;
+}
+
+realarray *Operator *(realarray2 *a, realarray *b)
+{
+ size_t n=checkArray(a);
+ size_t m=checkArray(b);
+ array *c=new array(n);
+ real *B=copyArrayC(b);
+ for(size_t i=0; i < n; ++i) {
+ array *ai=read<array*>(a,i);
+ if(checkArray(ai) != m) error(incommensurate);
+ real sum=0.0;
+ for(size_t j=0; j < m; ++j)
+ sum += read<real>(ai,j)*B[j];
+ (*c)[i]=sum;
+ }
+ delete[] B;
+ return c;
+}
+
+realarray2 *Operator *(realarray2 *a, realarray2 *b)
+{
+ size_t n=checkArray(a);
+
+ size_t nb=checkArray(b);
+ size_t na0=n == 0 ? 0 : checkArray(read<array*>(a,0));
+ if(na0 != nb)
+ error(incommensurate);
+
+ size_t nb0=nb == 0 ? 0 : checkArray(read<array*>(b,0));
+
+ array *c=new array(n);
+
+ real *A=copyArray2C(a,false);
+ real *B=copyArray2C(b,false);
+
+ for(size_t i=0; i < n; ++i) {
+ real *Ai=A+i*nb;
+ array *ci=new array(nb0);
+ (*c)[i]=ci;
+ for(size_t j=0; j < nb0; ++j) {
+ real sum=0.0;
+ size_t kj=j;
+ for(size_t k=0; k < nb; ++k, kj += nb0)
+ sum += Ai[k]*B[kj];
+ (*ci)[j]=sum;
+ }
+ }
+
+ delete[] B;
+ delete[] A;
+
+ return c;
+}
+
+triple Operator *(realarray2 *t, triple v)
+{
+ return *t*v;
+}
+
+pair project(triple v, realarray2 *t)
+{
+ size_t n=checkArray(t);
+ if(n != 4) error(incommensurate);
+ array *t0=read<array*>(t,0);
+ array *t1=read<array*>(t,1);
+ array *t3=read<array*>(t,3);
+ if(checkArray(t0) != 4 || checkArray(t1) != 4 || checkArray(t3) != 4)
+ error(incommensurate);
+
+ real x=v.getx();
+ real y=v.gety();
+ real z=v.getz();
+
+ real f=read<real>(t3,0)*x+read<real>(t3,1)*y+read<real>(t3,2)*z+
+ read<real>(t3,3);
+ if(f == 0.0) dividebyzero();
+ f=1.0/f;
+
+ return pair((read<real>(t0,0)*x+read<real>(t0,1)*y+read<real>(t0,2)*z+
+ read<real>(t0,3))*f,
+ (read<real>(t1,0)*x+read<real>(t1,1)*y+read<real>(t1,2)*z+
+ read<real>(t1,3))*f);
+}
+
+// Compute the dot product of vectors a and b.
+real dot(realarray *a, realarray *b)
+{
+ size_t n=checkArrays(a,b);
+ real sum=0.0;
+ for(size_t i=0; i < n; ++i)
+ sum += read<real>(a,i)*read<real>(b,i);
+ return sum;
+}
+
+// Solve the problem L\inv f, where f is an n vector and L is the n x n matrix
+//
+// [ b[0] c[0] a[0] ]
+// [ a[1] b[1] c[1] ]
+// [ a[2] b[2] c[2] ]
+// [ ... ]
+// [ c[n-1] a[n-1] b[n-1] ]
+realarray *tridiagonal(realarray *a, realarray *b, realarray *c, realarray *f)
+{
+ size_t n=checkArrays(a,b);
+ checkEqual(n,checkArray(c));
+ checkEqual(n,checkArray(f));
+
+ array *up=new array(n);
+ array& u=*up;
+
+ if(n == 0) return up;
+
+ // Special case: zero Dirichlet boundary conditions
+ if(read<real>(a,0) == 0.0 && read<real>(c,n-1) == 0.0) {
+ real temp=read<real>(b,0);
+ if(temp == 0.0) dividebyzero();
+ temp=1.0/temp;
+
+ real *work=new real[n];
+ u[0]=read<real>(f,0)*temp;
+ work[0]=-read<real>(c,0)*temp;
+
+ for(size_t i=1; i < n; i++) {
+ real temp=(read<real>(b,i)+read<real>(a,i)*work[i-1]);
+ if(temp == 0.0) {delete[] work; dividebyzero();}
+ temp=1.0/temp;
+ u[i]=(read<real>(f,i)-read<real>(a,i)*read<real>(u,i-1))*temp;
+ work[i]=-read<real>(c,i)*temp;
+ }
+
+ for(size_t i=n-1; i >= 1; i--)
+ u[i-1]=read<real>(u,i-1)+work[i-1]*read<real>(u,i);
+
+ delete[] work;
+ return up;
+ }
+
+ real binv=read<real>(b,0);
+ if(binv == 0.0) dividebyzero();
+ binv=1.0/binv;
+
+ if(n == 1) {u[0]=read<real>(f,0)*binv; return up;}
+ if(n == 2) {
+ real factor=(read<real>(b,0)*read<real>(b,1)-
+ read<real>(a,0)*read<real>(c,1));
+ if(factor== 0.0) dividebyzero();
+ factor=1.0/factor;
+ real temp=(read<real>(b,0)*read<real>(f,1)-
+ read<real>(c,1)*read<real>(f,0))*factor;
+ u[0]=(read<real>(b,1)*read<real>(f,0)-
+ read<real>(a,0)*read<real>(f,1))*factor;
+ u[1]=temp;
+ return up;
+ }
+
+ real *gamma=new real[n-2];
+ real *delta=new real[n-2];
+
+ gamma[0]=read<real>(c,0)*binv;
+ delta[0]=read<real>(a,0)*binv;
+ u[0]=read<real>(f,0)*binv;
+ real beta=read<real>(c,n-1);
+ real fn=read<real>(f,n-1)-beta*read<real>(u,0);
+ real alpha=read<real>(b,n-1)-beta*delta[0];
+
+ for(size_t i=1; i <= n-3; i++) {
+ real alphainv=read<real>(b,i)-read<real>(a,i)*gamma[i-1];
+ if(alphainv == 0.0) {delete[] gamma; delete[] delta; dividebyzero();}
+ alphainv=1.0/alphainv;
+ beta *= -gamma[i-1];
+ gamma[i]=read<real>(c,i)*alphainv;
+ u[i]=(read<real>(f,i)-read<real>(a,i)*read<real>(u,i-1))*alphainv;
+ fn -= beta*read<real>(u,i);
+ delta[i]=-read<real>(a,i)*delta[i-1]*alphainv;
+ alpha -= beta*delta[i];
+ }
+
+ real alphainv=read<real>(b,n-2)-read<real>(a,n-2)*gamma[n-3];
+ if(alphainv == 0.0) {delete[] gamma; delete[] delta; dividebyzero();}
+ alphainv=1.0/alphainv;
+ u[n-2]=(read<real>(f,n-2)-read<real>(a,n-2)*read<real>(u,n-3))
+ *alphainv;
+ beta=read<real>(a,n-1)-beta*gamma[n-3];
+ real dnm1=(read<real>(c,n-2)-read<real>(a,n-2)*delta[n-3])*alphainv;
+ real temp=alpha-beta*dnm1;
+ if(temp == 0.0) {delete[] gamma; delete[] delta; dividebyzero();}
+ u[n-1]=temp=(fn-beta*read<real>(u,n-2))/temp;
+ u[n-2]=read<real>(u,n-2)-dnm1*temp;
+
+ for(size_t i=n-2; i >= 1; i--)
+ u[i-1]=read<real>(u,i-1)-gamma[i-1]*read<real>(u,i)-delta[i-1]*temp;
+
+ delete[] delta;
+ delete[] gamma;
+
+ return up;
+}
+
+// Root solve by Newton-Raphson
+real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x,
+ bool verbose=false)
+{
+ static const real fuzz=1000.0*DBL_EPSILON;
+ Int i=0;
+ size_t oldPrec=0;
+ if(verbose)
+ oldPrec=cout.precision(DBL_DIG);
+
+ real diff=DBL_MAX;
+ real lastdiff;
+ do {
+ real x0=x;
+
+ Stack->push(x);
+ fprime->call(Stack);
+ real dfdx=pop<real>(Stack);
+
+ if(dfdx == 0.0) {
+ x=DBL_MAX;
+ break;
+ }
+
+ Stack->push(x);
+ f->call(Stack);
+ real fx=pop<real>(Stack);
+
+ x -= fx/dfdx;
+
+ lastdiff=diff;
+
+ if(verbose)
+ cout << "Newton-Raphson: " << x << endl;
+
+ diff=fabs(x-x0);
+ if(++i == iterations) {
+ x=DBL_MAX;
+ break;
+ }
+ } while (diff != 0.0 && (diff < lastdiff || diff > fuzz*fabs(x)));
+
+ if(verbose)
+ cout.precision(oldPrec);
+ return x;
+}
+
+// Root solve by Newton-Raphson bisection
+// cf. routine rtsafe (Press et al., Numerical Recipes, 1991).
+real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x1,
+ real x2, bool verbose=false)
+{
+ static const real fuzz=1000.0*DBL_EPSILON;
+ size_t oldPrec=0;
+ if(verbose)
+ oldPrec=cout.precision(DBL_DIG);
+
+ Stack->push(x1);
+ f->call(Stack);
+ real f1=pop<real>(Stack);
+ if(f1 == 0.0) return x1;
+
+ Stack->push(x2);
+ f->call(Stack);
+ real f2=pop<real>(Stack);
+ if(f2 == 0.0) return x2;
+
+ if((f1 > 0.0 && f2 > 0.0) || (f1 < 0.0 && f2 < 0.0)) {
+ ostringstream buf;
+ buf << "root not bracketed, f(x1)=" << f1 << ", f(x2)=" << f2 << endl;
+ error(buf);
+ }
+
+ real x=0.5*(x1+x2);
+ real dxold=fabs(x2-x1);
+ if(f1 > 0.0) {
+ real temp=x1;
+ x1=x2;
+ x2=temp;
+ }
+
+ if(verbose)
+ cout << "midpoint: " << x << endl;
+
+ real dx=dxold;
+ Stack->push(x);
+ f->call(Stack);
+ real y=pop<real>(Stack);
+
+ Stack->push(x);
+ fprime->call(Stack);
+ real dy=pop<real>(Stack);
+
+ Int j;
+ for(j=0; j < iterations; j++) {
+ if(((x-x2)*dy-y)*((x-x1)*dy-y) >= 0.0 || fabs(2.0*y) > fabs(dxold*dy)) {
+ dxold=dx;
+ dx=0.5*(x2-x1);
+ x=x1+dx;
+ if(verbose)
+ cout << "bisection: " << x << endl;
+ if(x1 == x) return x;
+ } else {
+ dxold=dx;
+ dx=y/dy;
+ real temp=x;
+ x -= dx;
+ if(verbose)
+ cout << "Newton-Raphson: " << x << endl;
+ if(temp == x) return x;
+ }
+ if(fabs(dx) < fuzz*fabs(x)) return x;
+
+ Stack->push(x);
+ f->call(Stack);
+ y=pop<real>(Stack);
+
+ Stack->push(x);
+ fprime->call(Stack);
+ dy=pop<real>(Stack);
+
+ if(y < 0.0) x1=x;
+ else x2=x;
+ }
+ if(verbose)
+ cout.precision(oldPrec);
+ return (j == iterations) ? DBL_MAX : x;
+}
+
+real simpson(callableReal *f, real a, real b, real acc=DBL_EPSILON,
+ real dxmax=0)
+{
+ real integral;
+ if(dxmax == 0) dxmax=b-a;
+ Func=f;
+ FuncStack=Stack;
+ if(!simpson(integral,wrapFunction,a,b,acc,dxmax))
+ error("nesting capacity exceeded in simpson");
+ return integral;
+}
+
+// Compute the fast Fourier transform of a pair array
+pairarray* :pairArrayFFT(pairarray *a, Int sign=1)
+{
+ unsigned n=(unsigned) checkArray(a);
+#ifdef HAVE_LIBFFTW3
+ array *c=new array(n);
+ if(n) {
+ Complex *f=FFTWComplex(n);
+ fft1d Forward(n,intcast(sign),f);
+
+ for(size_t i=0; i < n; i++) {
+ pair z=read<pair>(a,i);
+ f[i]=Complex(z.getx(),z.gety());
+ }
+ Forward.fft(f);
+
+ for(size_t i=0; i < n; i++) {
+ Complex z=f[i];
+ (*c)[i]=pair(z.real(),z.imag());
+ }
+ FFTWdelete(f);
+ }
+#else
+ unused(&n);
+ unused(&sign);
+ array *c=new array(0);
+#endif // HAVE_LIBFFTW3
+ return c;
+}
+
+Intarray2 *triangulate(pairarray *z)
+{
+ size_t nv=checkArray(z);
+// Call robust version of Gilles Dumoulin's port of Paul Bourke's
+// triangulation code.
+
+ XYZ *pxyz=new XYZ[nv+3];
+ ITRIANGLE *V=new ITRIANGLE[4*nv];
+
+ for(size_t i=0; i < nv; ++i) {
+ pair w=read<pair>(z,i);
+ pxyz[i].p[0]=w.getx();
+ pxyz[i].p[1]=w.gety();
+ pxyz[i].i=(Int) i;
+ }
+
+ Int ntri;
+ Triangulate((Int) nv,pxyz,V,ntri,true,false);
+
+ size_t nt=(size_t) ntri;
+ array *t=new array(nt);
+ for(size_t i=0; i < nt; ++i) {
+ array *ti=new array(3);
+ (*t)[i]=ti;
+ ITRIANGLE *Vi=V+i;
+ (*ti)[0]=pxyz[Vi->p1].i;
+ (*ti)[1]=pxyz[Vi->p2].i;
+ (*ti)[2]=pxyz[Vi->p3].i;
+ }
+
+ delete[] V;
+ delete[] pxyz;
+ return t;
+}
+
+
+// File operations
+
+bool ==(file *a, file *b)
+{
+ return a == b;
+}
+
+bool !=(file *a, file *b)
+{
+ return a != b;
+}
+
+file* :nullFile()
+{
+ return &camp::nullfile;
+}
+
+file* input(string name, bool check=true, string comment=commentchar)
+{
+ char c=comment.empty() ? (char) 0 : comment[0];
+ file *f=new ifile(name,c,check);
+ f->open();
+ return f;
+}
+
+file* output(string name, bool update=false, string comment=commentchar)
+{
+ file *f;
+ if(update) {
+ char c=comment.empty() ? (char) 0 : comment[0];
+ f=new iofile(name,c);
+ } else f=new ofile(name);
+ f->open();
+ if(update) f->seek(0,false);
+ return f;
+}
+
+file* xinput(string name, bool check=true)
+{
+#ifdef HAVE_RPC_RPC_H
+ file *f=new ixfile(name,check);
+ f->open();
+ return f;
+#else
+ ostringstream buf;
+ buf << name << ": XDR read support not enabled";
+ error(buf);
+ unused(&check); // Suppress unused variable warning
+#endif
+}
+
+file* xoutput(string name, bool update=false)
+{
+#ifdef HAVE_RPC_RPC_H
+ file *f;
+ if(update)
+ f=new ioxfile(name);
+ else f=new oxfile(name);
+ f->open();
+ if(update) f->seek(0,false);
+ return f;
+#else
+ ostringstream buf;
+ buf << name << ": XDR write support not enabled";
+ error(buf);
+ unused(&update); // Suppress unused variable warning
+#endif
+}
+
+file* binput(string name, bool check=true)
+{
+ file *f=new ibfile(name,check);
+ f->open();
+ return f;
+}
+
+file* boutput(string name, bool update=false)
+{
+ file *f;
+ if(update) f=new iobfile(name);
+ else f=new obfile(name);
+ f->open();
+ if(update) f->seek(0,false);
+ return f;
+}
+
+bool eof(file *File)
+{
+ return File->eof();
+}
+
+bool eol(file *File)
+{
+ return File->eol();
+}
+
+bool error(file *File)
+{
+ return File->error();
+}
+
+void clear(file *File)
+{
+ File->clear();
+}
+
+void close(file *File)
+{
+ File->close();
+}
+
+Int precision(file *File=NULL, Int digits=0)
+{
+ if(File == 0) File=&camp::Stdout;
+ return File->precision(digits);
+}
+
+void flush(file *File)
+{
+ File->flush();
+}
+
+string getc(file *File)
+{
+ char c=0;
+ if(File->isOpen()) File->read(c);
+ static char str[1];
+ str[0]=c;
+ return string(str);
+}
+
+Int tell(file *File)
+{
+ return File->tell();
+}
+
+void seek(file *File, Int pos)
+{
+ File->seek(pos,pos >= 0);
+}
+
+void seekeof(file *File)
+{
+ File->seek(0,false);
+}
+
+// Set file dimensions
+file* dimension(file *File, Int nx)
+{
+ File->dimension(nx);
+ return File;
+}
+
+file* dimension(file *File, Int nx, Int ny)
+{
+ File->dimension(nx,ny);
+ return File;
+}
+
+file* dimension(file *File, Int nx, Int ny, Int nz)
+{
+ File->dimension(nx,ny,nz);
+ return File;
+}
+
+// Set file to read comma-separated values
+file* csv(file *File, bool b=true)
+{
+ File->CSVMode(b);
+ return File;
+}
+
+// Set file to read whitespace-separated values
+file* word(file *File, bool b=true)
+{
+ File->WordMode(b);
+ return File;
+}
+
+// Set file to read arrays in line-at-a-time mode
+file* line(file *File, bool b=true)
+{
+ File->LineMode(b);
+ return File;
+}
+
+// Set file to read/write single-precision XDR values.
+file* single(file *File, bool b=true)
+{
+ File->SingleReal(b);
+ File->SingleInt(b);
+ return File;
+}
+
+// Set file to read/write single-precision real XDR values.
+file* single(file *File, real x, bool b=true)
+{
+ File->SingleReal(b);
+ unused(&x);
+ return File;
+}
+
+// Set file to read/write single-precision int XDR values.
+file* single(file *File, Int x, bool b=true)
+{
+ File->SingleInt(b);
+ unused(&x);
+ return File;
+}
+
+// Set file to read an array1 (1 Int size followed by a 1d array)
+file* read1(file *File)
+{
+ File->dimension(-2);
+ return File;
+}
+
+// Set file to read an array2 (2 Int sizes followed by a 2d array)
+file* read2(file *File)
+{
+ File->dimension(-2,-2);
+ return File;
+}
+
+// Set file to read an array3 (3 Int sizes followed by a 3d array)
+file* read3(file *File)
+{
+ File->dimension(-2,-2,-2);
+ return File;
+}
+
+// Return the last n lines of the history named name.
+stringarray* history(string name, Int n=1)
+{
+#if defined(HAVE_LIBREADLINE) && defined(HAVE_LIBCURSES)
+ bool newhistory=historyMap.find(name) == historyMap.end();
+
+ string filename;
+
+ if(newhistory) {
+ filename=historyfilename(name);
+ std::ifstream exists(filename.c_str());
+ if(!exists) return new array(0);
+ }
+
+ store_history(&history_save);
+ HISTORY_STATE& history=historyMap[name].state;
+ history_set_history_state(&history);
+
+ if(newhistory)
+ read_history(filename.c_str());
+
+ array *a=get_history(n);
+
+ store_history(&history);
+ history_set_history_state(&history_save);
+
+ return a;
+#else
+ unused(&n);
+ return new array(0);
+#endif
+}
+
+// Return the last n lines of the interactive history.
+stringarray* history(Int n=0)
+{
+#if defined(HAVE_LIBREADLINE) && defined(HAVE_LIBCURSES)
+ return get_history(n);
+#else
+ unused(&n);
+ return new array(0);
+#endif
+}
+
+// Prompt for a string using prompt, the GNU readline library, and a
+// local history named name.
+string readline(string prompt=emptystring, string name=emptystring,
+ bool tabcompletion=false)
+{
+ if(!isatty(STDIN_FILENO))
+ return emptystring;
+#if defined(HAVE_LIBREADLINE) && defined(HAVE_LIBCURSES)
+ init_readline(tabcompletion);
+
+ store_history(&history_save);
+ bool newhistory=historyMap.find(name) == historyMap.end();
+ historyState& h=historyMap[name];
+ HISTORY_STATE& history=h.state;
+ history_set_history_state(&history);
+
+ if(newhistory)
+ read_history(historyfilename(name).c_str());
+
+ static char *line=NULL;
+ /* Return the memory to the free pool
+ if the buffer has already been allocated. */
+ if(line) {
+ free(line);
+ line=NULL;
+ }
+
+ /* Get a line from the user. */
+ line=readline(prompt.c_str());
+
+ if(!line) cout << endl;
+
+ history_set_history_state(&history_save);
+
+ return line ? string(line) : emptystring;
+#else
+ cout << prompt;
+ string s;
+ getline(cin,s);
+ unused(&tabcompletion); // Avoid unused variable warning message.
+ return s;
+#endif
+}
+
+// Save a string in a local history named name.
+// If store=true, store the local history in the file historyfilename(name).
+void saveline(string name, string value, bool store=true)
+{
+#if defined(HAVE_LIBREADLINE) && defined(HAVE_LIBCURSES)
+ store_history(&history_save);
+ bool newhistory=historyMap.find(name) == historyMap.end();
+ historyState& h=historyMap[name];
+ h.store=store;
+ HISTORY_STATE& history=h.state;
+ history_set_history_state(&history);
+
+ if(newhistory)
+ read_history(historyfilename(name).c_str());
+
+ if(value != "") {
+ add_history(value.c_str());
+ if(store) {
+ std::ofstream hout(historyfilename(name).c_str(),std::ios::app);
+ hout << value << endl;
+ }
+ }
+
+ store_history(&history);
+ history_set_history_state(&history_save);
+#else
+ unused(&store);
+#endif
+}
+
+void generate_random_backtrace()
+{
+#if defined(USEGC) && defined(GC_DEBUG) && defined(GC_BACKTRACE)
+ GC_generate_random_backtrace();
+#else
+ error("generate_random_backtrace() requires ./configure --enable-gc-debug");
+#endif
+}
+
+void print_random_addresses(Int n=1)
+{
+#if defined(USEGC) && defined(GC_DEBUG) && defined(GC_BACKTRACE)
+ GC_gcollect();
+ for (Int i=0; i < n; ++i)
+ GC_debug_print_heap_obj_proc(GC_base(GC_generate_random_valid_address()));
+#else
+ error("print_random_addresses() requires ./configure --enable-gc-debug");
+ unused(&n); // Avoid unused variable warning message.
+#endif
+}