1// Asymptote mathematics routines
2
3int quadrant(real degrees)
4{
5  return floor(degrees/90) % 4;
6}
7
8// Roots of unity.
9pair unityroot(int n, int k=1)
10{
11  return expi(2pi*k/n);
12}
13
14real csc(real x) {return 1/sin(x);}
15real sec(real x) {return 1/cos(x);}
16real cot(real x) {return tan(pi/2-x);}
17
18real acsc(real x) {return asin(1/x);}
19real asec(real x) {return acos(1/x);}
20real acot(real x) {return pi/2-atan(x);}
21
22real frac(real x) {return x-(int)x;}
23
24pair exp(explicit pair z) {return exp(z.x)*expi(z.y);}
25pair log(explicit pair z) {return log(abs(z))+I*angle(z);}
26
27// Return an Nx by Ny unit square lattice with lower-left corner at (0,0).
28picture grid(int Nx, int Ny, pen p=currentpen)
29{
30  picture pic;
31  for(int i=0; i <= Nx; ++i) draw(pic,(i,0)--(i,Ny),p);
32  for(int j=0; j <= Ny; ++j) draw(pic,(0,j)--(Nx,j),p);
33  return pic;
34}
35
36bool polygon(path p)
37{
38  return cyclic(p) && piecewisestraight(p);
39}
40
41// Return the intersection time of the point on the line through p and q
42// that is closest to z.
43real intersect(pair p, pair q, pair z)
44{
45  pair u=q-p;
46  real denom=dot(u,u);
47  return denom == 0 ? infinity : dot(z-p,u)/denom;
48}
49
50// Return the intersection time of the extension of the line segment PQ
51// with the plane perpendicular to n and passing through Z.
52real intersect(triple P, triple Q, triple n, triple Z)
53{
54  real d=n.x*Z.x+n.y*Z.y+n.z*Z.z;
55  real denom=n.x*(Q.x-P.x)+n.y*(Q.y-P.y)+n.z*(Q.z-P.z);
56  return denom == 0 ? infinity : (d-n.x*P.x-n.y*P.y-n.z*P.z)/denom;
57}
58
59// Return any point on the intersection of the two planes with normals
60// n0 and n1 passing through points P0 and P1, respectively.
61// If the planes are parallel return (infinity,infinity,infinity).
62triple intersectionpoint(triple n0, triple P0, triple n1, triple P1)
63{
64  real Dx=n0.y*n1.z-n1.y*n0.z;
65  real Dy=n0.z*n1.x-n1.z*n0.x;
66  real Dz=n0.x*n1.y-n1.x*n0.y;
67  if(abs(Dx) > abs(Dy) && abs(Dx) > abs(Dz)) {
68    Dx=1/Dx;
69    real d0=n0.y*P0.y+n0.z*P0.z;
70    real d1=n1.y*P1.y+n1.z*P1.z+n1.x*(P1.x-P0.x);
71    real y=(d0*n1.z-d1*n0.z)*Dx;
72    real z=(d1*n0.y-d0*n1.y)*Dx;
73    return (P0.x,y,z);
74  } else if(abs(Dy) > abs(Dz)) {
75    Dy=1/Dy;
76    real d0=n0.z*P0.z+n0.x*P0.x;
77    real d1=n1.z*P1.z+n1.x*P1.x+n1.y*(P1.y-P0.y);
78    real z=(d0*n1.x-d1*n0.x)*Dy;
79    real x=(d1*n0.z-d0*n1.z)*Dy;
80    return (x,P0.y,z);
81  } else {
82    if(Dz == 0) return (infinity,infinity,infinity);
83    Dz=1/Dz;
84    real d0=n0.x*P0.x+n0.y*P0.y;
85    real d1=n1.x*P1.x+n1.y*P1.y+n1.z*(P1.z-P0.z);
86    real x=(d0*n1.y-d1*n0.y)*Dz;
87    real y=(d1*n0.x-d0*n1.x)*Dz;
88    return (x,y,P0.z);
89  }
90}
91
92// Given a real array a, return its partial sums.
93real[] partialsum(real[] a)
94{
95  real[] b=new real[a.length];
96  real sum=0;
97  for(int i=0; i < a.length; ++i) {
98    sum += a[i];
99    b[i]=sum;
100  }
101  return b;
102}
103
104// Given a real array a, return its partial dx-weighted sums.
105real[] partialsum(real[] a, real[] dx)
106{
107  real[] b=new real[a.length];
108  real sum=0;
109  for(int i=0; i < a.length; ++i) {
110    sum += a[i]*dx[i];
111    b[i]=sum;
112  }
113  return b;
114}
115
116// Given an integer array a, return its partial sums.
117int[] partialsum(int[] a)
118{
119  int[] b=new int[a.length];
120  int sum=0;
121  for(int i=0; i < a.length; ++i) {
122    sum += a[i];
123    b[i]=sum;
124  }
125  return b;
126}
127
128// Given an integer array a, return its partial dx-weighted sums.
129int[] partialsum(int[] a, int[] dx)
130{
131  int[] b=new int[a.length];
132  int sum=0;
133  for(int i=0; i < a.length; ++i) {
134    sum += a[i]*dx[i];
135    b[i]=sum;
136  }
137  return b;
138}
139
140// If strict=false, return whether i > j implies a[i] >= a[j]
141// If strict=true, return whether  i > j implies a[i] > a[j]
142bool increasing(real[] a, bool strict=false)
143{
144  real[] ap=copy(a);
145  ap.delete(0);
146  ap.push(0);
147  bool[] b=strict ? (ap > a) : (ap >= a);
148  b[a.length-1]=true;
149  return all(b);
150}
151
152// Return the first and last indices of consecutive true-element segments
153// of bool[] b.
154int[][] segmentlimits(bool[] b)
155{
156  int[][] segment;
157  bool[] n=copy(b);
158  n.delete(0);
159  n.push(!b[b.length-1]);
160  int[] edge=(b != n) ? sequence(1,b.length) : null;
161  edge.insert(0,0);
162  int stop=edge[0];
163  for(int i=1; i < edge.length; ++i) {
164    int start=stop;
165    stop=edge[i];
166    if(b[start])
167      segment.push(new int[] {start,stop-1});
168  }
169  return segment;
170}
171
172// Return the indices of consecutive true-element segments of bool[] b.
173int[][] segment(bool[] b)
174{
175  int[][] S=segmentlimits(b);
176  return sequence(new int[](int i) {
177      return sequence(S[i][0],S[i][1]);
178    },S.length);
179}
180
181// If the sorted array a does not contain x, insert it sequentially,
182// returning the index of x in the resulting array.
183int unique(real[] a, real x) {
184  int i=search(a,x);
185  if(i == -1 || x != a[i]) {
186    ++i;
187    a.insert(i,x);
188  }
189  return i;
190}
191
192int unique(string[] a, string x) {
193  int i=search(a,x);
194  if(i == -1 || x != a[i]) {
195    ++i;
196    a.insert(i,x);
197  }
198  return i;
199}
200
201bool lexorder(pair a, pair b) {
202  return a.x < b.x || (a.x == b.x && a.y < b.y);
203}
204
205bool lexorder(triple a, triple b) {
206  return a.x < b.x || (a.x == b.x && (a.y < b.y || (a.y == b.y && a.z < b.z)));
207}
208
209real[] zero(int n)
210{
211  return sequence(new real(int) {return 0;},n);
212}
213
214real[][] zero(int n, int m)
215{
216  real[][] M=new real[n][];
217  for(int i=0; i < n; ++i)
218    M[i]=sequence(new real(int) {return 0;},m);
219  return M;
220}
221
222bool square(real[][] m)
223{
224  int n=m.length;
225  for(int i=0; i < n; ++i)
226    if(m[i].length != n) return false;
227  return true;
228}
229
230bool rectangular(real[][] m)
231{
232  int n=m.length;
233  if(n > 0) {
234    int m0=m[0].length;
235    for(int i=1; i < n; ++i)
236      if(m[i].length != m0) return false;
237  }
238  return true;
239}
240
241bool rectangular(pair[][] m)
242{
243  int n=m.length;
244  if(n > 0) {
245    int m0=m[0].length;
246    for(int i=1; i < n; ++i)
247      if(m[i].length != m0) return false;
248  }
249  return true;
250}
251
252bool rectangular(triple[][] m)
253{
254  int n=m.length;
255  if(n > 0) {
256    int m0=m[0].length;
257    for(int i=1; i < n; ++i)
258      if(m[i].length != m0) return false;
259  }
260  return true;
261}
262
263// draw the (infinite) line going through P and Q, without altering the
264// size of picture pic.
265void drawline(picture pic=currentpicture, pair P, pair Q, pen p=currentpen)
266{
267  pic.add(new void (frame f, transform t, transform T, pair m, pair M) {
268      // Reduce the bounds by the size of the pen.
269      m -= min(p); M -= max(p);
270
271      // Calculate the points and direction vector in the transformed space.
272      t=t*T;
273      pair z=t*P;
274      pair v=t*Q-z;
275
276      // Handle horizontal and vertical lines.
277      if(v.x == 0) {
278        if(m.x <= z.x && z.x <= M.x)
279          draw(f,(z.x,m.y)--(z.x,M.y),p);
280      } else if(v.y == 0) {
281        if(m.y <= z.y && z.y <= M.y)
282          draw(f,(m.x,z.y)--(M.x,z.y),p);
283      } else {
284        // Calculate the maximum and minimum t values allowed for the
285        // parametric equation z + t*v
286        real mx=(m.x-z.x)/v.x, Mx=(M.x-z.x)/v.x;
287        real my=(m.y-z.y)/v.y, My=(M.y-z.y)/v.y;
288        real tmin=max(v.x > 0 ? mx : Mx, v.y > 0 ? my : My);
289        real tmax=min(v.x > 0 ? Mx : mx, v.y > 0 ? My : my);
290        if(tmin <= tmax)
291          draw(f,z+tmin*v--z+tmax*v,p);
292      }
293    },true);
294}
295
296real interpolate(real[] x, real[] y, real x0, int i)
297{
298  int n=x.length;
299  if(n == 0) abort("Zero data points in interpolate");
300  if(n == 1) return y[0];
301  if(i < 0) {
302    real dx=x[1]-x[0];
303    return y[0]+(y[1]-y[0])/dx*(x0-x[0]);
304  }
305  if(i >= n-1) {
306    real dx=x[n-1]-x[n-2];
307    return y[n-1]+(y[n-1]-y[n-2])/dx*(x0-x[n-1]);
308  }
309
310  real D=x[i+1]-x[i];
311  real B=(x0-x[i])/D;
312  real A=1.0-B;
313  return A*y[i]+B*y[i+1];
314}
315
316// Linearly interpolate data points (x,y) to (x0,y0), where the elements of
317// real[] x are listed in ascending order and return y0. Values outside the
318// available data range are linearly extrapolated using the first derivative
319// at the nearest endpoint.
320real interpolate(real[] x, real[] y, real x0)
321{
322  return interpolate(x,y,x0,search(x,x0));
323}
324
325private string nopoint="point not found";
326
327// Return the nth intersection time of path g with the vertical line through x.
328real time(path g, real x, int n=0)
329{
330  real[] t=times(g,x);
331  if(t.length <= n) abort(nopoint);
332  return t[n];
333}
334
335// Return the nth intersection time of path g with the horizontal line through
336// (0,z.y).
337real time(path g, explicit pair z, int n=0)
338{
339  real[] t=times(g,z);
340  if(t.length <= n) abort(nopoint);
341  return t[n];
342}
343
344// Return the nth y value of g at x.
345real value(path g, real x, int n=0)
346{
347  return point(g,time(g,x,n)).y;
348}
349
350// Return the nth x value of g at y=z.y.
351real value(path g, explicit pair z, int n=0)
352{
353  return point(g,time(g,(0,z.y),n)).x;
354}
355
356// Return the nth slope of g at x.
357real slope(path g, real x, int n=0)
358{
359  pair a=dir(g,time(g,x,n));
360  return a.y/a.x;
361}
362
363// Return the nth slope of g at y=z.y.
364real slope(path g, explicit pair z, int n=0)
365{
366  pair a=dir(g,time(g,(0,z.y),n));
367  return a.y/a.x;
368}
369
370// A quartic complex root solver based on these references:
371// http://planetmath.org/encyclopedia/GaloisTheoreticDerivationOfTheQuarticFormula.html
372// Neumark, S., Solution of Cubic and Quartic Equations, Pergamon Press
373// Oxford (1965).
374pair[] quarticroots(real a, real b, real c, real d, real e)
375{
376  real Fuzz=100000*realEpsilon;
377
378  // Remove roots at numerical infinity.
379  if(abs(a) <= Fuzz*(abs(b)+Fuzz*(abs(c)+Fuzz*(abs(d)+Fuzz*abs(e)))))
380    return cubicroots(b,c,d,e);
381
382  // Detect roots at numerical zero.
383  if(abs(e) <= Fuzz*(abs(d)+Fuzz*(abs(c)+Fuzz*(abs(b)+Fuzz*abs(a)))))
384    return cubicroots(a,b,c,d);
385
386  real ainv=1/a;
387  b *= ainv;
388  c *= ainv;
389  d *= ainv;
390  e *= ainv;
391
392  pair[] roots;
393  real[] T=cubicroots(1,-2c,c^2+b*d-4e,d^2+b^2*e-b*c*d);
394  if(T.length == 0) return roots;
395
396  real t0=T[0];
397  pair[] sum=quadraticroots((1,0),(b,0),(t0,0));
398  pair[] product=quadraticroots((1,0),(t0-c,0),(e,0));
399
400  if(abs(sum[0]*product[0]+sum[1]*product[1]+d) <
401     abs(sum[0]*product[1]+sum[1]*product[0]+d))
402    product=reverse(product);
403
404  for(int i=0; i < 2; ++i)
405    roots.append(quadraticroots((1,0),-sum[i],product[i]));
406
407  return roots;
408}
409
410pair[][] fft(pair[][] a, int sign=1)
411{
412  pair[][] A=new pair[a.length][];
413  int k=0;
414  for(pair[] v : a) {
415    A[k]=fft(v,sign);
416    ++k;
417  }
418  a=transpose(A);
419  k=0;
420  for(pair[] v : a) {
421    A[k]=fft(v,sign);
422    ++k;
423  }
424  return transpose(A);
425}
426
427// Given a matrix A with independent columns, return
428// the unique vector y minimizing |Ay - b|^2 (the L2 norm).
429// If the columns of A are not linearly independent,
430// throw an error (if warn == true) or return an empty array
431// (if warn == false).
432real[] leastsquares(real[][] A, real[] b, bool warn=true)
433{
434  real[] solution=solve(AtA(A),b*A,warn=false);
435  if (solution.length == 0 && warn)
436    abort("Cannot compute least-squares approximation for " +
437          "a matrix with linearly dependent columns.");
438  return solution;
439}
440
441// Namespace
442struct rootfinder_settings {
443  static real roottolerance=1e-4;
444}
445
446real findroot(real f(real), real a, real b,
447              real tolerance=rootfinder_settings.roottolerance,
448              real fa=f(a), real fb=f(b))
449{
450  return _findroot(f,a,b,tolerance,fa,fb);
451}
452