summaryrefslogtreecommitdiff
path: root/graphics/asymptote/runarray.in
diff options
context:
space:
mode:
Diffstat (limited to 'graphics/asymptote/runarray.in')
-rw-r--r--graphics/asymptote/runarray.in268
1 files changed, 134 insertions, 134 deletions
diff --git a/graphics/asymptote/runarray.in b/graphics/asymptote/runarray.in
index 293c628f18..db9e63d420 100644
--- a/graphics/asymptote/runarray.in
+++ b/graphics/asymptote/runarray.in
@@ -28,7 +28,7 @@ callableReal* => realRealFunction()
#ifdef HAVE_LIBFFTW3
#include "fftw++.h"
-static const char *rectangular="matrix must be rectangular";
+ static const char *rectangular="matrix must be rectangular";
#else
static const char *installFFTW=
"Please install fftw3, run ./configure, and recompile";
@@ -70,7 +70,7 @@ void outOfBounds(const char *op, size_t len, Int n)
error(buf);
}
-inline item& arrayRead(array *a, Int n)
+inline item& arrayRead(array *a, Int n)
{
size_t len=checkArray(a);
bool cyclic=a->cyclic();
@@ -83,7 +83,7 @@ inline item& arrayRead(array *a, Int n)
static array* deepArray(Int depth, Int *dims)
{
assert(depth > 0);
-
+
if (depth == 1) {
return new array(dims[0]);
} else {
@@ -161,7 +161,7 @@ array *copyArray(array *a)
{
size_t size=checkArray(a);
array *c=new array(size);
- for(size_t i=0; i < size; i++)
+ for(size_t i=0; i < size; i++)
(*c)[i]=(*a)[i];
return c;
}
@@ -175,7 +175,7 @@ array *copyArray2(array *a)
size_t aisize=checkArray(ai);
array *ci=new array(aisize);
(*c)[i]=ci;
- for(size_t j=0; j < aisize; j++)
+ for(size_t j=0; j < aisize; j++)
(*ci)[j]=(*ai)[j];
}
return c;
@@ -186,12 +186,12 @@ double *copyTripleArray2Components(array *a, size_t &N, GCPlacement placement)
size_t n=checkArray(a);
N=0;
for(size_t i=0; i < n; i++)
- N += checkArray(read<array*>(a,i));
-
+ N += checkArray(read<array*>(a,i));
+
double *A=(placement == NoGC) ? new double [3*N] :
new(placement) double[3*N];
double *p=A;
-
+
for(size_t i=0; i < n; i++) {
array *ai=read<array*>(a,i);
size_t m=checkArray(ai);
@@ -211,12 +211,12 @@ triple *copyTripleArray2C(array *a, size_t &N, GCPlacement placement)
size_t n=checkArray(a);
N=0;
for(size_t i=0; i < n; i++)
- N += checkArray(read<array*>(a,i));
-
+ N += checkArray(read<array*>(a,i));
+
triple *A=(placement == NoGC) ? new triple [N] :
new(placement) triple[N];
triple *p=A;
-
+
for(size_t i=0; i < n; i++) {
array *ai=read<array*>(a,i);
size_t m=checkArray(ai);
@@ -234,20 +234,20 @@ triple operator *(const array& t, const triple& v)
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 ||
+
+ 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+
@@ -260,14 +260,14 @@ template<class T>
array *mult(array *a, array *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)
+ if(na0 != nb)
error(incommensurate);
-
+
size_t nb0=nb == 0 ? 0 : checkArray(read<array*>(b,0));
-
+
array *c=new array(n);
T *A,*B;
@@ -286,10 +286,10 @@ array *mult(array *a, array *b)
(*ci)[j]=sum;
}
}
-
+
delete[] B;
delete[] A;
-
+
return c;
}
@@ -299,7 +299,7 @@ array *AtA(array *a)
{
size_t n=checkArray(a);
size_t m=n == 0 ? 0 : checkArray(read<array*>(a,0));
-
+
array *c=new array(m);
T *A;
@@ -317,12 +317,12 @@ array *AtA(array *a)
(*ci)[j]=sum;
}
}
-
+
delete[] A;
return c;
}
-double norm(double *a, size_t n)
+double norm(double *a, size_t n)
{
if(n == 0) return 0.0;
double M=fabs(a[0]);
@@ -331,7 +331,7 @@ double norm(double *a, size_t n)
return M;
}
-double norm(triple *a, size_t n)
+double norm(triple *a, size_t n)
{
if(n == 0) return 0.0;
double M=a[0].abs2();
@@ -350,7 +350,7 @@ void transpose(double *a, size_t n)
double temp=a[ij];
a[ij]=a[ji];
a[ji]=temp;
- }
+ }
}
}
@@ -377,16 +377,16 @@ void inverse(double *M, size_t n)
real a=M[0], b=M[1], c=M[2];
real d=M[3], e=M[4], f=M[5];
real g=M[6], h=M[7], i=M[8];
-
+
real A=e*i-f*h;
real B=f*g-d*i;
real C=d*h-e*g;
-
+
real det=a*A+b*B+c*C;
if(det == 0.0)
error(singular);
det=1.0/det;
-
+
M[0]=A*det; M[1]=(c*h-b*i)*det; M[2]=(b*f-c*e)*det;
M[3]=B*det; M[4]=(a*i-c*g)*det; M[5]=(c*d-a*f)*det;
M[6]=C*det; M[7]=(b*g-a*h)*det; M[8]=(a*e-b*d)*det;
@@ -394,10 +394,10 @@ void inverse(double *M, size_t n)
}
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++) {
@@ -422,7 +422,7 @@ void inverse(double *M, size_t n)
}
}
++(pivot[col]);
-
+
// Interchange rows, if needed, to put the pivot element on the diagonal.
double *acol=M+n*col;
if(row != col) {
@@ -433,8 +433,8 @@ void inverse(double *M, size_t n)
acol[k]=temp;
}
}
-
- Row[i]=row;
+
+ Row[i]=row;
Col[i]=col;
// Divide the pivot row by the pivot element.
@@ -445,12 +445,12 @@ void inverse(double *M, size_t n)
}
real pivinv=1.0/denom;
acol[col]=1.0;
- for(size_t k=0; k < n; k++)
+ for(size_t k=0; k < n; k++)
acol[k]=acol[k]*pivinv;
-
+
// Reduce all rows except for the pivoted one.
for(size_t k=0; k < n; k++) {
- if(k != col) {
+ if(k != col) {
double *ak=M+n*k;
real akcol=ak[col];
ak[col]=0.0;
@@ -459,7 +459,7 @@ void inverse(double *M, size_t n)
}
}
}
-
+
// Unscramble the inverse matrix in view of the column interchanges.
for(size_t k=n; k > 0;) {
k--;
@@ -550,7 +550,7 @@ Int LUdecompose(double *a, size_t n, size_t* index, bool warn=true)
swap *= -1;
vv[imax]=vv[j];
}
- if(index)
+ if(index)
index[j]=imax;
if(j != n) {
double denom=aj[j];
@@ -576,7 +576,7 @@ void dividebyzero(size_t i)
buf << "Divide by zero";
error(buf);
}
-
+
void integeroverflow(size_t i)
{
ostringstream buf;
@@ -643,7 +643,7 @@ array* :newAppendedArray(array* tail, Int n)
for (Int index = n-1; index >= 0; index--)
(*a)[index] = pop(Stack);
-
+
copy(tail->begin(), tail->end(), back_inserter(*a));
return a;
@@ -840,7 +840,7 @@ callable* :arrayAppend(array *a)
item :arrayPopHelper(array *a)
{
size_t asize=checkArray(a);
- if(asize == 0)
+ if(asize == 0)
error("cannot pop element from empty array");
return a->pop();
}
@@ -857,7 +857,7 @@ 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)
+ if(i < 0 || i > (Int) asize)
outOfBounds("inserting",asize,i);
(*a).insert((*a).begin()+i,(*x).begin(),(*x).end());
}
@@ -911,7 +911,7 @@ Intarray* complement(Intarray *a, Int n)
}
for(Int i=0; i < n; i++)
if(keep[i]) r->push(i);
-
+
delete[] keep;
return r;
}
@@ -972,7 +972,7 @@ Int :arraySearch(array *a, item key, callable *less)
size_t u=size-1;
if(!compareFunction(key,(*a)[u])) return Intcast(u);
size_t l=0;
-
+
while (l < u) {
size_t i=(l+u)/2;
if(compareFunction(key,(*a)[i])) u=i;
@@ -1064,7 +1064,7 @@ array* :array2Transpose(array *a)
// 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}.
+// 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;
@@ -1074,10 +1074,10 @@ array* :array3Transpose(array *a, array *perm)
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;
@@ -1088,13 +1088,13 @@ array* :array3Transpose(array *a, array *perm)
}
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);
@@ -1109,15 +1109,15 @@ array* :array3Transpose(array *a, array *perm)
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);
@@ -1127,9 +1127,9 @@ array* :array3Transpose(array *a, array *perm)
(*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]) {
@@ -1141,9 +1141,9 @@ array* :array3Transpose(array *a, array *perm)
}
}
}
-
- delete[] i;
- delete[] size;
+
+ delete[] i;
+ delete[] size;
return c;
}
@@ -1167,7 +1167,7 @@ Int find(boolarray *a, Int n=1)
return j;
}
-// Find all indices of true values in a boolean array.
+// Find all indices of true values in a boolean array.
Intarray *findall(boolarray *a)
{
size_t size=checkArray(a);
@@ -1230,24 +1230,24 @@ realarray2 *inverse(realarray2 *a)
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,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,b);
-
+
for(size_t i=0; i < n; ++i) {
size_t ip=index[i];
real sum=B[ip];
@@ -1257,7 +1257,7 @@ realarray *solve(realarray2 *a, realarray *b, bool warn=true)
sum -= Ai[j]*B[j];
B[i]=sum;
}
-
+
for(size_t i=n; i > 0;) {
--i;
real sum=B[i];
@@ -1266,14 +1266,14 @@ realarray *solve(realarray2 *a, realarray *b, bool warn=true)
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;
}
@@ -1283,23 +1283,23 @@ realarray *solve(realarray2 *a, realarray *b, bool warn=true)
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,*B;
copyArray2C(A,a);
copyArray2C(B,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;
@@ -1313,7 +1313,7 @@ realarray2 *solve(realarray2 *a, realarray2 *b, bool warn=true)
Bi[k]=sum;
}
}
-
+
for(size_t i=n; i > 0;) {
--i;
real *Ai=A+i*n;
@@ -1326,7 +1326,7 @@ realarray2 *solve(realarray2 *a, realarray2 *b, bool warn=true)
Bi[k]=sum/Ai[i];
}
}
-
+
for(size_t i=0; i < n; ++i) {
real *Bi=B+i*m;
array *xi=new array(m);
@@ -1334,11 +1334,11 @@ realarray2 *solve(realarray2 *a, realarray2 *b, bool warn=true)
for(size_t j=0; j < m; ++j)
(*xi)[j]=Bi[j];
}
-
+
delete[] index;
delete[] B;
delete[] A;
-
+
return x;
}
@@ -1348,14 +1348,14 @@ real determinant(realarray2 *a)
real *A;
copyArray2C(A,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;
}
@@ -1441,16 +1441,16 @@ pair project(triple v, realarray2 *t)
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+
@@ -1458,7 +1458,7 @@ pair project(triple v, realarray2 *t)
}
// Compute the dot product of vectors a and b.
-real dot(realarray *a, realarray *b)
+real dot(realarray *a, realarray *b)
{
size_t n=checkArrays(a,b);
real sum=0.0;
@@ -1468,7 +1468,7 @@ real dot(realarray *a, realarray *b)
}
// Compute the complex dot product of vectors a and b.
-pair dot(pairarray *a, pairarray *b)
+pair dot(pairarray *a, pairarray *b)
{
size_t n=checkArrays(a,b);
pair sum=zero;
@@ -1489,22 +1489,22 @@ 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();}
@@ -1515,15 +1515,15 @@ realarray *tridiagonal(realarray *a, realarray *b, realarray *c, realarray *f)
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)-
@@ -1537,10 +1537,10 @@ realarray *tridiagonal(realarray *a, realarray *b, realarray *c, realarray *f)
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;
@@ -1559,7 +1559,7 @@ realarray *tridiagonal(realarray *a, realarray *b, realarray *c, realarray *f)
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;
@@ -1571,13 +1571,13 @@ realarray *tridiagonal(realarray *a, realarray *b, realarray *c, realarray *f)
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;
}
@@ -1588,18 +1588,18 @@ real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x,
static const real fuzz=1000.0*DBL_EPSILON;
Int i=0;
size_t oldPrec=0;
- if(verbose)
+ 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;
@@ -1608,14 +1608,14 @@ real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x,
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;
@@ -1635,19 +1635,19 @@ real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x1,
{
static const real fuzz=1000.0*DBL_EPSILON;
size_t oldPrec=0;
- if(verbose)
+ 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;
@@ -1661,7 +1661,7 @@ real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x1,
x1=x2;
x2=temp;
}
-
+
if(verbose)
cout << "midpoint: " << x << endl;
@@ -1669,7 +1669,7 @@ real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x1,
Stack->push(x);
f->call(Stack);
real y=pop<real>(Stack);
-
+
Stack->push(x);
fprime->call(Stack);
real dy=pop<real>(Stack);
@@ -1693,11 +1693,11 @@ real newton(Int iterations=100, callableReal *f, callableReal *fprime, real x1,
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);
@@ -1723,10 +1723,10 @@ real _findroot(callableReal *f, real a, real b, real tolerance,
{
if(fa == 0.0) return a;
if(fb == 0.0) return b;
-
+
const char* oppsign="fa and fb must have opposite signs";
int sign;
-
+
if(fa < 0.0) {
if(fb < 0.0) error(oppsign);
sign=1;
@@ -1743,7 +1743,7 @@ real _findroot(callableReal *f, real a, real b, real tolerance,
while(b-a > tolerance) {
t=(a+b)*0.5;
-
+
Stack->push(t);
f->call(Stack);
ft=sign*pop<double>(Stack);
@@ -1757,11 +1757,11 @@ real _findroot(callableReal *f, real a, real b, real tolerance,
real q_A=2.0*(fa-2.0*ft+fb)*factor*factor;
real q_B=(fb-fa)*factor;
quadraticroots Q=quadraticroots(q_A,q_B,ft);
-
+
// If the interpolation somehow failed, continue on to the next binary
// search step. This may or may not be possible, depending on what
// theoretical guarantees are provided by the quadraticroots function.
-
+
real root;
bool found=Q.roots > 0;
if(found) {
@@ -1774,7 +1774,7 @@ real _findroot(callableReal *f, real a, real b, real tolerance,
}
}
}
-
+
if(found) {
if(ft > 0.0) {
b=t;
@@ -1783,7 +1783,7 @@ real _findroot(callableReal *f, real a, real b, real tolerance,
a=t;
fa=ft;
}
-
+
t=root;
// If the interpolated value is close to one edge of
@@ -1792,7 +1792,7 @@ real _findroot(callableReal *f, real a, real b, real tolerance,
real margin=(b-a)*1.0e-3;
if(t-a < margin) t=a+2.0*(t-a);
else if(b-t < margin) t=b-2.0*(b-t);
-
+
Stack->push(t);
f->call(Stack);
ft=sign*pop<double>(Stack);
@@ -1800,7 +1800,7 @@ real _findroot(callableReal *f, real a, real b, real tolerance,
if(ft == 0.0) return t;
}
}
-
+
if(ft > 0.0) {
b=t;
fb=ft;
@@ -1835,13 +1835,13 @@ pairarray* fft(pairarray *a, Int sign=1)
if(n) {
Complex *f=utils::ComplexAlign(n);
fftwpp::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());
@@ -1970,14 +1970,14 @@ Intarray2 *triangulate(pairarray *z)
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);
@@ -1991,7 +1991,7 @@ Intarray2 *triangulate(pairarray *z)
(*ti)[1]=pxyz[Vi->p2].i;
(*ti)[2]=pxyz[Vi->p3].i;
}
-
+
delete[] V;
delete[] pxyz;
return t;
@@ -2042,13 +2042,13 @@ 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);
@@ -2115,25 +2115,25 @@ realarray *_projection()
gl::projection P=gl::camera();
size_t k=0;
(*a)[k++]=P.orthographic ? 1.0 : 0.0;
-
+
triple camera=P.camera;
(*a)[k++]=camera.getx();
(*a)[k++]=camera.gety();
(*a)[k++]=camera.getz();
-
+
triple up=P.up;
(*a)[k++]=up.getx();
(*a)[k++]=up.gety();
(*a)[k++]=up.getz();
-
+
triple target=P.target;
(*a)[k++]=target.getx();
(*a)[k++]=target.gety();
(*a)[k++]=target.getz();
-
+
(*a)[k++]=P.zoom;
(*a)[k++]=P.angle;
-
+
(*a)[k++]=P.viewportshift.getx();
(*a)[k++]=P.viewportshift.gety();
#endif