/***** * 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 #include #include #include #include #include #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 #include #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 #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(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(a,i); return c; } array *copyArray2(array *a) { size_t size=checkArray(a); array *c=new array(size); for(size_t i=0; i < size; i++) { array *ai=read(a,i); size_t aisize=checkArray(ai); array *ci=new array(aisize); (*c)[i]=ci; for(size_t j=0; j < aisize; j++) (*ci)[j]=(*ai)[j]; } return c; } 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(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(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(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(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(ai,j); } else error(square ? "matrix must be square" : "matrix must be rectangular"); } return c; } triple *copyTripleArray2C(array *a, bool square=true, size_t dim2=0) { size_t n=checkArray(a); size_t m=(square || n == 0) ? n : checkArray(read(a,0)); if(n > 0 && dim2 && m != dim2) { ostringstream buf; buf << "second matrix dimension must be " << dim2 << endl; error(buf); } triple *c=new triple[n*m]; for(size_t i=0; i < n; i++) { array *ai=read(a,i); size_t aisize=checkArray(ai); if(aisize == m) { triple *ci=c+i*m; for(size_t j=0; j < m; j++) ci[j]=read(ai,j); } else error(square ? "matrix must be square" : "matrix must be rectangular"); } return c; } double *copyTripleArray2Components(array *a, bool square=true, size_t dim2=0) { size_t n=checkArray(a); size_t m=(square || n == 0) ? n : checkArray(read(a,0)); if(n > 0 && dim2 && m != dim2) { ostringstream buf; buf << "second matrix dimension must be " << dim2 << endl; error(buf); } size_t nm=n*m; double *cx=new double[3*nm]; double *cy=cx+nm; double *cz=cx+2*nm; for(size_t i=0; i < n; i++) { array *ai=read(a,i); size_t aisize=checkArray(ai); if(aisize == m) { double *xi=cx+i*m; double *yi=cy+i*m; double *zi=cz+i*m; for(size_t j=0; j < m; j++) { triple v=read(ai,j); xi[j]=v.getx(); yi[j]=v.gety(); zi[j]=v.getz(); } } else error(square ? "matrix must be square" : "matrix must be rectangular"); } return cx; } double norm(double *a, size_t n) { if(n == 0) return 0.0; double M=fabs(a[0]); for(size_t i=1; i < n; ++i) M=::max(M,fabs(a[i])); return M; } double norm(triple *a, size_t n) { if(n == 0) return 0.0; double M=a[0].abs2(); for(size_t i=1; i < n; ++i) M=::max(M,a[i].abs2()); return sqrt(M); } static 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(t,0); array *t1=read(t,1); array *t2=read(t,2); array *t3=read(t,3); if(checkArray(t0) != 4 || checkArray(t1) != 4 || checkArray(t2) != 4 || checkArray(t3) != 4) error(incommensurate); double x=v.getx(); double y=v.gety(); double z=v.getz(); double f=read(t3,0)*x+read(t3,1)*y+read(t3,2)*z+ read(t3,3); if(f == 0.0) run::dividebyzero(); f=1.0/f; return triple((read(t0,0)*x+read(t0,1)*y+read(t0,2)*z+ read(t0,3))*f, (read(t1,0)*x+read(t1,1)*y+read(t1,2)*z+ read(t1,3))*f, (read(t2,0)*x+read(t2,1)*y+read(t2,2)*z+ read(t2,3))*f); } triple multshiftless(const array& t, const triple& v) { size_t n=checkArray(&t); if(n != 4) error(incommensurate); array *t0=read(t,0); array *t1=read(t,1); array *t2=read(t,2); array *t3=read(t,3); if(checkArray(t0) != 4 || checkArray(t1) != 4 || checkArray(t2) != 4 || checkArray(t3) != 4) error(incommensurate); double x=v.getx(); double y=v.gety(); double z=v.getz(); double f=read(t3,0)*x+read(t3,1)*y+read(t3,2)*z+ read(t3,3); if(f == 0.0) run::dividebyzero(); f=1.0/f; return triple((read(t0,0)*x+read(t0,1)*y+read(t0,2)*z)*f, (read(t1,0)*x+read(t1,1)*y+read(t1,2)*z)*f, (read(t2,0)*x+read(t2,1)*y+read(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(s,NULL); string S=pop(s); vm::item it=pop(s); bool defaultfile=isdefault(it); camp::file *f=defaultfile ? &camp::Stdout : vm::get(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(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 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("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(curPos.filename()); Stack->push((Int) curPos.Line()); Stack->push((Int) curPos.Column()); Stack->push(r ? r : item(def)); atBreakpointFunction->call(Stack); // returns a string } else Stack->push(""); } } 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::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(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(FuncStack); } callable *compareFunc; bool compareFunction(const vm::item& i, const vm::item& j) { FuncStack->push(i); FuncStack->push(j); compareFunc->call(FuncStack); return pop(FuncStack); } void unused(void *) { } static const double twelvepercent=0.12; static const double tenpercent=0.1; pair readpair(stringstream& s, double hscale=1.0, double vscale=1.0) { double x,y; s >> y; s >> x; return pair(hscale*x,vscale*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 keep, double hscale=1.0, double vscale=1.0) { array *P=new array(0); ostringstream pipe; pipe << "'" << getSetting("gs") << "' -q -dNOPAUSE -dBATCH"; if(safe) pipe << " -dSAFER"; #ifdef __CYGWIN__ const string null="NUL"; #else const 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; if(out.empty() && !gs.running()) break; buf << out; } if(verbose > 2) cout << endl; mem::vector 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,hscale,vscale)-offset; node.straight=false; break; } case 'L': { pair point=readpair(buf,hscale,vscale)-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,hscale,vscale)-offset; pair pre=readpair(buf,hscale,vscale)-offset; node.post=readpair(buf,hscale,vscale)-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,hscale,vscale); 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; } void cannotread(const string& s) { ostringstream buf; buf << "Cannot read from " << s << endl; error(buf); } void cannotwrite(const string& s) { ostringstream buf; buf << "Cannot write to " << s << endl; error(buf); } // 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; bool pdf=settings::pdf(getSetting("tex")); texfile tex(texname,b,true); tex.miniprologue(); tex.setfont(p); if(!pdf) { 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.verbatimline(*s+"%"); tex.epilogue(true); tex.close(); // TODO: Put in common subproutine (cf. teprocess). string aux=auxname(prefix,"aux"); unlink(aux.c_str()); string program=texprogram(); ostringstream cmd; bool context=settings::context(getSetting("tex")); cmd << program << (context ? " --nonstopmode '" : " \\nonstopmode\\input '") << texname << "'"; bool quiet=verbose <= 2; int status=System(cmd,quiet ? 1 : 0,true,"texpath",texpathmessage()); if(!status && getSetting("twice")) status=System(cmd,quiet ? 1 : 0,true,"texpath",texpathmessage()); if(status) { if(quiet) { ostringstream cmd; cmd << program << (context ? " --scrollmode '" : " \\scrollmode\\input '") << texname << "'"; System(cmd,0); } } string pdfname; if(!status) { if(pdf) { pdfname=auxname(prefix,"pdf"); std::ofstream ps(psname.c_str()); if(!ps) cannotwrite(psname); ps << ASYo << newl << "/ASY1 true def" << newl << "/fill {ASY1 {(Z) " << currentpoint << "/ASY1 false def} if " << pathforall << " (M) " << currentpoint << "currentpoint newpath moveto } bind def" << newl << "/stroke {ASY1 {(Z) " << currentpoint << "/ASY1 false def} if strokepath " << pathforall << " (M) " << currentpoint << "currentpoint newpath moveto } bind def" << newl; ostringstream cmd; cmd << "'" << getSetting("gs") << "' -q -dNOCACHE -dNOPAUSE -dBATCH"; if(safe) cmd << " -dSAFER"; cmd << " -sDEVICE=epswrite -sOutputFile=- '" << pdfname << "'"; iopipestream gs(cmd.str().c_str(),"gs","Ghostscript"); gs.block(false); while(true) { string line; gs >> line; if(line.empty() && !gs.running()) break; ps << line; } ps.close(); } else { ostringstream cmd; cmd << "'" << getSetting("dvips") << "' -R -Pdownload35 -D600 " << getSetting("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("keep"); if(!keep) { // Delete temporary files. unlink(texname.c_str()); if(!getSetting("keepaux")) unlink(aux.c_str()); unlink(auxname(prefix,"log").c_str()); if(pdf) unlink(pdfname.c_str()); else unlink(dviname.c_str()); if(context) { unlink(auxname(prefix,"top").c_str()); unlink(auxname(prefix,"tua").c_str()); unlink(auxname(prefix,"tui").c_str()); } } return readpath(psname,keep,twelvepercent, pdf ? twelvepercent : -twelvepercent); } patharray *textpath(string *s, pen p=CURRENTPEN) { array *P=new array(0); if(s->empty()) return P; string prefix=outname(); spaceToUnderscore(prefix); string outputname=auxname(prefix,getSetting("textoutputtype")); string textname=auxname(prefix,getSetting("textextension")); std::ofstream text(textname.c_str()); if(!text) cannotwrite(textname); text << getSetting("textprologue") << newl << p.Font() << newl << *s << newl << getSetting("textepilogue") << endl; text.close(); string psname=auxname(prefix,"ps"); std::ofstream ps(psname.c_str()); if(!ps) cannotwrite(psname); ps << ASYo << newl << "/ASY1 true def" << newl << "/fill {ASY1 {(Z) " << currentpoint << "/ASY1 false def} if " << pathforall << " (M) " << currentpoint << "currentpoint newpath moveto } bind def" << newl << "/stroke {ASY1 {(Z) " << currentpoint << "/ASY1 false def} if strokepath " << pathforall << " (M) " << currentpoint << "currentpoint newpath moveto } bind def" << endl; ostringstream pipe; pipe << getSetting("textcommand") << " '" << textname << "'"; iopipestream typesetter(pipe.str().c_str()); ostringstream pipe2; pipe2 << "'" << getSetting("gs") << "' -q -dNOCACHE -dNOPAUSE -dBATCH"; if(safe) pipe2 << " -dSAFER"; pipe2 << " -sDEVICE=epswrite -sOutputFile=- -"; iopipestream gs(pipe2.str().c_str(),"gs","Ghostscript"); gs.block(false); // TODO: Simplify by connecting the pipes directly. while(true) { string out; if(typesetter.isopen()) { typesetter >> out; if(!out.empty()) gs << out; else if(!typesetter.running()) { typesetter.pipeclose(); gs.eof(); } } string out2; gs >> out2; if(out2.empty() && !gs.running()) break; ps << out2; } ps.close(); if(verbose > 2) cout << endl; bool keep=getSetting("keep"); if(!keep) // Delete temporary files. unlink(textname.c_str()); return readpath(psname,keep,tenpercent,tenpercent); } 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); ps.epilogue(); ps.close(); return readpath(psname,getSetting("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(g,j); assert(P); Int stop=Min(P->size(),in-k); mem::vector& 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("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("tex")); } bool pdf() { return pdf(getSetting("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(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, realarray *background, 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), copyArrayC(background),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("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 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("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("convert") << "' -density " << res << "x" << res << " -transparent white"; hint=converter="convert"; } } picture *group=new picture; xform->call(Stack); transform t=pop(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("keep")) { for(mem::list::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(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(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 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; } 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 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()); 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 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 nodes(n); for(size_t i=0; i < n; ++i) { nodes[i].pre=read(pre,i); nodes[i].point=read(point,i); nodes[i].post=read(post,i); nodes[i].straight=read(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 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 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()); return V; } realarray2* intersections(path3 p, triplearray2 *P, real fuzz=-1) { triple *A=copyTripleArray2C(P,true,4); if(fuzz <= 0) fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())), norm(A,16)); std::vector T,U,V; intersections(T,U,V,p,A,fuzz); delete[] A; size_t n=T.size(); array *W=new array(n); for(size_t i=0; i < n; ++i) { array *Wi=new array(3); (*W)[i]=Wi; (*Wi)[0]=T[i]; (*Wi)[1]=U[i]; (*Wi)[2]=V[i]; } return W; // Sorting will done in asy. } 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(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(0)); if (n==1) { v.push_back(&curlout); v.push_back(&curlin); } else for (size_t i=1; iread(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) { if(g.size() > 0) 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); } triple minbound(triplearray2 *P, triple b) { real *A=copyTripleArray2Components(P,true,4); b=triple(bound(A,::min,b.getx(),sqrtFuzz*norm(A,16)), bound(A+16,::min,b.gety(),sqrtFuzz*norm(A+16,16)), bound(A+32,::min,b.getz(),sqrtFuzz*norm(A+32,16))); delete[] A; return b; } triple maxbound(triplearray2 *P, triple b) { real *A=copyTripleArray2Components(P,true,4); b=triple(bound(A,::max,b.getx(),sqrtFuzz*norm(A,16)), bound(A+16,::max,b.gety(),sqrtFuzz*norm(A+16,16)), bound(A+32,::max,b.getz(),sqrtFuzz*norm(A+32,16))); delete[] A; return b; } pair minbound(triplearray2 *P, realarray2 *t, pair b) { triple *A=copyTripleArray2C(P,true,4); real *T=copyArray2C(t,4); real fuzz=sqrtFuzz*norm(A,16); b=pair(bound(A,::min,xproject,T,b.getx(),fuzz), bound(A,::min,yproject,T,b.gety(),fuzz)); delete[] T; delete[] A; return b; } pair maxbound(triplearray2 *P, realarray2 *t, pair b) { triple *A=copyTripleArray2C(P,true,4); real *T=copyArray2C(t,4); real fuzz=sqrtFuzz*norm(A,16); b=pair(bound(A,::max,xproject,T,b.getx(),fuzz), bound(A,::max,yproject,T,b.gety(),fuzz)); 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(a,0); size_t m=checkArray(a0); if(m == 0) return 0.0; triple a00=vm::read(a0,0); real M=0.0; for(size_t i=0; i < n; ++i) { vm::array *ai=vm::read(a,i); size_t m=checkArray(ai); for(size_t j=0; j < m; ++j) { real a=(vm::read(ai,j)-a00).abs2(); if(a > M) M=a; } } 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(translate,i); checkArray(a); } const char *p=S->c_str(); ostringstream buf; while(*p) { for(size_t i=0; i < size;) { array *a=read(translate,i); string* from=read(a,0); size_t len=from->length(); if(strncmp(p,from->c_str(),len) != 0) {i++; continue;} buf << read(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, bool warn=true) { if(!warn && v.getx() == 0.0 && v.gety() == 0.0 && v.getz() == 0.0) return 0.0; return v.polar(); } real azimuth(triple v, bool warn=true) { if(!warn && v.getx() == 0.0 && v.gety() == 0.0) return 0.0; return v.azimuth(); } real colatitude(triple v, bool warn=true) { if(!warn && v.getx() == 0.0 && v.gety() == 0.0 && v.getz() == 0.0) return 0.0; return degrees(v.polar()); } real latitude(triple v, bool warn=true) { if(!warn && v.getx() == 0.0 && v.gety() == 0.0 && v.getz() == 0.0) return 0.0; 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::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("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("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(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; ipush((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(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(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((*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(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(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(i)); } array *result=new array(resultSize); size_t ri=0; for (size_t i=0; i < numArgs; ++i) { array *arg=a->read(i); size_t size=checkArray(arg); for (size_t j=0; j < size; ++j) { (*result)[ri]=(*arg)[j]; ++ri; } } 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(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(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(perm,i); size_t P=(size_t) p; if(p < 0 || P >= DIM) { ostringstream buf; buf << "permutation index out of range: " << p; error(buf); } size[P]=P; } for(size_t i=0; i < DIM; ++i) if(size[i] == DIM) error("permutation indices must be distinct"); static const char *rectangular= "3D transpose implemented for rectangular matrices only"; size_t isize=size[0]=checkArray(a); array *a0=read(a,0); size[1]=checkArray(a0); array *a00=read(a0,0); size[2]=checkArray(a00); for(size_t i=0; i < isize; i++) { array *ai=read(a,i); size_t jsize=checkArray(ai); if(jsize != size[1]) error(rectangular); for(size_t j=0; j < jsize; j++) { array *aij=read(ai,j); if(checkArray(aij) != size[2]) error(rectangular); } } size_t perm0=(size_t) read(perm,0); size_t perm1=(size_t) read(perm,1); size_t perm2=(size_t) read(perm,2); size_t sizep0=size[perm0]; size_t sizep1=size[perm1]; size_t sizep2=size[perm2]; array *c=new array(sizep0); for(size_t i=0; i < sizep0; ++i) { array *ci=new array(sizep1); (*c)[i]=ci; for(size_t j=0; j < sizep1; ++j) { array *cij=new array(sizep2); (*ci)[j]=cij; } } size_t* i=new size_t[DIM]; for(i[0]=0; i[0] < size[0]; ++i[0]) { array *a0=read(a,i[0]); for(i[1]=0; i[1] < size[1]; ++i[1]) { array *a1=read(a0,i[1]); for(i[2]=0; i[2] < size[2]; ++i[2]) { array *c0=read(c,i[perm0]); array *c1=read(c0,i[perm1]); (*c1)[i[perm2]]=read(a1,i[2]); } } } delete [] i; delete [] size; 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(a,i)) { n--; if(n == 0) {j=(Int) i; break;} } if(n < 0) for(size_t i=size; i > 0;) if(read(a,--i)) { n++; if(n == 0) {j=(Int) i; break;} } 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(a,0)); if(n0 != checkArray(read(b,0))) return false; for(size_t i=0; i < n; ++i) { array *ai=read(a,i); array *bi=read(b,i); for(size_t j=0; j < n0; ++j) { if(read(ai,j) != read(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(a,i) ? (*b)[i] : (*c)[i]; } else { r->clear(); if(b) { checkArrays(a,b); for(size_t i=0; i < size; i++) if(read(a,i)) r->push((*b)[i]); } else if(c) { checkArrays(a,c); for(size_t i=0; i < size; i++) if(!read(a,i)) r->push((*c)[i]); } } 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(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(a,j); if(pivot[j] != 1) { for(size_t k=0; k < n; k++) { if(pivot[k] == 0) { real temp=fabs(read(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(a,col); if(row != col) { array *arow=read(a,row); for(size_t l=0; l < n; l++) { real temp=read(arow,l); (*arow)[l]=read(acol,l); (*acol)[l]=temp; } } Row[i]=row; Col[i]=col; // Divide the pivot row by the pivot element. real denom=read(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(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(a,k); real akcol=read(ak,col); (*ak)[col]=0.0; for(size_t l=0; l < n; l++) (*ak)[l]=read(ak,l)-read(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(a,k); real temp=read(ak,r); (*ak)[r]=read(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(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(a,i); if(checkArray(ai) != m) error(incommensurate); real sum=0.0; for(size_t j=0; j < m; ++j) sum += read(ai,j)*B[j]; (*c)[i]=sum; } delete[] B; 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(a,0)); if(na0 != nb) error(incommensurate); size_t nb0=nb == 0 ? 0 : checkArray(read(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(t,0); array *t1=read(t,1); array *t3=read(t,3); if(checkArray(t0) != 4 || checkArray(t1) != 4 || checkArray(t3) != 4) error(incommensurate); real x=v.getx(); real y=v.gety(); real z=v.getz(); real f=read(t3,0)*x+read(t3,1)*y+read(t3,2)*z+ read(t3,3); if(f == 0.0) dividebyzero(); f=1.0/f; return pair((read(t0,0)*x+read(t0,1)*y+read(t0,2)*z+ read(t0,3))*f, (read(t1,0)*x+read(t1,1)*y+read(t1,2)*z+ read(t1,3))*f); } // 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(a,i)*read(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(a,0) == 0.0 && read(c,n-1) == 0.0) { real temp=read(b,0); if(temp == 0.0) dividebyzero(); temp=1.0/temp; real *work=new real[n]; u[0]=read(f,0)*temp; work[0]=-read(c,0)*temp; for(size_t i=1; i < n; i++) { real temp=(read(b,i)+read(a,i)*work[i-1]); if(temp == 0.0) {delete[] work; dividebyzero();} temp=1.0/temp; u[i]=(read(f,i)-read(a,i)*read(u,i-1))*temp; work[i]=-read(c,i)*temp; } for(size_t i=n-1; i >= 1; i--) u[i-1]=read(u,i-1)+work[i-1]*read(u,i); delete[] work; return up; } real binv=read(b,0); if(binv == 0.0) dividebyzero(); binv=1.0/binv; if(n == 1) {u[0]=read(f,0)*binv; return up;} if(n == 2) { real factor=(read(b,0)*read(b,1)- read(a,0)*read(c,1)); if(factor== 0.0) dividebyzero(); factor=1.0/factor; real temp=(read(b,0)*read(f,1)- read(c,1)*read(f,0))*factor; u[0]=(read(b,1)*read(f,0)- read(a,0)*read(f,1))*factor; u[1]=temp; return up; } real *gamma=new real[n-2]; real *delta=new real[n-2]; gamma[0]=read(c,0)*binv; delta[0]=read(a,0)*binv; u[0]=read(f,0)*binv; real beta=read(c,n-1); real fn=read(f,n-1)-beta*read(u,0); real alpha=read(b,n-1)-beta*delta[0]; for(size_t i=1; i <= n-3; i++) { real alphainv=read(b,i)-read(a,i)*gamma[i-1]; if(alphainv == 0.0) {delete[] gamma; delete[] delta; dividebyzero();} alphainv=1.0/alphainv; beta *= -gamma[i-1]; gamma[i]=read(c,i)*alphainv; u[i]=(read(f,i)-read(a,i)*read(u,i-1))*alphainv; fn -= beta*read(u,i); delta[i]=-read(a,i)*delta[i-1]*alphainv; alpha -= beta*delta[i]; } real alphainv=read(b,n-2)-read(a,n-2)*gamma[n-3]; if(alphainv == 0.0) {delete[] gamma; delete[] delta; dividebyzero();} alphainv=1.0/alphainv; u[n-2]=(read(f,n-2)-read(a,n-2)*read(u,n-3)) *alphainv; beta=read(a,n-1)-beta*gamma[n-3]; real dnm1=(read(c,n-2)-read(a,n-2)*delta[n-3])*alphainv; real temp=alpha-beta*dnm1; if(temp == 0.0) {delete[] gamma; delete[] delta; dividebyzero();} u[n-1]=temp=(fn-beta*read(u,n-2))/temp; u[n-2]=read(u,n-2)-dnm1*temp; for(size_t i=n-2; i >= 1; i--) u[i-1]=read(u,i-1)-gamma[i-1]*read(u,i)-delta[i-1]*temp; delete[] delta; delete[] gamma; 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(Stack); if(dfdx == 0.0) { x=DBL_MAX; break; } Stack->push(x); f->call(Stack); real fx=pop(Stack); x -= fx/dfdx; lastdiff=diff; if(verbose) cout << "Newton-Raphson: " << x << endl; diff=fabs(x-x0); if(++i == iterations) { x=DBL_MAX; break; } } while (diff != 0.0 && (diff < lastdiff || diff > fuzz*fabs(x))); if(verbose) cout.precision(oldPrec); 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(Stack); if(f1 == 0.0) return x1; Stack->push(x2); f->call(Stack); real f2=pop(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(Stack); Stack->push(x); fprime->call(Stack); real dy=pop(Stack); Int j; for(j=0; j < iterations; j++) { if(((x-x2)*dy-y)*((x-x1)*dy-y) >= 0.0 || fabs(2.0*y) > fabs(dxold*dy)) { dxold=dx; dx=0.5*(x2-x1); x=x1+dx; if(verbose) cout << "bisection: " << x << endl; if(x1 == x) 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(Stack); Stack->push(x); fprime->call(Stack); dy=pop(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(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(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 }