diff options
Diffstat (limited to 'Build/source/utils/asymptote/runtime.in')
-rw-r--r-- | Build/source/utils/asymptote/runtime.in | 5003 |
1 files changed, 1 insertions, 5002 deletions
diff --git a/Build/source/utils/asymptote/runtime.in b/Build/source/utils/asymptote/runtime.in index d4e80bd8b9a..08aabf2ea3e 100644 --- a/Build/source/utils/asymptote/runtime.in +++ b/Build/source/utils/asymptote/runtime.in @@ -66,21 +66,9 @@ stringarray2* => stringArray2() #include "path.h" #include "path3.h" #include "pen.h" +#include "drawpath.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" @@ -96,18 +84,7 @@ stringarray2* => stringArray2() #include "lexical.h" #include "process.h" #include "arrayop.h" -#include "predicates.h" -#include "Delaunay.h" -#ifdef HAVE_LIBFFTW3 -#include "fftw++.h" -#endif - -#if defined(HAVE_LIBREADLINE) && defined(HAVE_LIBCURSES) -#include <readline/readline.h> -#include <readline/history.h> -#endif - #if defined(USEGC) && defined(GC_DEBUG) && defined(GC_BACKTRACE) extern "C" { void *GC_generate_random_valid_address(void); @@ -127,10 +104,6 @@ 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; @@ -180,27 +153,11 @@ 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()); @@ -213,326 +170,9 @@ namespace vm { namespace run { -const char *invalidargument="invalid argument"; const char *arrayempty="cannot take min or max of empty array"; const char *noruntime="no runtime environment for embedded eval"; -// 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"; - error(buf); - } - return size; -} - -double *copyArrayC(array *a, size_t dim=0) -{ - size_t size=checkdimension(a,dim); - double *c=new double[size]; - for(size_t i=0; i < size; i++) - c[i]=read<double>(a,i); - return c; -} - -triple *copyTripleArrayC(array *a, size_t dim=0) -{ - size_t size=checkdimension(a,dim); - triple *c=new triple[size]; - for(size_t i=0; i < size; i++) - c[i]=read<triple>(a,i); - return c; -} - -array *copyArray2(array *a) -{ - size_t size=checkArray(a); - array *c=new array(size); - for(size_t i=0; i < size; i++) { - array *ai=read<array*>(a,i); - size_t aisize=checkArray(ai); - array *ci=new array(aisize); - (*c)[i]=ci; - for(size_t j=0; j < aisize; j++) - (*ci)[j]=(*ai)[j]; - } - return c; -} - -array *copyArray3(array *a) -{ - size_t size=checkArray(a); - array *c=new array(size); - for(size_t i=0; i < size; i++) { - array *ai=read<array*>(a,i); - size_t aisize=checkArray(ai); - array *ci=new array(aisize); - (*c)[i]=ci; - for(size_t j=0; j < aisize; j++) { - array *aij=read<array*>(ai,j); - size_t aijsize=checkArray(aij); - array *cij=new array(aijsize); - (*ci)[j]=cij; - for(size_t k=0; k < aijsize; k++) - (*cij)[k]=(*aij)[k]; - } - } - return c; -} - -double *copyArray2C(array *a, bool square=true, size_t dim2=0) -{ - size_t n=checkArray(a); - size_t m=(square || n == 0) ? n : checkArray(read<array*>(a,0)); - if(n > 0 && dim2 && m != dim2) { - ostringstream buf; - buf << "second matrix dimension must be " << dim2; - error(buf); - } - - double *c=new double[n*m]; - for(size_t i=0; i < n; i++) { - array *ai=read<array*>(a,i); - size_t aisize=checkArray(ai); - if(aisize == m) { - double *ci=c+i*m; - for(size_t j=0; j < m; j++) - ci[j]=read<double>(ai,j); - } else - error(square ? "matrix must be square" : "matrix must be rectangular"); - } - return c; -} - -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<array*>(a,0)); - if(n > 0 && dim2 && m != dim2) { - ostringstream buf; - buf << "second matrix dimension must be " << dim2; - error(buf); - } - - triple *c=new triple[n*m]; - for(size_t i=0; i < n; i++) { - array *ai=read<array*>(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<triple>(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<array*>(a,0)); - if(n > 0 && dim2 && m != dim2) { - ostringstream buf; - buf << "second matrix dimension must be " << dim2; - 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<array*>(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<triple>(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; -} - -static const char *incommensurate="Incommensurate matrices"; -static const char *singular="Singular matrix"; -static size_t *pivot,*Row,*Col; - -triple operator *(const array& t, const triple& v) -{ - size_t n=checkArray(&t); - if(n != 4) error(incommensurate); - array *t0=read<array*>(t,0); - array *t1=read<array*>(t,1); - array *t2=read<array*>(t,2); - array *t3=read<array*>(t,3); - - if(checkArray(t0) != 4 || checkArray(t1) != 4 || - checkArray(t2) != 4 || checkArray(t3) != 4) - error(incommensurate); - - double x=v.getx(); - double y=v.gety(); - double z=v.getz(); - - double f=read<real>(t3,0)*x+read<real>(t3,1)*y+read<real>(t3,2)*z+ - read<real>(t3,3); - if(f == 0.0) run::dividebyzero(); - f=1.0/f; - - return triple((read<real>(t0,0)*x+read<real>(t0,1)*y+read<real>(t0,2)*z+ - read<real>(t0,3))*f, - (read<real>(t1,0)*x+read<real>(t1,1)*y+read<real>(t1,2)*z+ - read<real>(t1,3))*f, - (read<real>(t2,0)*x+read<real>(t2,1)*y+read<real>(t2,2)*z+ - read<real>(t2,3))*f); -} - -triple multshiftless(const array& t, const triple& v) -{ - size_t n=checkArray(&t); - if(n != 4) error(incommensurate); - array *t0=read<array*>(t,0); - array *t1=read<array*>(t,1); - array *t2=read<array*>(t,2); - array *t3=read<array*>(t,3); - - if(checkArray(t0) != 4 || checkArray(t1) != 4 || - checkArray(t2) != 4 || checkArray(t3) != 4) - error(incommensurate); - - double x=v.getx(); - double y=v.gety(); - double z=v.getz(); - - double f=read<real>(t3,0)*x+read<real>(t3,1)*y+read<real>(t3,2)*z+ - read<real>(t3,3); - if(f == 0.0) run::dividebyzero(); - f=1.0/f; - - return triple((read<real>(t0,0)*x+read<real>(t0,1)*y+read<real>(t0,2)*z)*f, - (read<real>(t1,0)*x+read<real>(t1,1)*y+read<real>(t1,2)*z)*f, - (read<real>(t2,0)*x+read<real>(t2,1)*y+read<real>(t2,2)*z)*f); -} - -static inline void inverseAllocate(size_t n) -{ - pivot=new size_t[n]; - Row=new size_t[n]; - Col=new size_t[n]; -} - -static inline void inverseDeallocate() -{ - delete[] pivot; - delete[] Row; - delete[] Col; -} - void writestring(stack *s) { callable *suffix=pop<callable *>(s,NULL); @@ -550,547 +190,18 @@ void writestring(stack *s) } } -void checkSquare(array *a) -{ - size_t n=checkArray(a); - for(size_t i=0; i < n; i++) - if(checkArray(read<array*>(a,i)) != n) - error("matrix a must be square"); -} - -// Crout's algorithm for computing the LU decomposition of a square matrix. -// cf. routine ludcmp (Press et al., Numerical Recipes, 1991). -Int LUdecompose(double *a, size_t n, size_t* index, bool warn=true) -{ - double *vv=new double[n]; - Int swap=1; - for(size_t i=0; i < n; ++i) { - double big=0.0; - double *ai=a+i*n; - for(size_t j=0; j < n; ++j) { - double temp=fabs(ai[j]); - if(temp > big) big=temp; - } - if(big == 0.0) { - delete[] vv; - if(warn) error(singular); - else return 0; - } - vv[i]=1.0/big; - } - for(size_t j=0; j < n; ++j) { - for(size_t i=0; i < j; ++i) { - double *ai=a+i*n; - double sum=ai[j]; - for(size_t k=0; k < i; ++k) { - sum -= ai[k]*a[k*n+j]; - } - ai[j]=sum; - } - double big=0.0; - size_t imax=j; - for(size_t i=j; i < n; ++i) { - double *ai=a+i*n; - double sum=ai[j]; - for(size_t k=0; k < j; ++k) - sum -= ai[k]*a[k*n+j]; - ai[j]=sum; - double temp=vv[i]*fabs(sum); - if(temp >= big) { - big=temp; - imax=i; - } - } - double *aj=a+j*n; - double *aimax=a+imax*n; - if(j != imax) { - for(size_t k=0; k < n; ++k) { - double temp=aimax[k]; - aimax[k]=aj[k]; - aj[k]=temp; - } - swap *= -1; - vv[imax]=vv[j]; - } - if(index) - index[j]=imax; - if(j != n) { - double denom=aj[j]; - if(denom == 0.0) { - delete[] vv; - if(warn) error(singular); - else return 0; - } - for(size_t i=j+1; i < n; ++i) - a[i*n+j] /= denom; - } - } - delete[] vv; - return swap; -} - -void dividebyzero(size_t i) -{ - ostringstream buf; - if(i > 0) buf << "array element " << i << ": "; - buf << "Divide by zero"; - error(buf); -} - -void integeroverflow(size_t i) -{ - ostringstream buf; - if(i > 0) buf << "array element " << i << ": "; - buf << "Integer overflow"; - error(buf); -} - -#if defined(HAVE_LIBREADLINE) && defined(HAVE_LIBCURSES) -struct historyState { - bool store; - HISTORY_STATE state; -}; - -typedef mem::map<CONST string, historyState> historyMap_t; -historyMap_t historyMap; -static HISTORY_STATE history_save; - -// Store a deep copy of the current readline history in dest. -void store_history(HISTORY_STATE *dest) -{ - HISTORY_STATE *src=history_get_history_state(); - if(src) { - *dest=*src; - for(Int i=0; i < src->length; ++i) - dest->entries[i]=src->entries[i]; - free(src); - } -} - -stringarray* get_history(Int n) -{ - int N=intcast(n); - if(N <= 0) N=history_length; - else N=Min(N,history_length); - array *a=new array((size_t) N); - int offset=history_length-N+1; - for(int i=0; i < N; ++i) { - HIST_ENTRY *last=history_get(offset+i); - string s=last ? last->line : ""; - (*a)[i]=s; - } - return a; -} - -string historyfilename(const string &name) -{ - return historyname+"_"+name; -} -#endif - -#if defined(HAVE_LIBREADLINE) && defined(HAVE_LIBCURSES) -int readline_startup_hook() -{ -#ifdef __CYGWIN__ - rl_set_key("\\M-[3~",rl_delete,rl_get_keymap()); - rl_set_key("\\M-[2~",rl_overwrite_mode,rl_get_keymap()); -#endif - return 0; -} - -void init_readline(bool tabcompletion=true) -{ - static bool first=true; - if(first) { - first=false; -#ifdef __CYGWIN__ - rl_startup_hook=readline_startup_hook; -#endif - } - rl_bind_key('\t',tabcompletion ? rl_complete : rl_insert); -} -#endif - -void cleanup() -{ -#if defined(HAVE_LIBREADLINE) && defined(HAVE_LIBCURSES) - store_history(&history_save); - int nlines=intcast(getSetting<Int>("historylines")); - for(historyMap_t::iterator h=historyMap.begin(); h != historyMap.end(); - ++h) { - history_set_history_state(&h->second.state); - stifle_history(nlines); - if(h->second.store) write_history(historyfilename(h->first).c_str()); - } - history_set_history_state(&history_save); -#endif -} - -void purge(Int divisor=0) -{ -#ifdef USEGC - if(divisor > 0) GC_set_free_space_divisor((GC_word) divisor); - GC_gcollect(); -#endif -} - -void updateFunction(stack *Stack) -{ - callable *atUpdateFunction=processData().atUpdateFunction; - if(atUpdateFunction && !nullfunc::instance()->compare(atUpdateFunction)) - atUpdateFunction->call(Stack); -} - -void exitFunction(stack *Stack) -{ - callable *atExitFunction=processData().atExitFunction; - if(atExitFunction && !nullfunc::instance()->compare(atExitFunction)) - atExitFunction->call(Stack); -} - default_t def; string emptystring; -array *emptyarray=new array(0); -string commentchar="#"; pair zero; -void breakpoint(stack *Stack, runnable *r) -{ - callable *atBreakpointFunction=processData().atBreakpointFunction; - if(atBreakpointFunction && - !nullfunc::instance()->compare(atBreakpointFunction)) { - position curPos=getPos(); - Stack->push<string>(curPos.filename()); - Stack->push<Int>((Int) curPos.Line()); - Stack->push<Int>((Int) curPos.Column()); - Stack->push(r ? r : item(def)); - atBreakpointFunction->call(Stack); // returns a string - } else Stack->push<string>(""); -} - -} - -namespace types { -extern const char *names[]; -} - -void checkformat(const char *ptr, bool intformat) -{ - while(*ptr != '\0') { - if(*ptr != '%') /* While we have regular characters, print them. */ - ptr++; - else { /* We've got a format specifier. */ - ptr++; - - while(*ptr && strchr ("-+ #0'I", *ptr)) /* Move past flags. */ - *ptr++; - - if(*ptr == '*') - ptr++; - else while(isdigit(*ptr)) /* Handle explicit numeric value. */ - ptr++; - - if(*ptr == '.') { - *ptr++; /* Go past the period. */ - if(*ptr == '*') { - ptr++; - } else - while(isdigit(*ptr)) /* Handle explicit numeric value. */ - *ptr++; - } - while(*ptr && strchr ("hlL", *ptr)) - *ptr++; - - if(*ptr == '%') ++ptr; - else if(*ptr != '\0') { - if(intformat) { - switch(*ptr) { - case 'd': - case 'i': - case 'o': - case 'u': - case 'x': - case 'X': - case 'c': - break; - default: - ostringstream buf; - buf << "Invalid format '" << *ptr << "' for type " - << types::names[types::ty_Int]; - error(buf); - break; - } - } else { - switch(*ptr) { - case 'f': - case 'F': - case 'e': - case 'E': - case 'g': - case 'G': - break; - default: - ostringstream buf; - buf << "Invalid format '" << *ptr << "' for type " - << types::names[types::ty_real]; - error(buf); - break; - } - } - } - } /* End of else statement */ - } -} - -// Return an angle in the interval [0,360). -inline double principalBranch(double deg) -{ - if(deg < 0) deg += 360; - return deg; } static string defaulttransparency=string("Compatible"); -static const string defaulttimeformat=string("%a %b %d %T %Z %Y"); -#ifdef HAVE_STRFTIME -static const size_t nTime=256; -static char Time[nTime]; -#endif - -void clear(string file, Int line, bool warn=false) -{ - bpinfo bp(file,line); - for(mem::list<bpinfo>::iterator p=bplist.begin(); p != bplist.end(); ++p) { - if(*p == bp) { - cout << "cleared breakpoint at " << file << ": " << line << endl; - bplist.remove(bp); - return; - } - } - if(warn) - cout << "No such breakpoint at " << file << ": " << line << endl; -} - -Int windingnumber(array *p, camp::pair z) -{ - size_t size=checkArray(p); - Int count=0; - for(size_t i=0; i < size; i++) - count += read<path *>(p,i)->windingnumber(z); - return count; -} - -string convertname(string name, const string& format, bool check=true) -{ - if(name.empty()) - return buildname(outname(),format,""); - else - if(check) checkLocal(name); - return format.empty() ? name : format+":"+name; -} - -callable *Func; -stack *FuncStack; -double wrapFunction(double x) -{ - FuncStack->push(x); - Func->call(FuncStack); - return pop<double>(FuncStack); -} - -callable *compareFunc; -bool compareFunction(const vm::item& i, const vm::item& j) -{ - FuncStack->push(i); - FuncStack->push(j); - compareFunc->call(FuncStack); - return pop<bool>(FuncStack); -} void unused(void *) { } -pair readpair(stringstream& s, 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 ASYx="/ASYx {( ) print ASYX sub 12 string cvs print} bind def"; -string ASYy="/ASYy {( ) print ASYY sub 12 string cvs print} bind def"; -string pathforall="{(M) print ASYy ASYx} {(L) print ASYy ASYx} {(C) print ASYy ASYx ASYy ASYx ASYy ASYx} {(c) print} pathforall"; -string currentpoint="print currentpoint ASYy ASYx "; -string ASYinit="/ASYX currentpoint pop def /ASYY currentpoint exch pop def "; -string ASY1="ASY1 {"+ASYinit+"/ASY1 false def} if "; - -void showpath(std::ofstream& ps) -{ - ps << ASYx << newl - << ASYy << newl - << "/ASY1 true def" << newl - << "/fill {" << ASY1 - << pathforall << " (M) " << currentpoint - << "currentpoint newpath moveto } bind def" << newl - << "/stroke {" << ASY1 << "strokepath " - << pathforall << " (M) " << currentpoint - << "currentpoint newpath moveto } bind def" << endl; -} - -array *readpath(const string& psname, bool keep, - double hscale=1.0, double vsign=1.0) -{ - double vscale=vsign*hscale; - array *P=new array(0); - mem::vector<string> cmd; - cmd.push_back(getSetting<string>("gs")); - cmd.push_back("-q"); - cmd.push_back("-dNOPAUSE"); - cmd.push_back("-dBATCH"); - if(safe) cmd.push_back("-dSAFER"); -#ifdef __CYGWIN__ - const string null="NUL"; -#else - const string null="/dev/null"; -#endif - cmd.push_back("-sDEVICE=epswrite"); - cmd.push_back("-sOutputFile="+null); - cmd.push_back(psname); - iopipestream gs(cmd,"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<solvedKnot> nodes; - solvedKnot node; - bool cyclic=false; - bool active=false; - - 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); - node.straight=false; - break; - } - case 'L': - { - pair point=readpair(buf,hscale,vscale); - 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); - pair pre=readpair(buf,hscale,vscale); - node.post=readpair(buf,hscale,vscale); - node.straight=false; - nodes.push_back(node); - active=true; - node.pre=pre; - node.point=point; - break; - } - case 'c': - { - cyclic=true; - 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; - error(buf); -} - -void cannotwrite(const string& s) -{ - ostringstream buf; - buf << "Cannot write to " << s; - error(buf); -} - // Autogenerated routines: @@ -1145,57 +256,6 @@ bool :isDefault(item 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) @@ -1214,689 +274,6 @@ path :guideToPath(guide *g) } -// 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(); -} - -realarray *texsize(string *s, pen p=CURRENTPEN) -{ - texinit(); - processDataStruct &pd=processData(); - - string texengine=getSetting<string>("tex"); - const char **abort=texabort(texengine); - setpen(pd.tex,texengine,p); - - double width,height,depth; - if(!texbounds(width,height,depth,pd.tex,*s,abort,false,true)) - return new array(0); - - array *t=new array(3); - (*t)[0]=width; - (*t)[1]=height; - (*t)[2]=depth; - return t; -} - -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; - string texengine=getSetting<string>("tex"); - bool pdf=settings::pdf(texengine); - texfile tex(texname,b,true); - tex.miniprologue(); - tex.setfont(p); - - if(!pdf) { - tex.verbatimline("\\special{ps:"); - tex.verbatimline(ASYx); - tex.verbatimline(ASYy); - tex.verbatimline("/ASY1 true def"); - tex.verbatimline("/v {"+ASY1+"neg exch 4 copy 4 2 roll 2 copy 6 2 roll 2 copy (M) print ASYy ASYx (L) print ASYy add ASYx (L) print add ASYy add ASYx (L) print add ASYy ASYx (c) print} bind def"); - tex.verbatimline("/show {"+ASY1+ - "currentpoint newpath moveto false charpath "+ - pathforall+"} bind def}"); - } - tex.verbatimline(*s+"%"); - tex.epilogue(true); - tex.close(); - - int status=opentex(texname,prefix); - - string pdfname; - if(!status) { - if(pdf) { - pdfname=auxname(prefix,"pdf"); - std::ofstream ps(psname.c_str()); - if(!ps) cannotwrite(psname); - - showpath(ps); - - mem::vector<string> cmd; - cmd.push_back(getSetting<string>("gs")); - cmd.push_back("-q"); - cmd.push_back("-dNOCACHE"); - cmd.push_back("-dNOPAUSE"); - cmd.push_back("-dBATCH"); - if(safe) cmd.push_back("-dSAFER"); - cmd.push_back("-sDEVICE=epswrite"); - cmd.push_back("-sOutputFile=-"); - cmd.push_back(pdfname); - iopipestream gs(cmd,"gs","Ghostscript"); - - gs.block(false); - while(true) { - string line; - gs >> line; - if(line.empty() && !gs.running()) break; - ps << line; - } - ps.close(); - } else { - mem::vector<string> dcmd; - dcmd.push_back(getSetting<string>("dvips")); - dcmd.push_back("-R"); - dcmd.push_back("-Pdownload35"); - dcmd.push_back("-D600"); - push_split(dcmd,getSetting<string>("dvipsOptions")); - if(verbose <= 2) dcmd.push_back("-q"); - dcmd.push_back("-o"+psname); - dcmd.push_back(dviname); - status=System(dcmd,0,true,"dvips"); - } - } - - if(status != 0) - error("texpath failed"); - - bool keep=getSetting<bool>("keep"); - if(!keep) { // Delete temporary files. - unlink(texname.c_str()); - if(!getSetting<bool>("keepaux")) - unlink(auxname(prefix,"aux").c_str()); - unlink(auxname(prefix,"log").c_str()); - if(pdf) - unlink(pdfname.c_str()); - else - unlink(dviname.c_str()); - if(settings::context(texengine)) { - unlink(auxname(prefix,"top").c_str()); - unlink(auxname(prefix,"tua").c_str()); - unlink(auxname(prefix,"tui").c_str()); - } - } - return pdf ? readpath(psname,keep,0.1) : readpath(psname,keep,0.12,-1.0); -} - -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<string>("textoutformat")); - - string textname=auxname(prefix,getSetting<string>("textextension")); - std::ofstream text(textname.c_str()); - - if(!text) cannotwrite(textname); - - text << getSetting<string>("textprologue") << newl - << p.Font() << newl - << *s << newl - << getSetting<string>("textepilogue") << endl; - text.close(); - - string psname=auxname(prefix,"ps"); - std::ofstream ps(psname.c_str()); - if(!ps) cannotwrite(psname); - - showpath(ps); - - mem::vector<string> cmd; - cmd.push_back(getSetting<string>("textcommand")); - push_split(cmd,getSetting<string>("textcommandOptions")); - cmd.push_back(textname); - iopipestream typesetter(cmd); - - mem::vector<string> cmd2; - cmd2.push_back(getSetting<string>("gs")); - cmd2.push_back("-q"); - cmd2.push_back("-dNOCACHE"); - cmd2.push_back("-dNOPAUSE"); - cmd2.push_back("-dBATCH"); - if(safe) cmd2.push_back("-dSAFER"); - cmd2.push_back("-sDEVICE=epswrite"); - cmd2.push_back("-sOutputFile=-"); - cmd2.push_back("-"); - iopipestream gs(cmd2,"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<bool>("keep"); - if(!keep) // Delete temporary files. - unlink(textname.c_str()); - return readpath(psname,keep,0.1); -} - -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(ASYx); - ps.verbatimline(ASYy); - ps.verbatimline("/stroke {"+ASYinit+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<bool>("keep")); -} - -void _draw(picture *f, path g, pen p) -{ - f->append(new drawPath(g,p)); -} - -void fill(picture *f, patharray *g, pen p=CURRENTPEN, bool copy=true) -{ - array *(*copyarray)(array *a)=copy ? copyArray: nop; - f->append(new drawFill(*copyarray(g),false,p)); -} - -void latticeshade(picture *f, patharray *g, bool stroke=false, - pen fillrule=CURRENTPEN, penarray2 *p, bool copy=true) -{ - array *(*copyarray)(array *a)=copy ? copyArray: nop; - f->append(new drawLatticeShade(*copyarray(g),stroke,fillrule,*copyarray(p))); -} - -void axialshade(picture *f, patharray *g, bool stroke=false, pen pena, pair a, - pen penb, pair b, bool copy=true) -{ - array *(*copyarray)(array *a)=copy ? copyArray: nop; - f->append(new drawAxialShade(*copyarray(g),stroke,pena,a,penb,b)); -} - -void radialshade(picture *f, patharray *g, bool stroke=false, pen pena, - pair a, real ra, pen penb, pair b, real rb, bool copy=true) -{ - array *(*copyarray)(array *a)=copy ? copyArray: nop; - f->append(new drawRadialShade(*copyarray(g),stroke,pena,a,ra,penb,b,rb)); -} - -void gouraudshade(picture *f, patharray *g, bool stroke=false, - pen fillrule=CURRENTPEN, penarray *p, pairarray *z, - Intarray *edges, bool copy=true) -{ - array *(*copyarray)(array *a)=copy ? copyArray: nop; - checkArrays(p,z); - checkArrays(z,edges); - f->append(new drawGouraudShade(*copyarray(g),stroke,fillrule,*copyarray(p), - *copyarray(z),*copyarray(edges))); -} - -void gouraudshade(picture *f, patharray *g, bool stroke=false, - pen fillrule=CURRENTPEN, penarray *p, Intarray *edges, - bool copy=true) -{ - array *(*copyarray)(array *a)=copy ? copyArray: nop; - size_t n=checkArrays(p,edges); - size_t m=checkArray(g); - array *z=new array(n); - Int k=0; - Int in=(Int) n; - for(size_t j=0; j < m; ++j) { - path *P=read<path *>(g,j); - assert(P); - Int stop=Min(P->size(),in-k); - mem::vector<solvedKnot>& nodes=P->Nodes(); - for(Int i=0; i < stop; ++i) - (*z)[k++]=nodes[i].point; - } - checkArrays(p,z); - - f->append(new drawGouraudShade(*copyarray(g),stroke,fillrule,*copyarray(p), - *z,*copyarray(edges))); -} - -void tensorshade(picture *f, patharray *g, bool stroke=false, - pen fillrule=CURRENTPEN, penarray2 *p, patharray *b=NULL, - pairarray2 *z=emptyarray, bool copy=true) -{ - array *(*copyarray)(array *a)=copy ? copyArray: nop; - array *(*copyarray2)(array *a)=copy ? copyArray2: nop; - if(b == NULL) b=g; - size_t n=checkArrays(p,b); - size_t nz=checkArray(z); - if(nz != 0) - checkEqual(nz,n); - f->append(new drawTensorShade(*copyarray(g),stroke,fillrule,*copyarray2(p), - *copyarray(b),*copyarray2(z))); -} - -void functionshade(picture *f, patharray *g, bool stroke=false, - pen fillrule=CURRENTPEN, string shader=emptystring, - bool copy=true) -{ - array *(*copyarray)(array *a)=copy ? copyArray: nop; - f->append(new drawFunctionShade(*copyarray(g),stroke,fillrule,shader)); -} - -// Clip a picture to a superpath using the given fill rule. -// Subsequent additions to the picture will not be affected by the clipping. -void clip(picture *f, patharray *g, bool stroke=false, - pen fillrule=CURRENTPEN, bool copy=true) -{ - array *(*copyarray)(array *a)=copy ? copyArray: nop; - drawClipBegin *begin=new drawClipBegin(*copyarray(g),stroke,fillrule,true); - f->enclose(begin,new drawClipEnd(true,begin)); -} - -void beginclip(picture *f, patharray *g, bool stroke=false, - pen fillrule=CURRENTPEN, bool copy=true) -{ - array *(*copyarray)(array *a)=copy ? copyArray: nop; - f->append(new drawClipBegin(*copyarray(g),stroke,fillrule,false)); -} - -void endclip(picture *f) -{ - f->append(new drawClipEnd(false)); -} - -void gsave(picture *f) -{ - f->append(new drawGsave()); -} - -void grestore(picture *f) -{ - f->append(new drawGrestore()); -} - -void begingroup(picture *f) -{ - f->append(new drawBegin()); -} - -void endgroup(picture *f) -{ - f->append(new drawEnd()); -} - -void add(picture *dest, picture *src) -{ - dest->add(*src); -} - -void prepend(picture *dest, picture *src) -{ - dest->prepend(*src); -} - -void postscript(picture *f, string s) -{ - f->append(new drawVerbatim(PostScript,s)); -} - -void tex(picture *f, string s) -{ - f->append(new drawVerbatim(TeX,s)); -} - -void postscript(picture *f, string s, pair min, pair max) -{ - f->append(new drawVerbatim(PostScript,s,min,max)); -} - -void tex(picture *f, string s, pair min, pair max) -{ - f->append(new drawVerbatim(TeX,s,min,max)); -} - -void texpreamble(string s) -{ - string t=s+"\n"; - processDataStruct &pd=processData(); - pd.TeXpipepreamble.push_back(t); - pd.TeXpreamble.push_back(t); -} - -void deletepreamble() -{ - if(getSetting<bool>("inlinetex")) { - unlink(auxname(outname(),"pre").c_str()); - } -} - -void _labelpath(picture *f, string s, string size, path g, string justify, - pair offset, pen p) -{ - f->append(new drawLabelPath(s,size,g,justify,offset,p)); -} - -void texreset() -{ - processDataStruct &pd=processData(); - pd.TeXpipepreamble.clear(); - pd.TeXpreamble.clear(); - pd.tex.pipeclose(); -} - -void layer(picture *f) -{ - f->append(new drawLayer()); -} - -void newpage(picture *f) -{ - f->append(new drawNewPage()); -} - -void _image(picture *f, realarray2 *data, pair initial, pair final, - penarray *palette=NULL, transform t=identity, bool copy=true, - bool antialias=false) -{ - array *(*copyarray)(array *a)=copy ? copyArray: nop; - array *(*copyarray2)(array *a)=copy ? copyArray2: nop; - f->append(new drawImage(*copyarray2(data),*copyarray(palette), - t*matrix(initial,final),antialias)); -} - -void _image(picture *f, penarray2 *data, pair initial, pair final, - transform t=identity, bool copy=true, bool antialias=false) -{ - array *(*copyarray2)(array *a)=copy ? copyArray2: nop; - f->append(new drawImage(*copyarray2(data),t*matrix(initial,final),antialias)); -} - -string nativeformat() -{ - return nativeformat(); -} - -bool latex() -{ - return latex(getSetting<string>("tex")); -} - -bool pdf() -{ - return pdf(getSetting<string>("tex")); -} - -void shipout(string prefix=emptystring, picture *f, picture *preamble=NULL, - string format=emptystring, bool wait=false, bool view=true, - callableTransform *xform) -{ - if(prefix.empty()) prefix=outname(); - - picture *result=new picture; - unsigned level=0; - picture::nodelist::iterator p; - for(p = f->nodes.begin(); p != f->nodes.end(); ++p) { - xform->call(Stack); - transform t=pop<transform>(Stack); - static transform Zero=transform(0.0,0.0,0.0,0.0,0.0,0.0); - bool Delete=(t == Zero); - picture *group=new picture; - assert(*p); - if((*p)->endgroup()) error(nobegin); - if((*p)->begingroup()) { - ++level; - while(p != f->nodes.end() && level) { - if(!Delete) { - drawElement *e=t.isIdentity() ? *p : (*p)->transformed(t); - group->append(e); - } - ++p; - if(p == f->nodes.end()) break; - assert(*p); - if((*p)->begingroup()) ++level; - if((*p)->endgroup()) { - if(level) --level; - else error(nobegin); - } - } - } - if(p == f->nodes.end()) break; - assert(*p); - if(!Delete) { - drawElement *e=t.isIdentity() ? *p : (*p)->transformed(t); - group->append(e); - result->add(*group); - } - } - - result->shipout(preamble,prefix,format,0.0,wait,view); -} - -void shipout3(string prefix, picture *f, string format=emptystring, - real width, real height, real angle, real zoom, - triple m, triple M, pair shift, 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,zoom,m,M,shift, - 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(); - const string xformat="png"; - - static long arg_max=sysconf(_SC_ARG_MAX); - const unsigned maxargs=::min(arg_max/(prefix.size()+xformat.size()+25ul), - 256ul); - - cout << maxargs << newl; - - string preformat=nativeformat(); - const string Done="Done"; - const string Error="Error"; - - mem::vector<string> 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; - - if(magnification > 0.0) { - mem::list<string> nameStack; - string outname; - unsigned arg=0; - unsigned batch=0; - for(picture::nodelist::iterator p=f->nodes.begin();;) { - if(p == f->nodes.end()) break; - if(arg == 0) { - cmd.clear(); - ostringstream buf; - buf << batch << "_"; - outname=buildname(prefix+buf.str()+"%d",xformat,""); - converter="gs"; - hint="Ghostscript"; - cmd.push_back(getSetting<string>(converter)); - cmd.push_back("-q"); - cmd.push_back("-dNOPAUSE"); - cmd.push_back("-dBATCH"); - cmd.push_back("-sDEVICE=pngalpha"); - cmd.push_back("-dEPSCrop"); - if(safe) - cmd.push_back("-dSAFER"); - cmd.push_back("-r"+String(res)+"x"+String(res)); - cmd.push_back("-sOutputFile="+outname); - } - - picture *group=new picture; - xform->call(Stack); - transform t=pop<transform>(Stack); - assert(*p); - if((*p)->endgroup()) { - cout << Error << endl; - error(nobegin); - } - if((*p)->begingroup()) { - ++level; - while(p != f->nodes.end() && level) { - drawElement *e=t.isIdentity() ? *p : (*p)->transformed(t); - group->append(e); - ++p; - if(p == f->nodes.end()) break; - assert(*p); - if((*p)->begingroup()) ++level; - if((*p)->endgroup()) { - if(level) --level; - else { - cout << Error << endl; - error(nobegin); - } - } - } - } - if(p != f->nodes.end()) { - assert(*p); - drawElement *e=t.isIdentity() ? *p : (*p)->transformed(t); - group->append(e); - bbox b; - ostringstream buf; - buf << prefix << "_" << n; - group->shipout(preamble,buf.str(),preformat,magnification,false,false); - string Preformat=group->Transparency() ? "pdf" : preformat; - string name=buildname(buf.str(),Preformat); - nameStack.push_back(name); - cmd.push_back(name); - b=group->bounds(); - b *= magnification; - - cout << b << newl; - ++n; - ++p; - ++arg; - } - - if(p == f->nodes.end() || arg >= maxargs) { - arg=0; - ++batch; - cout.flush(); - - int status=System(cmd,0,true,converter,hint); - if(status) { - cout << Error << endl; - error("deconstruct failed"); - } - } - } - - if(!getSetting<bool>("keep")) { - for(mem::list<string>::iterator p=nameStack.begin(); - p != nameStack.end(); ++p) - unlink(p->c_str()); - } - - cout << Done << endl; - } -} - -void purge(Int divisor=0) -{ - purge(divisor); -} - - // Pen operations pen :newPen() @@ -2349,589 +726,6 @@ void list(string *s, bool imports=false) } -// Path operations - -path :nullPath() -{ - return nullpath; -} - -bool ==(path a, path b) -{ - return a == b; -} - -bool !=(path a, path b) -{ - return !(a == b); -} - -pair point(path p, Int t) -{ - return p.point((Int) t); -} - -pair point(path p, real t) -{ - return p.point(t); -} - -pair precontrol(path p, Int t) -{ - return p.precontrol((Int) t); -} - -pair precontrol(path p, real t) -{ - return p.precontrol(t); -} - -pair postcontrol(path p, Int t) -{ - return p.postcontrol((Int) t); -} - -pair postcontrol(path p, real t) -{ - return p.postcontrol(t); -} - -pair dir(path p, Int t, Int sign=0, bool normalize=true) -{ - return p.dir(t,sign,normalize); -} - -pair dir(path p, real t, bool normalize=true) -{ - return p.dir(t,normalize); -} - -pair accel(path p, Int t, Int sign=0) -{ - return p.accel(t,sign); -} - -pair accel(path p, real t) -{ - return p.accel(t); -} - -real radius(path p, real t) -{ - pair v=p.dir(t,false); - pair a=p.accel(t); - real d=dot(a,v); - real v2=v.abs2(); - real a2=a.abs2(); - real denom=v2*a2-d*d; - real r=v2*sqrt(v2); - return denom > 0 ? r/sqrt(denom) : 0.0; -} - -path reverse(path p) -{ - return p.reverse(); -} - -path subpath(path p, Int a, Int b) -{ - return p.subpath((Int) a, (Int) b); -} - -path subpath(path p, real a, real b) -{ - return p.subpath(a,b); -} - -path nurb(pair z0, pair z1, pair z2, pair z3, - real w0, real w1, real w2, real w3, Int m) -{ - return nurb(z0,z1,z2,z3,w0,w1,w2,w3,m); -} - -Int length(path p) -{ - return p.length(); -} - -bool cyclic(path p) -{ - return p.cyclic(); -} - -bool straight(path p, Int t) -{ - return p.straight(t); -} - -path unstraighten(path p) -{ - return p.unstraighten(); -} - -bool piecewisestraight(path p) -{ - return p.piecewisestraight(); -} - -real arclength(path p) -{ - return p.arclength(); -} - -real arctime(path p, real dval) -{ - return p.arctime(dval); -} - -real dirtime(path p, pair z) -{ - return p.directiontime(z); -} - -realarray* intersect(path p, path q, real fuzz=-1) -{ - bool exact=fuzz <= 0.0; - if(fuzz < 0) - fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())), - ::max(length(q.max()),length(q.min()))); - std::vector<real> S,T; - real s,t; - if(intersections(s,t,S,T,p,q,fuzz,true,exact)) { - array *V=new array(2); - (*V)[0]=s; - (*V)[1]=t; - return V; - } - return new array(0); -} - -realarray2* intersections(path p, path q, real fuzz=-1) -{ - bool exact=fuzz <= 0.0; - if(fuzz < 0.0) - fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())), - ::max(length(q.max()),length(q.min()))); - real s,t; - std::vector<real> S,T; - intersections(s,t,S,T,p,q,fuzz,false,true); - size_t n=S.size(); - if(n == 0 && !exact) { - if(intersections(s,t,S,T,p,q,fuzz,true,false)) { - array *V=new array(1); - array *Vi=new array(2); - (*V)[0]=Vi; - (*Vi)[0]=s; - (*Vi)[1]=t; - return V; - } - } - array *V=new array(n); - for(size_t i=0; i < n; ++i) { - array *Vi=new array(2); - (*V)[i]=Vi; - (*Vi)[0]=S[i]; - (*Vi)[1]=T[i]; - } - stable_sort(V->begin(),V->end(),run::compare2<real>()); - return V; -} - -realarray* intersections(path p, explicit pair a, explicit pair b, real fuzz=-1) -{ - if(fuzz < 0) - fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())), - ::max(length(a),length(b))); - std::vector<real> S; - intersections(S,p,a,b,fuzz); - sort(S.begin(),S.end()); - size_t n=S.size(); - array *V=new array(n); - for(size_t i=0; i < n; ++i) - (*V)[i]=S[i]; - return V; -} - -// Return the intersection point of the extensions of the line segments -// PQ and pq. -pair extension(pair P, pair Q, pair p, pair q) -{ - pair ac=P-Q; - pair bd=q-p; - real det=ac.getx()*bd.gety()-ac.gety()*bd.getx(); - if(det == 0) return pair(infinity,infinity); - return P+((p.getx()-P.getx())*bd.gety()-(p.gety()-P.gety())*bd.getx())*ac/det; -} - -Int size(path p) -{ - return p.size(); -} - -path &(path p, path q) -{ - return camp::concat(p,q); -} - -pair min(path p) -{ - return p.min(); -} - -pair max(path p) -{ - return p.max(); -} - -realarray *mintimes(path p) -{ - array *V=new array(2); - pair z=p.mintimes(); - (*V)[0]=z.getx(); - (*V)[1]=z.gety(); - return V; -} - -realarray *maxtimes(path p) -{ - array *V=new array(2); - pair z=p.maxtimes(); - (*V)[0]=z.getx(); - (*V)[1]=z.gety(); - return V; -} - -real relativedistance(real theta, real phi, real t, bool atleast) -{ - return camp::velocity(theta,phi,tension(t,atleast)); -} - -Int windingnumber(patharray *p, pair z) -{ - return windingnumber(p,z); -} - -bool inside(explicit patharray *g, pair z, pen fillrule=CURRENTPEN) -{ - return fillrule.inside(windingnumber(g,z)); -} - -bool inside(path g, pair z, pen fillrule=CURRENTPEN) -{ - return fillrule.inside(g.windingnumber(z)); -} - -// Determine the side of a--b that c lies on -// (negative=left, zero=on line, positive=right). -real side(pair a, pair b, pair c) -{ - return orient2d(a,b,c); -} - -// Determine the side of the counterclockwise circle through a,b,c that d -// lies on (negative=inside, 0=on circle, positive=right). -real incircle(pair a, pair b, pair c, pair d) -{ - return incircle(a.getx(),a.gety(),b.getx(),b.gety(),c.getx(),c.gety(), - d.getx(),d.gety()); -} - - -// Path3 operations - -path3 path3(triplearray *pre, triplearray *point, triplearray *post, - boolarray *straight, bool cyclic) -{ - size_t n=checkArrays(pre,point); - checkEqual(n,checkArray(post)); - checkEqual(n,checkArray(straight)); - mem::vector<solvedKnot3> nodes(n); - for(size_t i=0; i < n; ++i) { - nodes[i].pre=read<triple>(pre,i); - nodes[i].point=read<triple>(point,i); - nodes[i].post=read<triple>(post,i); - nodes[i].straight=read<bool>(straight,i); - } - - return path3(nodes,(Int) n,cyclic); -} - -path3 :nullPath3() -{ - return nullpath3; -} - -bool ==(path3 a, path3 b) -{ - return a == b; -} - -bool !=(path3 a, path3 b) -{ - return !(a == b); -} - -triple point(path3 p, Int t) -{ - return p.point((Int) t); -} - -triple point(path3 p, real t) -{ - return p.point(t); -} - -triple precontrol(path3 p, Int t) -{ - return p.precontrol((Int) t); -} - -triple precontrol(path3 p, real t) -{ - return p.precontrol(t); -} - -triple postcontrol(path3 p, Int t) -{ - return p.postcontrol((Int) t); -} - -triple postcontrol(path3 p, real t) -{ - return p.postcontrol(t); -} - -triple dir(path3 p, Int t, Int sign=0, bool normalize=true) -{ - return p.dir(t,sign,normalize); -} - -triple dir(path3 p, real t, bool normalize=true) -{ - return p.dir(t,normalize); -} - -triple accel(path3 p, Int t, Int sign=0) -{ - return p.accel(t,sign); -} - -triple accel(path3 p, real t) -{ - return p.accel(t); -} - -real radius(path3 p, real t) -{ - triple v=p.dir(t,false); - triple a=p.accel(t); - real d=dot(a,v); - real v2=v.abs2(); - real a2=a.abs2(); - real denom=v2*a2-d*d; - real r=v2*sqrt(v2); - return denom > 0 ? r/sqrt(denom) : 0.0; -} - -real radius(triple z0, triple c0, triple c1, triple z1, real t) -{ - triple v=(3.0*(z1-z0)+9.0*(c0-c1))*t*t+(6.0*(z0+c1)-12.0*c0)*t+3.0*(c0-z0); - triple a=6.0*(z1-z0+3.0*(c0-c1))*t+6.0*(z0+c1)-12.0*c0; - real d=dot(a,v); - real v2=v.abs2(); - real a2=a.abs2(); - real denom=v2*a2-d*d; - real r=v2*sqrt(v2); - return denom > 0 ? r/sqrt(denom) : 0.0; -} - -path3 reverse(path3 p) -{ - return p.reverse(); -} - -path3 subpath(path3 p, Int a, Int b) -{ - return p.subpath((Int) a, (Int) b); -} - -path3 subpath(path3 p, real a, real b) -{ - return p.subpath(a,b); -} - -Int length(path3 p) -{ - return p.length(); -} - -bool cyclic(path3 p) -{ - return p.cyclic(); -} - -bool straight(path3 p, Int t) -{ - return p.straight(t); -} - -// Return the component of vector v perpendicular to a unit vector u. -triple perp(triple v, triple u) -{ - return perp(v,u); -} - -// Return the maximum perpendicular deviation of segment i of path3 g -// from a straight line. -real straightness(path3 p, Int t) -{ - if(p.straight(t)) return 0; - triple z0=p.point(t); - triple u=unit(p.point(t+1)-z0); - return ::max(length(perp(p.postcontrol(t)-z0,u)), - length(perp(p.precontrol(t+1)-z0,u))); -} - -// Return the maximum perpendicular deviation of z0..controls c0 and c1..z1 -// from a straight line. -real straightness(triple z0, triple c0, triple c1, triple z1) -{ - triple u=unit(z1-z0); - return ::max(length(perp(c0-z0,u)),length(perp(c1-z0,u))); -} - -bool piecewisestraight(path3 p) -{ - return p.piecewisestraight(); -} - -real arclength(path3 p) -{ - return p.arclength(); -} - -real arctime(path3 p, real dval) -{ - return p.arctime(dval); -} - -realarray* intersect(path3 p, path3 q, real fuzz=-1) -{ - bool exact=fuzz <= 0.0; - if(fuzz < 0) - fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())), - ::max(length(q.max()),length(q.min()))); - - std::vector<real> S,T; - real s,t; - if(intersections(s,t,S,T,p,q,fuzz,true,exact)) { - array *V=new array(2); - (*V)[0]=s; - (*V)[1]=t; - return V; - } else - return new array(0); -} - -realarray2* intersections(path3 p, path3 q, real fuzz=-1) -{ - bool exact=fuzz <= 0.0; - if(fuzz < 0) - fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())), - ::max(length(q.max()),length(q.min()))); - bool single=!exact; - - real s,t; - std::vector<real> S,T; - bool found=intersections(s,t,S,T,p,q,fuzz,single,exact); - if(!found) return new array(0); - array *V; - if(single) { - V=new array(1); - array *Vi=new array(2); - (*V)[0]=Vi; - (*Vi)[0]=s; - (*Vi)[1]=t; - } else { - size_t n=S.size(); - V=new array(n); - for(size_t i=0; i < n; ++i) { - array *Vi=new array(2); - (*V)[i]=Vi; - (*Vi)[0]=S[i]; - (*Vi)[1]=T[i]; - } - } - stable_sort(V->begin(),V->end(),run::compare2<real>()); - return V; -} - -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<real> 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() @@ -3181,498 +975,7 @@ guide *reverse(guide *g) } 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 PRCshininess, real granularity, triple normal, - bool lighton, penarray *colors) -{ - f->append(new drawSurface(*g,straight,*p,opacity,shininess,PRCshininess, - granularity,normal,lighton,*colors)); -} - -triple min3(picture *f) -{ - return f->bounds3().Min(); -} - -triple max3(picture *f) -{ - return f->bounds3().Max(); -} - -pair minratio(picture *f) -{ - return f->ratio(::min); -} - -pair maxratio(picture *f) -{ - return f->ratio(::max); -} - -triple minbezier(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 maxbezier(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 minratio(triplearray2 *P, pair b) -{ - triple *A=copyTripleArray2C(P,true,4); - real fuzz=sqrtFuzz*norm(A,16); - b=pair(bound(A,::min,xratio,b.getx(),fuzz), - bound(A,::min,yratio,b.gety(),fuzz)); - delete[] A; - return b; -} - -pair maxratio(triplearray2 *P, pair b) -{ - triple *A=copyTripleArray2C(P,true,4); - real fuzz=sqrtFuzz*norm(A,16); - b=pair(bound(A,::max,xratio,b.getx(),fuzz), - bound(A,::max,yratio,b.gety(),fuzz)); - delete[] A; - return b; -} - -pair minratio(path3 g) -{ - return g.ratio(::min); -} - -pair maxratio(path3 g) -{ - return g.ratio(::max); -} - -real norm(realarray *a) -{ - size_t n=checkArray(a); - real M=0.0; - for(size_t i=0; i < n; ++i) { - real x=fabs(vm::read<real>(a,i)); - if(x > M) M=x; - } - return M; -} - -real norm(realarray2 *a) -{ - size_t n=checkArray(a); - real M=0.0; - for(size_t i=0; i < n; ++i) { - vm::array *ai=vm::read<vm::array*>(a,i); - size_t m=checkArray(ai); - for(size_t j=0; j < m; ++j) { - real a=fabs(vm::read<real>(ai,j)); - if(a > M) M=a; - } - } - return M; -} - -real norm(triplearray2 *a) -{ - size_t n=checkArray(a); - real M=0.0; - for(size_t i=0; i < n; ++i) { - vm::array *ai=vm::read<vm::array*>(a,i); - size_t m=checkArray(ai); - for(size_t j=0; j < m; ++j) { - real a=vm::read<triple>(ai,j).abs2(); - if(a > M) M=a; - } - } - return sqrt(M); -} - -real change2(triplearray2 *a) -{ - size_t n=checkArray(a); - if(n == 0) return 0.0; - - vm::array *a0=vm::read<vm::array*>(a,0); - size_t m=checkArray(a0); - if(m == 0) return 0.0; - triple a00=vm::read<triple>(a0,0); - real M=0.0; - - for(size_t i=0; i < n; ++i) { - vm::array *ai=vm::read<vm::array*>(a,i); - size_t m=checkArray(ai); - for(size_t j=0; j < m; ++j) { - real a=(vm::read<triple>(ai,j)-a00).abs2(); - if(a > M) M=a; - } - } - return M; -} - -bool is3D(picture *f) -{ - return f->have3D(); -} - -pair bezier(pair a, pair b, pair c, pair d, real t) -{ - real onemt=1-t; - real onemt2=onemt*onemt; - return onemt2*onemt*a+t*(3.0*(onemt2*b+t*onemt*c)+t*t*d); -} - -pair bezierP(pair a, pair b, pair c, pair d, real t) -{ - return 3.0*(t*t*(d-a+3.0*(b-c))+t*(2.0*(a+c)-4.0*b)+b-a); -} - -pair bezierPP(pair a, pair b, pair c, pair d, real t) -{ - return 6.0*(t*(d-a+3.0*(b-c))+a+c-2.0*b); -} - -pair bezierPPP(pair a, pair b, pair c, pair d) -{ - return 6.0*(d-a+3.0*(b-c)); -} - -triple bezier(triple a, triple b, triple c, triple d, real t) -{ - real onemt=1-t; - real onemt2=onemt*onemt; - return onemt2*onemt*a+t*(3.0*(onemt2*b+t*onemt*c)+t*t*d); -} - -triple bezierP(triple a, triple b, triple c, triple d, real t) -{ - return 3.0*(t*t*(d-a+3.0*(b-c))+t*(2.0*(a+c)-4.0*b)+b-a); -} - -triple bezierPP(triple a, triple b, triple c, triple d, real t) -{ - return 6.0*(t*(d-a+3.0*(b-c))+a+c-2.0*b); -} - -triple bezierPPP(triple a, triple b, triple c, triple d) -{ - return 6.0*(d-a+3.0*(b-c)); -} - -// String operations - -string :emptyString() -{ - return emptystring; -} - -Int length(string *s) -{ - return (Int) s->length(); -} - -Int find(string *s, string t, Int pos=0) -{ - size_t n=s->find(t,pos); - return n == string::npos ? (Int) -1 : (Int) n; -} - -Int rfind(string *s, string t, Int pos=-1) -{ - size_t n=s->rfind(t,pos); - return n == string::npos ? (Int) -1 : (Int) n; -} - -string reverse(string s) -{ - reverse(s.begin(),s.end()); - return s; -} - -string insert(string s, Int pos, string t) -{ - if ((size_t) pos < s.length()) - return s.insert(pos,t); - return s; -} - -string substr(string* s, Int pos, Int n=-1) -{ - if ((size_t) pos < s->length()) - return s->substr(pos,n); - return emptystring; -} - -string erase(string s, Int pos, Int n) -{ - if ((size_t) pos < s.length()) - return s.erase(pos,n); - return s; -} - -string downcase(string s) -{ - std::transform(s.begin(),s.end(),s.begin(),tolower); - return s; -} - -string upcase(string s) -{ - std::transform(s.begin(),s.end(),s.begin(),toupper); - return s; -} - -// returns a string constructed by translating all occurrences of the string -// from in an array of string pairs {from,to} to the string to in string s. -string replace(string *S, stringarray2 *translate) -{ - size_t size=checkArray(translate); - for(size_t i=0; i < size; i++) { - array *a=read<array*>(translate,i); - checkArray(a); - } - const char *p=S->c_str(); - ostringstream buf; - while(*p) { - for(size_t i=0; i < size;) { - array *a=read<array*>(translate,i); - string* from=read<string*>(a,0); - size_t len=from->length(); - if(strncmp(p,from->c_str(),len) != 0) {i++; continue;} - buf << read<string>(a,1); - p += len; - if(*p == 0) return buf.str(); - i=0; - } - buf << *(p++); - } - return buf.str(); -} - -string format(string *format, Int x) -{ - const char *f=format->c_str(); - - checkformat(f,true); - - Int size=snprintf(NULL,0,f,x)+1; - if(size < 1) size=255; // Workaround for non-C99 compliant systems. - char *buf=new char[size]; - snprintf(buf,size,f,x); - string s=string(buf); - delete[] buf; - return s; -} - -string format(string *format, real x, string locale=emptystring) -{ - bool tex=getSetting<string>("tex") != "none"; - bool texify=false; - 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; - char prev=0; - while(*p != 0) { - char curr=*p; - if(tex && curr == '$' && prev != '\\') texify=true; - prev=curr; - if(curr == '%') { - 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(texify && (*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() { @@ -3689,185 +992,6 @@ realarray *_cputime() } -// 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) @@ -4026,2129 +1150,4 @@ 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 dir(explicit triple z) -{ - return unit(z); -} - -triple expi(real polar, real azimuth) -{ - return expi(polar,azimuth); -} - -triple dir(real colatitude, real longitude) -{ - return expi(radians(colatitude),radians(longitude)); -} - - -// System routines - -void atupdate(callable *f) -{ - processData().atUpdateFunction=f; -} - -callable *atupdate() -{ - return processData().atUpdateFunction; -} - -void atexit(callable *f) -{ - processData().atExitFunction=f; -} - -callable *atexit() -{ - return processData().atExitFunction; -} - -void atbreakpoint(callableBp *f) -{ - processData().atBreakpointFunction=f; -} - -void breakpoint(runnable *s=NULL) -{ - breakpoint(Stack,s); -} - -string locatefile(string file) -{ - return locateFile(file); -} - -void stop(string file, Int line, runnable *s=NULL) -{ - file=locateFile(file); - clear(file,line); - cout << "setting breakpoint at " << file << ": " << line << endl; - bplist.push_back(bpinfo(file,line,s)); -} - -void breakpoints() -{ - for(mem::list<bpinfo>::iterator p=bplist.begin(); p != bplist.end(); ++p) - cout << p->f.name() << ": " << p->f.line() << endl; -} - -void clear(string file, Int line) -{ - file=locateFile(file); - clear(file,line,true); -} - -void clear() -{ - bplist.clear(); -} - -void warn(string s) -{ - Warn(s); -} - -void nowarn(string s) -{ - noWarn(s); -} - -void warning(string s, string t, bool position=false) -{ - if(settings::warn(s)) { - em.warning(position ? getPos() : nullPos,s); - em << t; - } -} - -// 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) -{ - string name=convertname(file,format); - mem::vector<string> cmd; - cmd.push_back(getSetting<string>("convert")); - push_split(cmd,args); - cmd.push_back(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()) { - mem::vector<string> cmd; - cmd.push_back(getSetting<string>("animate")); - push_split(cmd,args); - cmd.push_back(name); - return System(cmd,0,false,"animate","your animated GIF viewer"); - } -#endif - return 0; -} - -// Delete file named s. -Int delete(string *s) -{ - checkLocal(*s); - Int rc=unlink(s->c_str()); - if(rc == 0 && verbose > 0) - cout << "Deleted " << *s << endl; - return rc; -} - -// Rename file "from" to file "to". -Int rename(string *from, string *to) -{ - checkLocal(*from); - checkLocal(*to); - Int rc=rename(from->c_str(),to->c_str()); - if(rc == 0 && verbose > 0) - cout << "Renamed " << *from << " to " << *to << endl; - return rc; -} - - -// Array operations - -// Create an empty array. -array* :emptyArray() -{ - return new array(0); -} - -// Create a new array (technically a vector). -// This array will be multidimensional. First the number of dimensions -// is popped off the stack, followed by each dimension in reverse order. -// The array itself is technically a one dimensional array of one -// dimension arrays and so on. -array* :newDeepArray(Int depth) -{ - assert(depth > 0); - - Int *dims = new Int[depth]; - - for (Int index = depth-1; index >= 0; index--) { - Int i=pop<Int>(Stack); - if(i < 0) error("cannot create a negative length array"); - dims[index]=i; - } - - array *a=deepArray(depth, dims); - delete[] dims; - return a; -} - -// Creates an array with elements already specified. First, the number -// of elements is popped off the stack, followed by each element in -// reverse order. -array* :newInitializedArray(Int n) -{ - assert(n >= 0); - - array *a = new array(n); - - for (Int index = n-1; index >= 0; index--) - (*a)[index] = pop(Stack); - - return a; -} - -// Similar to newInitializedArray, but after the n elements, append another -// array to it. -array* :newAppendedArray(array* tail, Int n) -{ - assert(n >= 0); - - array *a = new array(n); - - for (Int index = n-1; index >= 0; index--) - (*a)[index] = pop(Stack); - - copy(tail->begin(), tail->end(), back_inserter(*a)); - - return a; -} - -// The function T[] array(int n, T value, int depth=0) produces a array of n -// copies of x, where each copy is copied up to depth. -array* :newDuplicateArray(Int n, item value, Int depth=Int_MAX) -{ - if(n < 0) error("cannot create a negative length array"); - if(depth < 0) error("cannot copy to a negative depth"); - - return new array(n, value, depth); -} - -// Read an element from an array. Checks for initialization & bounds. -item :arrayRead(array *a, Int n) -{ - item& i=arrayRead(a,n); - if (i.empty()) { - ostringstream buf; - buf << "read uninitialized value from array at index " << n; - error(buf); - } - return i; -} - -// Slice a substring from an array. -item :arraySliceRead(array *a, Int left, Int right) -{ - checkArray(a); - return a->slice(left, right); -} - -// Slice a substring from an array. This implements the cases a[i:] and a[:] -// where the endpoint is not given, and assumed to be the length of the array. -item :arraySliceReadToEnd(array *a, Int left) -{ - size_t len=checkArray(a); - return a->slice(left, (Int)len); -} - -// Read an element from an array of arrays. Check bounds and initialize -// as necessary. -item :arrayArrayRead(array *a, Int n) -{ - item& i=arrayRead(a,n); - if (i.empty()) i=new array(0); - return i; -} - -// Write an element to an array. Increase size if necessary. -item :arrayWrite(item value, array *a, Int n) -{ - size_t len=checkArray(a); - bool cyclic=a->cyclic(); - if(cyclic && len > 0) n=imod(n,len); - else { - if(cyclic) outOfBounds("writing cyclic",len,n); - if(n < 0) outOfBounds("writing",len,n); - if(len <= (size_t) n) - a->resize(n+1); - } - (*a)[n] = value; - return value; -} - -array * :arraySliceWrite(array *src, array *dest, Int left, Int right) -{ - checkArray(src); - checkArray(dest); - dest->setSlice(left, right, src); - return src; -} - -array * :arraySliceWriteToEnd(array *src, array *dest, Int left) -{ - checkArray(src); - size_t len=checkArray(dest); - dest->setSlice(left, (Int) len, src); - return src; -} - -// Returns the length of an array. -Int :arrayLength(array *a) -{ - return (Int) checkArray(a); -} - -// Returns an array of integers representing the keys of the array. -array * :arrayKeys(array *a) -{ - size_t size=checkArray(a); - - array *keys=new array(); - for (size_t i=0; i<size; ++i) { - item& cell = (*a)[i]; - if (!cell.empty()) - keys->push((Int)i); - } - - return keys; -} - -// Return the cyclic flag for an array. -bool :arrayCyclicFlag(array *a) -{ - checkArray(a); - return a->cyclic(); -} - -bool :arraySetCyclicFlag(bool b, array *a) -{ - checkArray(a); - a->cyclic(b); - return b; -} - -// Check to see if an array element is initialized. -bool :arrayInitializedHelper(Int n, array *a) -{ - size_t len=checkArray(a); - bool cyclic=a->cyclic(); - if(cyclic && len > 0) n=imod(n,len); - else if(n < 0 || n >= (Int) len) return false; - item&i=(*a)[(unsigned) n]; - return !i.empty(); -} - -// Returns the initialize method for an array. -callable* :arrayInitialized(array *a) -{ - return new thunk(new bfunc(arrayInitializedHelper),a); -} - -// The helper function for the cyclic method that sets the cyclic flag. -void :arrayCyclicHelper(bool b, array *a) -{ - checkArray(a); - a->cyclic(b); -} - -// Set the cyclic flag for an array. -callable* :arrayCyclic(array *a) -{ - return new thunk(new bfunc(arrayCyclicHelper),a); -} - -// The helper function for the push method that does the actual operation. -item :arrayPushHelper(item x, array *a) -{ - checkArray(a); - a->push(x); - return x; -} - -// Returns the push method for an array. -callable* :arrayPush(array *a) -{ - return new thunk(new bfunc(arrayPushHelper),a); -} - -// The helper function for the append method that appends b to a. -void :arrayAppendHelper(array *b, array *a) -{ - checkArray(a); - size_t size=checkArray(b); - for(size_t i=0; i < size; i++) - a->push((*b)[i]); -} - -// Returns the append method for an array. -callable* :arrayAppend(array *a) -{ - return new thunk(new bfunc(arrayAppendHelper),a); -} - -// The helper function for the pop method. -item :arrayPopHelper(array *a) -{ - size_t asize=checkArray(a); - if(asize == 0) - error("cannot pop element from empty array"); - return a->pop(); -} - -// Returns the pop method for an array. -callable* :arrayPop(array *a) -{ - return new thunk(new bfunc(arrayPopHelper),a); -} - -// The helper function for the insert method. -item :arrayInsertHelper(Int i, array *x, array *a) -{ - size_t asize=checkArray(a); - checkArray(x); - if(a->cyclic() && asize > 0) i=imod(i,asize); - if(i < 0 || i > (Int) asize) - outOfBounds("inserting",asize,i); - (*a).insert((*a).begin()+i,(*x).begin(),(*x).end()); -} - -// Returns the insert method for an array. -callable* :arrayInsert(array *a) -{ - return new thunk(new bfunc(arrayInsertHelper),a); -} - -// Returns the delete method for an array. -callable* :arrayDelete(array *a) -{ - return new thunk(new bfunc(arrayDeleteHelper),a); -} - -bool :arrayAlias(array *a, array *b) -{ - return a==b; -} - -// Return array formed by indexing array a with elements of integer array b -array* :arrayIntArray(array *a, array *b) -{ - size_t asize=checkArray(a); - size_t bsize=checkArray(b); - array *r=new array(bsize); - bool cyclic=a->cyclic(); - for(size_t i=0; i < bsize; i++) { - Int index=read<Int>(b,i); - if(cyclic && asize > 0) index=imod(index,asize); - else - if(index < 0 || index >= (Int) asize) - outOfBounds("reading",asize,index); - (*r)[i]=(*a)[index]; - } - return r; -} - -// returns the complement of the integer array a in {0,2,...,n-1}, -// so that b[complement(a,b.length)] yields the complement of b[a]. -Intarray* complement(Intarray *a, Int n) -{ - size_t asize=checkArray(a); - array *r=new array(0); - bool *keep=new bool[n]; - for(Int i=0; i < n; ++i) keep[i]=true; - for(size_t i=0; i < asize; ++i) { - Int j=read<Int>(a,i); - if(j >= 0 && j < n) keep[j]=false; - } - for(Int i=0; i < n; i++) - if(keep[i]) r->push(i); - - delete[] keep; - return r; -} - -// Generate the sequence {f(i) : i=0,1,...n-1} given a function f and integer n -Intarray* :arraySequence(callable *f, Int n) -{ - if(n < 0) n=0; - array *a=new array(n); - for(Int i=0; i < n; ++i) { - Stack->push(i); - f->call(Stack); - (*a)[i]=pop(Stack); - } - return a; -} - -// Return the array {0,1,...n-1} -Intarray *sequence(Int n) -{ - if(n < 0) n=0; - array *a=new array(n); - for(Int i=0; i < n; ++i) { - (*a)[i]=i; - } - return a; -} - -// Apply a function to each element of an array -array* :arrayFunction(callable *f, array *a) -{ - size_t size=checkArray(a); - array *b=new array(size); - for(size_t i=0; i < size; ++i) { - Stack->push((*a)[i]); - f->call(Stack); - (*b)[i]=pop(Stack); - } - return b; -} - -array* :arraySort(array *a, callable *f) -{ - array *c=copyArray(a); - compareFunc=f; - FuncStack=Stack; - stable_sort(c->begin(),c->end(),compareFunction); - return c; -} - -bool all(boolarray *a) -{ - size_t size=checkArray(a); - bool c=true; - for(size_t i=0; i < size; i++) - if(!get<bool>((*a)[i])) {c=false; break;} - return c; -} - -boolarray* !(boolarray* a) -{ - size_t size=checkArray(a); - array *c=new array(size); - for(size_t i=0; i < size; i++) - (*c)[i]=!read<bool>(a,i); - return c; -} - -Int sum(boolarray *a) -{ - size_t size=checkArray(a); - Int sum=0; - for(size_t i=0; i < size; i++) - sum += read<bool>(a,i) ? 1 : 0; - return sum; -} - -array* :arrayCopy(array *a) -{ - return copyArray(a); -} - -array* :arrayConcat(array *a) -{ - // a is an array of arrays to be concatenated together. - // The signature is - // T[] concat(... T[][] a); - - size_t numArgs=checkArray(a); - size_t resultSize=0; - for (size_t i=0; i < numArgs; ++i) { - resultSize += checkArray(a->read<array *>(i)); - } - - array *result=new array(resultSize); - - size_t ri=0; - for (size_t i=0; i < numArgs; ++i) { - array *arg=a->read<array *>(i); - size_t size=checkArray(arg); - - for (size_t j=0; j < size; ++j) { - (*result)[ri]=(*arg)[j]; - ++ri; - } - } - - return result; -} - -array* :array2Copy(array *a) -{ - return copyArray2(a); -} - -array* :array3Copy(array *a) -{ - return copyArray3(a); -} - -array* :array2Transpose(array *a) -{ - size_t asize=checkArray(a); - array *c=new array(0); - for(size_t i=0; i < asize; i++) { - size_t ip=i+1; - array *ai=read<array*>(a,i); - size_t aisize=checkArray(ai); - size_t csize=checkArray(c); - if(csize < aisize) { - c->resize(aisize); - for(size_t j=csize; j < aisize; j++) { - (*c)[j]=new array(ip); - } - } - for(size_t j=0; j < aisize; j++) { - array *cj=read<array*>(c,j); - if(checkArray(cj) < ip) cj->resize(ip); - (*cj)[i]=(*ai)[j]; - } - } - return c; -} - -// a is a rectangular 3D array; perm is an Int array indicating the type of -// permutation (021 or 120, etc; original is 012). -// Transpose by sending respective members to the permutated locations: -// return the array obtained by putting a[i][j][k] into position perm{ijk}. -array* :array3Transpose(array *a, array *perm) -{ - const size_t DIM=3; - - if(checkArray(perm) != DIM) { - ostringstream buf; - buf << "permutation array must have length " << DIM; - error(buf); - } - - size_t* size=new size_t[DIM]; - for(size_t i=0; i < DIM; ++i) size[i]=DIM; - - for(size_t i=0; i < DIM; ++i) { - Int p=read<Int>(perm,i); - size_t P=(size_t) p; - if(p < 0 || P >= DIM) { - ostringstream buf; - buf << "permutation index out of range: " << p; - error(buf); - } - size[P]=P; - } - - for(size_t i=0; i < DIM; ++i) - if(size[i] == DIM) error("permutation indices must be distinct"); - - static const char *rectangular= - "3D transpose implemented for rectangular matrices only"; - - size_t isize=size[0]=checkArray(a); - array *a0=read<array*>(a,0); - size[1]=checkArray(a0); - array *a00=read<array*>(a0,0); - size[2]=checkArray(a00); - for(size_t i=0; i < isize; i++) { - array *ai=read<array*>(a,i); - size_t jsize=checkArray(ai); - if(jsize != size[1]) error(rectangular); - for(size_t j=0; j < jsize; j++) { - array *aij=read<array*>(ai,j); - if(checkArray(aij) != size[2]) error(rectangular); - } - } - - size_t perm0=(size_t) read<Int>(perm,0); - size_t perm1=(size_t) read<Int>(perm,1); - size_t perm2=(size_t) read<Int>(perm,2); - - size_t sizep0=size[perm0]; - size_t sizep1=size[perm1]; - size_t sizep2=size[perm2]; - - array *c=new array(sizep0); - for(size_t i=0; i < sizep0; ++i) { - array *ci=new array(sizep1); - (*c)[i]=ci; - for(size_t j=0; j < sizep1; ++j) { - array *cij=new array(sizep2); - (*ci)[j]=cij; - } - } - - size_t* i=new size_t[DIM]; - - for(i[0]=0; i[0] < size[0]; ++i[0]) { - array *a0=read<array*>(a,i[0]); - for(i[1]=0; i[1] < size[1]; ++i[1]) { - array *a1=read<array*>(a0,i[1]); - for(i[2]=0; i[2] < size[2]; ++i[2]) { - array *c0=read<array*>(c,i[perm0]); - array *c1=read<array*>(c0,i[perm1]); - (*c1)[i[perm2]]=read<real>(a1,i[2]); - } - } - } - - delete [] i; - delete [] size; - - return c; -} - -// In a boolean array, find the index of the nth true value or -1 if not found -// If n is negative, search backwards. -Int find(boolarray *a, Int n=1) -{ - size_t size=checkArray(a); - Int j=-1; - if(n > 0) - for(size_t i=0; i < size; i++) - if(read<bool>(a,i)) { - n--; if(n == 0) {j=(Int) i; break;} - } - if(n < 0) - for(size_t i=size; i > 0;) - if(read<bool>(a,--i)) { - n++; if(n == 0) {j=(Int) i; break;} - } - return j; -} - -// construct vector obtained by replacing those elements of b for which the -// corresponding elements of a are false by the corresponding element of c. -array* :arrayConditional(array *a, array *b, array *c) -{ - size_t size=checkArray(a); - array *r=new array(size); - if(b && c) { - checkArrays(a,b); - checkArrays(b,c); - for(size_t i=0; i < size; i++) - (*r)[i]=read<bool>(a,i) ? (*b)[i] : (*c)[i]; - } else { - r->clear(); - if(b) { - checkArrays(a,b); - for(size_t i=0; i < size; i++) - if(read<bool>(a,i)) r->push((*b)[i]); - } else if(c) { - checkArrays(a,c); - for(size_t i=0; i < size; i++) - if(!read<bool>(a,i)) r->push((*c)[i]); - } - } - return r; -} - -// Return an n x n identity matrix. -realarray2 *identity(Int n) -{ - return Identity(n); -} - -// Return the diagonal matrix with diagonal entries given by a. -realarray2* :diagonal(realarray *a) -{ - size_t n=checkArray(a); - array *c=new array(n); - for(size_t i=0; i < n; ++i) { - array *ci=new array(n); - (*c)[i]=ci; - for(size_t j=0; j < i; ++j) - (*ci)[j]=0.0; - (*ci)[i]=read<real>(a,i); - for(size_t j=i+1; j < n; ++j) - (*ci)[j]=0.0; - } - return c; -} - -// Return the inverse of an n x n matrix a using Gauss-Jordan elimination. -realarray2 *inverse(realarray2 *a) -{ - a=copyArray2(a); - size_t n=checkArray(a); - checkSquare(a); - - inverseAllocate(n); - - for(size_t i=0; i < n; i++) - pivot[i]=0; - - size_t col=0, row=0; - // This is the main loop over the columns to be reduced. - for(size_t i=0; i < n; i++) { - real big=0.0; - // This is the outer loop of the search for a pivot element. - for(size_t j=0; j < n; j++) { - array *aj=read<array*>(a,j); - if(pivot[j] != 1) { - for(size_t k=0; k < n; k++) { - if(pivot[k] == 0) { - real temp=fabs(read<real>(aj,k)); - if(temp >= big) { - big=temp; - row=j; - col=k; - } - } else if(pivot[k] > 1) { - inverseDeallocate(); - error(singular); - } - } - } - } - ++(pivot[col]); - - // Interchange rows, if needed, to put the pivot element on the diagonal. - array *acol=read<array*>(a,col); - if(row != col) { - array *arow=read<array*>(a,row); - for(size_t l=0; l < n; l++) { - real temp=read<real>(arow,l); - (*arow)[l]=read<real>(acol,l); - (*acol)[l]=temp; - } - } - - Row[i]=row; - Col[i]=col; - - // Divide the pivot row by the pivot element. - real denom=read<real>(acol,col); - if(denom == 0.0) { - inverseDeallocate(); - error(singular); - } - real pivinv=1.0/denom; - (*acol)[col]=1.0; - for(size_t l=0; l < n; l++) - (*acol)[l]=read<real>(acol,l)*pivinv; - - // Reduce all rows except for the pivoted one. - for(size_t k=0; k < n; k++) { - if(k != col) { - array *ak=read<array*>(a,k); - real akcol=read<real>(ak,col); - (*ak)[col]=0.0; - for(size_t l=0; l < n; l++) - (*ak)[l]=read<real>(ak,l)-read<real>(acol,l)*akcol; - } - } - } - - // Unscramble the inverse matrix in view of the column interchanges. - for(size_t l=n; l > 0;) { - l--; - size_t r=Row[l]; - size_t c=Col[l]; - if(r != c) { - for(size_t k=0; k < n; k++) { - array *ak=read<array*>(a,k); - real temp=read<real>(ak,r); - (*ak)[r]=read<real>(ak,c); - (*ak)[c]=temp; - } - } - } - inverseDeallocate(); - return a; -} - -// Solve the linear equation ax=b by LU decomposition, returning the -// solution x, where a is an n x n matrix and b is an array of length n. -// If no solution exists, return an empty array. -realarray *solve(realarray2 *a, realarray *b, bool warn=true) -{ - size_t n=checkArray(a); - - if(n == 0) return new array(0); - - size_t m=checkArray(b); - if(m != n) error(incommensurate); - - real *A=copyArray2C(a); - size_t *index=new size_t[n]; - - if(LUdecompose(A,n,index,warn) == 0) - return new array(0); - - array *x=new array(n); - - real *B=copyArrayC(b); - - for(size_t i=0; i < n; ++i) { - size_t ip=index[i]; - real sum=B[ip]; - B[ip]=B[i]; - real *Ai=A+i*n; - for(size_t j=0; j < i; ++j) - sum -= Ai[j]*B[j]; - B[i]=sum; - } - - for(size_t i=n; i > 0;) { - --i; - real sum=B[i]; - real *Ai=A+i*n; - for(size_t j=i+1; j < n; ++j) - sum -= Ai[j]*B[j]; - B[i]=sum/Ai[i]; - } - - for(size_t i=0; i < n; ++i) - (*x)[i]=B[i]; - - delete[] index; - delete[] B; - delete[] A; - - return x; -} - -// Solve the linear equation ax=b by LU decomposition, returning the -// solution x, where a is an n x n matrix and b is an n x m matrix. -// If no solution exists, return an empty array. -realarray2 *solve(realarray2 *a, realarray2 *b, bool warn=true) -{ - size_t n=checkArray(a); - - if(n == 0) return new array(0); - - if(checkArray(b) != n) error(incommensurate); - size_t m=checkArray(read<array*>(b,0)); - - real *A=copyArray2C(a); - real *B=copyArray2C(b,false); - - size_t *index=new size_t[n]; - - if(LUdecompose(A,n,index,warn) == 0) - return new array(0); - - array *x=new array(n); - - for(size_t i=0; i < n; ++i) { - real *Ai=A+i*n; - real *Bi=B+i*m; - real *Bip=B+index[i]*m; - for(size_t k=0; k < m; ++k) { - real sum=Bip[k]; - Bip[k]=Bi[k]; - size_t jk=k; - for(size_t j=0; j < i; ++j, jk += m) - sum -= Ai[j]*B[jk]; - Bi[k]=sum; - } - } - - for(size_t i=n; i > 0;) { - --i; - real *Ai=A+i*n; - real *Bi=B+i*m; - for(size_t k=0; k < m; ++k) { - real sum=Bi[k]; - size_t jk=(i+1)*m+k; - for(size_t j=i+1; j < n; ++j, jk += m) - sum -= Ai[j]*B[jk]; - Bi[k]=sum/Ai[i]; - } - } - - for(size_t i=0; i < n; ++i) { - real *Bi=B+i*m; - array *xi=new array(m); - (*x)[i]=xi; - for(size_t j=0; j < m; ++j) - (*xi)[j]=Bi[j]; - } - - delete[] index; - delete[] B; - delete[] A; - - return x; -} - -// Compute the determinant of an n x n matrix. -real determinant(realarray2 *a) -{ - real *A=copyArray2C(a); - size_t n=checkArray(a); - - real det=LUdecompose(A,n,NULL,false); - size_t n1=n+1; - for(size_t i=0; i < n; ++i) - det *= A[i*n1]; - - delete[] A; - - return det; -} - -realarray *Operator *(realarray2 *a, realarray *b) -{ - size_t n=checkArray(a); - size_t m=checkArray(b); - array *c=new array(n); - real *B=copyArrayC(b); - for(size_t i=0; i < n; ++i) { - array *ai=read<array*>(a,i); - if(checkArray(ai) != m) error(incommensurate); - real sum=0.0; - for(size_t j=0; j < m; ++j) - sum += read<real>(ai,j)*B[j]; - (*c)[i]=sum; - } - delete[] B; - return c; -} - -realarray *Operator *(realarray *a, realarray2 *b) -{ - size_t n=checkArray(a); - if(n != checkArray(b)) error(incommensurate); - real *A=copyArrayC(a); - - array **B=new array*[n]; - array *bk=read<array *>(b,0); - B[0]=bk; - size_t m=bk->size(); - for(size_t k=1; k < n; k++) { - array *bk=read<array *>(b,k); - if(bk->size() != m) error(incommensurate); - B[k]=bk; - } - array *c=new array(m); - - for(size_t i=0; i < m; ++i) { - real sum=0.0; - for(size_t k=0; k < n; ++k) - sum += A[k]*read<real>(B[k],i); - (*c)[i]=sum; - } - delete[] B; - delete[] A; - return c; -} - -realarray2 *Operator *(realarray2 *a, realarray2 *b) -{ - size_t n=checkArray(a); - - size_t nb=checkArray(b); - size_t na0=n == 0 ? 0 : checkArray(read<array*>(a,0)); - if(na0 != nb) - error(incommensurate); - - size_t nb0=nb == 0 ? 0 : checkArray(read<array*>(b,0)); - - array *c=new array(n); - - real *A=copyArray2C(a,false); - real *B=copyArray2C(b,false); - - for(size_t i=0; i < n; ++i) { - real *Ai=A+i*nb; - array *ci=new array(nb0); - (*c)[i]=ci; - for(size_t j=0; j < nb0; ++j) { - real sum=0.0; - size_t kj=j; - for(size_t k=0; k < nb; ++k, kj += nb0) - sum += Ai[k]*B[kj]; - (*ci)[j]=sum; - } - } - - delete[] B; - delete[] A; - - return c; -} - -triple Operator *(realarray2 *t, triple v) -{ - return *t*v; -} - -pair project(triple v, realarray2 *t) -{ - size_t n=checkArray(t); - if(n != 4) error(incommensurate); - array *t0=read<array*>(t,0); - array *t1=read<array*>(t,1); - array *t3=read<array*>(t,3); - if(checkArray(t0) != 4 || checkArray(t1) != 4 || checkArray(t3) != 4) - error(incommensurate); - - real x=v.getx(); - real y=v.gety(); - real z=v.getz(); - - real f=read<real>(t3,0)*x+read<real>(t3,1)*y+read<real>(t3,2)*z+ - read<real>(t3,3); - if(f == 0.0) dividebyzero(); - f=1.0/f; - - return pair((read<real>(t0,0)*x+read<real>(t0,1)*y+read<real>(t0,2)*z+ - read<real>(t0,3))*f, - (read<real>(t1,0)*x+read<real>(t1,1)*y+read<real>(t1,2)*z+ - read<real>(t1,3))*f); -} - -// Compute the dot product of vectors a and b. -real dot(realarray *a, realarray *b) -{ - size_t n=checkArrays(a,b); - real sum=0.0; - for(size_t i=0; i < n; ++i) - sum += read<real>(a,i)*read<real>(b,i); - return sum; -} - -// Solve the problem L\inv f, where f is an n vector and L is the n x n matrix -// -// [ b[0] c[0] a[0] ] -// [ a[1] b[1] c[1] ] -// [ a[2] b[2] c[2] ] -// [ ... ] -// [ c[n-1] a[n-1] b[n-1] ] -realarray *tridiagonal(realarray *a, realarray *b, realarray *c, realarray *f) -{ - size_t n=checkArrays(a,b); - checkEqual(n,checkArray(c)); - checkEqual(n,checkArray(f)); - - array *up=new array(n); - array& u=*up; - - if(n == 0) return up; - - // Special case: zero Dirichlet boundary conditions - if(read<real>(a,0) == 0.0 && read<real>(c,n-1) == 0.0) { - real temp=read<real>(b,0); - if(temp == 0.0) dividebyzero(); - temp=1.0/temp; - - real *work=new real[n]; - u[0]=read<real>(f,0)*temp; - work[0]=-read<real>(c,0)*temp; - - for(size_t i=1; i < n; i++) { - real temp=(read<real>(b,i)+read<real>(a,i)*work[i-1]); - if(temp == 0.0) {delete[] work; dividebyzero();} - temp=1.0/temp; - u[i]=(read<real>(f,i)-read<real>(a,i)*read<real>(u,i-1))*temp; - work[i]=-read<real>(c,i)*temp; - } - - for(size_t i=n-1; i >= 1; i--) - u[i-1]=read<real>(u,i-1)+work[i-1]*read<real>(u,i); - - delete[] work; - return up; - } - - real binv=read<real>(b,0); - if(binv == 0.0) dividebyzero(); - binv=1.0/binv; - - if(n == 1) {u[0]=read<real>(f,0)*binv; return up;} - if(n == 2) { - real factor=(read<real>(b,0)*read<real>(b,1)- - read<real>(a,0)*read<real>(c,1)); - if(factor== 0.0) dividebyzero(); - factor=1.0/factor; - real temp=(read<real>(b,0)*read<real>(f,1)- - read<real>(c,1)*read<real>(f,0))*factor; - u[0]=(read<real>(b,1)*read<real>(f,0)- - read<real>(a,0)*read<real>(f,1))*factor; - u[1]=temp; - return up; - } - - real *gamma=new real[n-2]; - real *delta=new real[n-2]; - - gamma[0]=read<real>(c,0)*binv; - delta[0]=read<real>(a,0)*binv; - u[0]=read<real>(f,0)*binv; - real beta=read<real>(c,n-1); - real fn=read<real>(f,n-1)-beta*read<real>(u,0); - real alpha=read<real>(b,n-1)-beta*delta[0]; - - for(size_t i=1; i <= n-3; i++) { - real alphainv=read<real>(b,i)-read<real>(a,i)*gamma[i-1]; - if(alphainv == 0.0) {delete[] gamma; delete[] delta; dividebyzero();} - alphainv=1.0/alphainv; - beta *= -gamma[i-1]; - gamma[i]=read<real>(c,i)*alphainv; - u[i]=(read<real>(f,i)-read<real>(a,i)*read<real>(u,i-1))*alphainv; - fn -= beta*read<real>(u,i); - delta[i]=-read<real>(a,i)*delta[i-1]*alphainv; - alpha -= beta*delta[i]; - } - - real alphainv=read<real>(b,n-2)-read<real>(a,n-2)*gamma[n-3]; - if(alphainv == 0.0) {delete[] gamma; delete[] delta; dividebyzero();} - alphainv=1.0/alphainv; - u[n-2]=(read<real>(f,n-2)-read<real>(a,n-2)*read<real>(u,n-3)) - *alphainv; - beta=read<real>(a,n-1)-beta*gamma[n-3]; - real dnm1=(read<real>(c,n-2)-read<real>(a,n-2)*delta[n-3])*alphainv; - real temp=alpha-beta*dnm1; - if(temp == 0.0) {delete[] gamma; delete[] delta; dividebyzero();} - u[n-1]=temp=(fn-beta*read<real>(u,n-2))/temp; - u[n-2]=read<real>(u,n-2)-dnm1*temp; - - for(size_t i=n-2; i >= 1; i--) - u[i-1]=read<real>(u,i-1)-gamma[i-1]*read<real>(u,i)-delta[i-1]*temp; - - delete[] delta; - delete[] gamma; - - return up; -} - -// Root solve by Newton-Raphson -real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x, - bool verbose=false) -{ - static const real fuzz=1000.0*DBL_EPSILON; - Int i=0; - size_t oldPrec=0; - if(verbose) - oldPrec=cout.precision(DBL_DIG); - - real diff=DBL_MAX; - real lastdiff; - do { - real x0=x; - - Stack->push(x); - fprime->call(Stack); - real dfdx=pop<real>(Stack); - - if(dfdx == 0.0) { - x=DBL_MAX; - break; - } - - Stack->push(x); - f->call(Stack); - real fx=pop<real>(Stack); - - x -= fx/dfdx; - - lastdiff=diff; - - if(verbose) - cout << "Newton-Raphson: " << x << endl; - - diff=fabs(x-x0); - if(++i == iterations) { - x=DBL_MAX; - break; - } - } while (diff != 0.0 && (diff < lastdiff || diff > fuzz*fabs(x))); - if(verbose) - cout.precision(oldPrec); - return x; -} - -// Root solve by Newton-Raphson bisection -// cf. routine rtsafe (Press et al., Numerical Recipes, 1991). -real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x1, - real x2, bool verbose=false) -{ - static const real fuzz=1000.0*DBL_EPSILON; - size_t oldPrec=0; - if(verbose) - oldPrec=cout.precision(DBL_DIG); - - Stack->push(x1); - f->call(Stack); - real f1=pop<real>(Stack); - if(f1 == 0.0) return x1; - - Stack->push(x2); - f->call(Stack); - real f2=pop<real>(Stack); - if(f2 == 0.0) return x2; - - if((f1 > 0.0 && f2 > 0.0) || (f1 < 0.0 && f2 < 0.0)) { - ostringstream buf; - buf << "root not bracketed, f(x1)=" << f1 << ", f(x2)=" << f2 << endl; - error(buf); - } - - real x=0.5*(x1+x2); - real dxold=fabs(x2-x1); - if(f1 > 0.0) { - real temp=x1; - x1=x2; - x2=temp; - } - - if(verbose) - cout << "midpoint: " << x << endl; - - real dx=dxold; - Stack->push(x); - f->call(Stack); - real y=pop<real>(Stack); - - Stack->push(x); - fprime->call(Stack); - real dy=pop<real>(Stack); - - Int j; - for(j=0; j < iterations; j++) { - if(((x-x2)*dy-y)*((x-x1)*dy-y) >= 0.0 || fabs(2.0*y) > fabs(dxold*dy)) { - dxold=dx; - dx=0.5*(x2-x1); - x=x1+dx; - if(verbose) - cout << "bisection: " << x << endl; - if(x1 == x) return x; - } else { - dxold=dx; - dx=y/dy; - real temp=x; - x -= dx; - if(verbose) - cout << "Newton-Raphson: " << x << endl; - if(temp == x) return x; - } - if(fabs(dx) < fuzz*fabs(x)) return x; - - Stack->push(x); - f->call(Stack); - y=pop<real>(Stack); - - Stack->push(x); - fprime->call(Stack); - dy=pop<real>(Stack); - - if(y < 0.0) x1=x; - else x2=x; - } - if(verbose) - cout.precision(oldPrec); - return (j == iterations) ? DBL_MAX : x; -} - -real simpson(callableReal *f, real a, real b, real acc=DBL_EPSILON, - real dxmax=0) -{ - real integral; - if(dxmax == 0) dxmax=b-a; - Func=f; - FuncStack=Stack; - if(!simpson(integral,wrapFunction,a,b,acc,dxmax)) - error("nesting capacity exceeded in simpson"); - return integral; -} - -// Compute the fast Fourier transform of a pair array -pairarray* :pairArrayFFT(pairarray *a, Int sign=1) -{ - unsigned n=(unsigned) checkArray(a); -#ifdef HAVE_LIBFFTW3 - array *c=new array(n); - if(n) { - Complex *f=FFTWComplex(n); - fft1d Forward(n,intcast(sign),f); - - for(size_t i=0; i < n; i++) { - pair z=read<pair>(a,i); - f[i]=Complex(z.getx(),z.gety()); - } - Forward.fft(f); - - for(size_t i=0; i < n; i++) { - Complex z=f[i]; - (*c)[i]=pair(z.real(),z.imag()); - } - FFTWdelete(f); - } -#else - unused(&n); - unused(&sign); - array *c=new array(0); -#endif // HAVE_LIBFFTW3 - return c; -} - -Intarray2 *triangulate(pairarray *z) -{ - size_t nv=checkArray(z); -// Call robust version of Gilles Dumoulin's port of Paul Bourke's -// triangulation code. - - XYZ *pxyz=new XYZ[nv+3]; - ITRIANGLE *V=new ITRIANGLE[4*nv]; - - for(size_t i=0; i < nv; ++i) { - pair w=read<pair>(z,i); - pxyz[i].p[0]=w.getx(); - pxyz[i].p[1]=w.gety(); - pxyz[i].i=(Int) i; - } - - Int ntri; - Triangulate((Int) nv,pxyz,V,ntri,true,false); - - size_t nt=(size_t) ntri; - array *t=new array(nt); - for(size_t i=0; i < nt; ++i) { - array *ti=new array(3); - (*t)[i]=ti; - ITRIANGLE *Vi=V+i; - (*ti)[0]=pxyz[Vi->p1].i; - (*ti)[1]=pxyz[Vi->p2].i; - (*ti)[2]=pxyz[Vi->p3].i; - } - - delete[] V; - delete[] pxyz; - return t; -} - - -// File operations - -bool ==(file *a, file *b) -{ - return a == b; -} - -bool !=(file *a, file *b) -{ - return a != b; -} - -file* :nullFile() -{ - return &camp::nullfile; -} - -file* input(string name, bool check=true, string comment=commentchar) -{ - char c=comment.empty() ? (char) 0 : comment[0]; - file *f=new ifile(name,c,check); - f->open(); - return f; -} - -file* output(string name, bool update=false, string comment=commentchar) -{ - file *f; - if(update) { - char c=comment.empty() ? (char) 0 : comment[0]; - f=new iofile(name,c); - } else f=new ofile(name); - f->open(); - if(update) f->seek(0,false); - return f; -} - -file* xinput(string name, bool check=true) -{ -#ifdef HAVE_RPC_RPC_H - file *f=new ixfile(name,check); - f->open(); - return f; -#else - ostringstream buf; - buf << name << ": XDR read support not enabled"; - error(buf); - unused(&check); // Suppress unused variable warning -#endif -} - -file* xoutput(string name, bool update=false) -{ -#ifdef HAVE_RPC_RPC_H - file *f; - if(update) - f=new ioxfile(name); - else f=new oxfile(name); - f->open(); - if(update) f->seek(0,false); - return f; -#else - ostringstream buf; - buf << name << ": XDR write support not enabled"; - error(buf); - unused(&update); // Suppress unused variable warning -#endif -} - -file* binput(string name, bool check=true) -{ - file *f=new ibfile(name,check); - f->open(); - return f; -} - -file* boutput(string name, bool update=false) -{ - file *f; - if(update) f=new iobfile(name); - else f=new obfile(name); - f->open(); - if(update) f->seek(0,false); - return f; -} - -bool eof(file *f) -{ - return f->eof(); -} - -bool eol(file *f) -{ - return f->eol(); -} - -bool error(file *f) -{ - return f->error(); -} - -void clear(file *f) -{ - f->clear(); -} - -void close(file *f) -{ - f->close(); -} - -Int precision(file *f=NULL, Int digits=0) -{ - if(f == 0) f=&camp::Stdout; - return f->precision(digits); -} - -void flush(file *f) -{ - f->flush(); -} - -string getc(file *f) -{ - char c=0; - if(f->isOpen()) f->read(c); - static char str[1]; - str[0]=c; - return string(str); -} - -Int tell(file *f) -{ - return f->tell(); -} - -void seek(file *f, Int pos) -{ - f->seek(pos,pos >= 0); -} - -void seekeof(file *f) -{ - f->seek(0,false); -} - -string :namePart(file f) -{ - return f.filename(); -} - -string :modePart(file f) -{ - return f.FileMode(); -} - -// Set file dimensions -file* :dimensionSetHelper(Int nx=-1, Int ny=-1, Int nz=-1, file *f) -{ - f->dimension(nx,ny,nz); - return f; -} - -callable* :dimensionSet(file *f) -{ - return new thunk(new bfunc(dimensionSetHelper),f); -} - -array * :dimensionPart(file f) -{ - array *a=new array(3); - (*a)[0]=f.Nx(); - (*a)[1]=f.Ny(); - (*a)[2]=f.Nz(); - return a; -} - -// Set file f to read arrays in line-at-a-time mode -file* :lineSetHelper(bool b=true, file *f) -{ - f->LineMode(b); - return f; -} - -callable* :lineSet(file *f) -{ - return new thunk(new bfunc(lineSetHelper),f); -} - -bool :linePart(file f) -{ - return f.LineMode(); -} - -// Set file to read comma-separated values -file* :csvSetHelper(bool b=true, file *f) -{ - f->CSVMode(b); - return f; -} - -callable* :csvSet(file *f) -{ - return new thunk(new bfunc(csvSetHelper),f); -} - -bool :csvPart(file f) -{ - return f.CSVMode(); -} - -// Set file to read whitespace-separated values -file* :wordSetHelper(bool b=true, file *f) -{ - f->WordMode(b); - return f; -} - -callable* :wordSet(file *f) -{ - return new thunk(new bfunc(wordSetHelper),f); -} - -bool :wordPart(file f) -{ - return f.WordMode(); -} - -// Set file to read/write single precision real XDR values. -file* :singlerealSetHelper(bool b=true, file *f) -{ - f->SingleReal(b); - return f; -} - -callable* :singlerealSet(file *f) -{ - return new thunk(new bfunc(singlerealSetHelper),f); -} - -bool :singlerealPart(file f) -{ - return f.SingleReal(); -} - -// Set file to read/write single precision int XDR values. -file* :singleintSetHelper(bool b=true, file *f) -{ - f->SingleInt(b); - return f; -} - -callable* :singleintSet(file *f) -{ - return new thunk(new bfunc(singleintSetHelper),f); -} - -bool :singleintPart(file f) -{ - return f.SingleInt(); -} - -// Set file to read/write signed int XDR values. -file* :signedintSetHelper(bool b=true, file *f) -{ - f->SignedInt(b); - return f; -} - -callable* :signedintSet(file *f) -{ - return new thunk(new bfunc(signedintSetHelper),f); -} - -bool :signedintPart(file f) -{ - return f.SignedInt(); -} - -// Set file to read an arrayi (i int sizes followed by an i-dimensional array) -file* :readSetHelper(Int i, file *f) -{ - switch(i) { - case 1: - f->dimension(-2); - break; - - case 2: - f->dimension(-2,-2); - break; - - case 3: - f->dimension(-2,-2,-2); - break; - - default: - f->dimension(); - } - - return f; -} - -callable* :readSet(file *f) -{ - return new thunk(new bfunc(readSetHelper),f); -} - -// 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) || getSetting<bool>("interactive"))) - 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 -} |