diff options
Diffstat (limited to 'Build/source/utils/asymptote/runtime.in')
-rw-r--r-- | Build/source/utils/asymptote/runtime.in | 5818 |
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 +} |