summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/runtime.in
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/utils/asymptote/runtime.in')
-rw-r--r--Build/source/utils/asymptote/runtime.in5003
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
-}