From e0c6872cf40896c7be36b11dcc744620f10adf1d Mon Sep 17 00:00:00 2001
From: Norbert Preining <norbert@preining.info>
Date: Mon, 2 Sep 2019 13:46:59 +0900
Subject: Initial commit

---
 graphics/asymptote/base/math.asy | 451 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 451 insertions(+)
 create mode 100644 graphics/asymptote/base/math.asy

(limited to 'graphics/asymptote/base/math.asy')

diff --git a/graphics/asymptote/base/math.asy b/graphics/asymptote/base/math.asy
new file mode 100644
index 0000000000..9d10e2b779
--- /dev/null
+++ b/graphics/asymptote/base/math.asy
@@ -0,0 +1,451 @@
+// Asymptote mathematics routines
+
+int quadrant(real degrees)
+{
+  return floor(degrees/90) % 4;
+}
+
+// Roots of unity.
+pair unityroot(int n, int k=1)
+{
+  return expi(2pi*k/n);
+}
+
+real csc(real x) {return 1/sin(x);}
+real sec(real x) {return 1/cos(x);}
+real cot(real x) {return tan(pi/2-x);}
+
+real acsc(real x) {return asin(1/x);}
+real asec(real x) {return acos(1/x);}
+real acot(real x) {return pi/2-atan(x);}
+
+real frac(real x) {return x-(int)x;}
+
+pair exp(explicit pair z) {return exp(z.x)*expi(z.y);}
+pair log(explicit pair z) {return log(abs(z))+I*angle(z);}
+
+// Return an Nx by Ny unit square lattice with lower-left corner at (0,0).
+picture grid(int Nx, int Ny, pen p=currentpen)
+{
+  picture pic;
+  for(int i=0; i <= Nx; ++i) draw(pic,(i,0)--(i,Ny),p);
+  for(int j=0; j <= Ny; ++j) draw(pic,(0,j)--(Nx,j),p);
+  return pic; 
+}
+
+bool polygon(path p)
+{
+  return cyclic(p) && piecewisestraight(p);
+}
+
+// Return the intersection time of the point on the line through p and q
+// that is closest to z.
+real intersect(pair p, pair q, pair z)
+{
+  pair u=q-p;
+  real denom=dot(u,u);
+  return denom == 0 ? infinity : dot(z-p,u)/denom;
+}
+
+// Return the intersection time of the extension of the line segment PQ
+// with the plane perpendicular to n and passing through Z.
+real intersect(triple P, triple Q, triple n, triple Z)
+{
+  real d=n.x*Z.x+n.y*Z.y+n.z*Z.z;
+  real denom=n.x*(Q.x-P.x)+n.y*(Q.y-P.y)+n.z*(Q.z-P.z);
+  return denom == 0 ? infinity : (d-n.x*P.x-n.y*P.y-n.z*P.z)/denom;
+}
+                    
+// Return any point on the intersection of the two planes with normals
+// n0 and n1 passing through points P0 and P1, respectively.
+// If the planes are parallel return (infinity,infinity,infinity).
+triple intersectionpoint(triple n0, triple P0, triple n1, triple P1)
+{
+  real Dx=n0.y*n1.z-n1.y*n0.z;
+  real Dy=n0.z*n1.x-n1.z*n0.x;
+  real Dz=n0.x*n1.y-n1.x*n0.y;
+  if(abs(Dx) > abs(Dy) && abs(Dx) > abs(Dz)) {
+    Dx=1/Dx;
+    real d0=n0.y*P0.y+n0.z*P0.z;
+    real d1=n1.y*P1.y+n1.z*P1.z+n1.x*(P1.x-P0.x);
+    real y=(d0*n1.z-d1*n0.z)*Dx;
+    real z=(d1*n0.y-d0*n1.y)*Dx;
+    return (P0.x,y,z);
+  } else if(abs(Dy) > abs(Dz)) {
+    Dy=1/Dy;
+    real d0=n0.z*P0.z+n0.x*P0.x;
+    real d1=n1.z*P1.z+n1.x*P1.x+n1.y*(P1.y-P0.y);
+    real z=(d0*n1.x-d1*n0.x)*Dy;
+    real x=(d1*n0.z-d0*n1.z)*Dy;
+    return (x,P0.y,z);
+  } else {
+    if(Dz == 0) return (infinity,infinity,infinity);
+    Dz=1/Dz;
+    real d0=n0.x*P0.x+n0.y*P0.y;
+    real d1=n1.x*P1.x+n1.y*P1.y+n1.z*(P1.z-P0.z);
+    real x=(d0*n1.y-d1*n0.y)*Dz;
+    real y=(d1*n0.x-d0*n1.x)*Dz;
+    return (x,y,P0.z);
+  }
+}
+
+// Given a real array a, return its partial sums.
+real[] partialsum(real[] a)
+{
+  real[] b=new real[a.length];
+  real sum=0;
+  for(int i=0; i < a.length; ++i) {
+    sum += a[i];
+    b[i]=sum;
+  }
+  return b;
+}
+
+// Given a real array a, return its partial dx-weighted sums.
+real[] partialsum(real[] a, real[] dx)
+{
+  real[] b=new real[a.length];
+  real sum=0;
+  for(int i=0; i < a.length; ++i) {
+    sum += a[i]*dx[i];
+    b[i]=sum;
+  }
+  return b;
+}
+
+// Given an integer array a, return its partial sums.
+int[] partialsum(int[] a)
+{
+  int[] b=new int[a.length];
+  int sum=0;
+  for(int i=0; i < a.length; ++i) {
+    sum += a[i];
+    b[i]=sum;
+  }
+  return b;
+}
+
+// Given an integer array a, return its partial dx-weighted sums.
+int[] partialsum(int[] a, int[] dx)
+{
+  int[] b=new int[a.length];
+  int sum=0;
+  for(int i=0; i < a.length; ++i) {
+    sum += a[i]*dx[i];
+    b[i]=sum;
+  }
+  return b;
+}
+
+// If strict=false, return whether i > j implies a[i] >= a[j]
+// If strict=true, return whether  i > j implies a[i] > a[j]
+bool increasing(real[] a, bool strict=false)
+{
+  real[] ap=copy(a);
+  ap.delete(0);
+  ap.push(0);
+  bool[] b=strict ? (ap > a) : (ap >= a);
+  b[a.length-1]=true;
+  return all(b);
+}
+
+// Return the first and last indices of consecutive true-element segments
+// of bool[] b.
+int[][] segmentlimits(bool[] b)
+{
+  int[][] segment;
+  bool[] n=copy(b);
+  n.delete(0);
+  n.push(!b[b.length-1]);
+  int[] edge=(b != n) ? sequence(1,b.length) : null;
+  edge.insert(0,0);
+  int stop=edge[0];
+  for(int i=1; i < edge.length; ++i) {
+    int start=stop;
+    stop=edge[i];
+    if(b[start])
+      segment.push(new int[] {start,stop-1});
+  }
+  return segment;
+}
+
+// Return the indices of consecutive true-element segments of bool[] b.
+int[][] segment(bool[] b)
+{
+  int[][] S=segmentlimits(b);
+  return sequence(new int[](int i) {
+      return sequence(S[i][0],S[i][1]);
+    },S.length);
+}
+
+// If the sorted array a does not contain x, insert it sequentially,
+// returning the index of x in the resulting array.
+int unique(real[] a, real x) {
+  int i=search(a,x);
+  if(i == -1 || x != a[i]) {
+    ++i;
+    a.insert(i,x);
+  }
+  return i;
+}
+
+int unique(string[] a, string x) {
+  int i=search(a,x);
+  if(i == -1 || x != a[i]) {
+    ++i;
+    a.insert(i,x);
+  }
+  return i;
+}
+
+bool lexorder(pair a, pair b) {
+  return a.x < b.x || (a.x == b.x && a.y < b.y);
+}
+
+bool lexorder(triple a, triple b) {
+  return a.x < b.x || (a.x == b.x && (a.y < b.y || (a.y == b.y && a.z < b.z)));
+}
+
+real[] zero(int n)
+{
+  return sequence(new real(int) {return 0;},n);
+}
+
+real[][] zero(int n, int m)
+{
+  real[][] M=new real[n][];
+  for(int i=0; i < n; ++i)
+    M[i]=sequence(new real(int) {return 0;},m);
+  return M;
+}
+
+bool square(real[][] m)
+{
+  int n=m.length;
+  for(int i=0; i < n; ++i)
+    if(m[i].length != n) return false;
+  return true;
+}
+
+bool rectangular(real[][] m)
+{
+  int n=m.length;
+  if(n > 0) {
+    int m0=m[0].length;
+    for(int i=1; i < n; ++i)
+      if(m[i].length != m0) return false;
+  }
+  return true;
+}
+
+bool rectangular(pair[][] m)
+{
+  int n=m.length;
+  if(n > 0) {
+    int m0=m[0].length;
+    for(int i=1; i < n; ++i)
+      if(m[i].length != m0) return false;
+  }
+  return true;
+}
+
+bool rectangular(triple[][] m)
+{
+  int n=m.length;
+  if(n > 0) {
+    int m0=m[0].length;
+    for(int i=1; i < n; ++i)
+      if(m[i].length != m0) return false;
+  }
+  return true;
+}
+
+// draw the (infinite) line going through P and Q, without altering the
+// size of picture pic.
+void drawline(picture pic=currentpicture, pair P, pair Q, pen p=currentpen)
+{
+  pic.add(new void (frame f, transform t, transform T, pair m, pair M) {
+      // Reduce the bounds by the size of the pen.
+      m -= min(p); M -= max(p);
+
+      // Calculate the points and direction vector in the transformed space.
+      t=t*T;
+      pair z=t*P;
+      pair v=t*Q-z;
+
+      // Handle horizontal and vertical lines.
+      if(v.x == 0) {
+        if(m.x <= z.x && z.x <= M.x)
+          draw(f,(z.x,m.y)--(z.x,M.y),p);
+      } else if(v.y == 0) {
+        if(m.y <= z.y && z.y <= M.y)
+          draw(f,(m.x,z.y)--(M.x,z.y),p);
+      } else {
+        // Calculate the maximum and minimum t values allowed for the
+        // parametric equation z + t*v
+        real mx=(m.x-z.x)/v.x, Mx=(M.x-z.x)/v.x;
+        real my=(m.y-z.y)/v.y, My=(M.y-z.y)/v.y;
+        real tmin=max(v.x > 0 ? mx : Mx, v.y > 0 ? my : My);
+        real tmax=min(v.x > 0 ? Mx : mx, v.y > 0 ? My : my);
+        if(tmin <= tmax)
+          draw(f,z+tmin*v--z+tmax*v,p);
+      }
+    },true);
+}
+
+real interpolate(real[] x, real[] y, real x0, int i) 
+{
+  int n=x.length;
+  if(n == 0) abort("Zero data points in interpolate");
+  if(n == 1) return y[0];
+  if(i < 0) {
+    real dx=x[1]-x[0];
+    return y[0]+(y[1]-y[0])/dx*(x0-x[0]);
+  }
+  if(i >= n-1) {
+    real dx=x[n-1]-x[n-2];
+    return y[n-1]+(y[n-1]-y[n-2])/dx*(x0-x[n-1]);
+  }
+
+  real D=x[i+1]-x[i];
+  real B=(x0-x[i])/D;
+  real A=1.0-B;
+  return A*y[i]+B*y[i+1];
+}
+
+// Linearly interpolate data points (x,y) to (x0,y0), where the elements of
+// real[] x are listed in ascending order and return y0. Values outside the
+// available data range are linearly extrapolated using the first derivative
+// at the nearest endpoint.
+real interpolate(real[] x, real[] y, real x0) 
+{
+  return interpolate(x,y,x0,search(x,x0));
+}
+
+private string nopoint="point not found";
+
+// Return the nth intersection time of path g with the vertical line through x.
+real time(path g, real x, int n=0)
+{
+  real[] t=times(g,x);
+  if(t.length <= n) abort(nopoint);
+  return t[n];
+}
+
+// Return the nth intersection time of path g with the horizontal line through
+// (0,z.y).
+real time(path g, explicit pair z, int n=0)
+{
+  real[] t=times(g,z);
+  if(t.length <= n) abort(nopoint);
+  return t[n];
+}
+
+// Return the nth y value of g at x.
+real value(path g, real x, int n=0)
+{
+  return point(g,time(g,x,n)).y;
+}
+
+// Return the nth x value of g at y=z.y.
+real value(path g, explicit pair z, int n=0)
+{
+  return point(g,time(g,(0,z.y),n)).x;
+}
+
+// Return the nth slope of g at x.
+real slope(path g, real x, int n=0)
+{
+  pair a=dir(g,time(g,x,n));
+  return a.y/a.x;
+}
+
+// Return the nth slope of g at y=z.y.
+real slope(path g, explicit pair z, int n=0)
+{
+  pair a=dir(g,time(g,(0,z.y),n));
+  return a.y/a.x;
+}
+
+// A quartic complex root solver based on these references:
+// http://planetmath.org/encyclopedia/GaloisTheoreticDerivationOfTheQuarticFormula.html
+// Neumark, S., Solution of Cubic and Quartic Equations, Pergamon Press
+// Oxford (1965).
+pair[] quarticroots(real a, real b, real c, real d, real e)
+{
+  real Fuzz=100000*realEpsilon;
+
+  // Remove roots at numerical infinity.
+  if(abs(a) <= Fuzz*(abs(b)+Fuzz*(abs(c)+Fuzz*(abs(d)+Fuzz*abs(e)))))
+    return cubicroots(b,c,d,e);
+  
+  // Detect roots at numerical zero.
+  if(abs(e) <= Fuzz*(abs(d)+Fuzz*(abs(c)+Fuzz*(abs(b)+Fuzz*abs(a)))))
+    return cubicroots(a,b,c,d);
+
+  real ainv=1/a;
+  b *= ainv;
+  c *= ainv;
+  d *= ainv;
+  e *= ainv;
+  
+  pair[] roots;
+  real[] T=cubicroots(1,-2c,c^2+b*d-4e,d^2+b^2*e-b*c*d);
+  if(T.length == 0) return roots;
+
+  real t0=T[0];
+  pair[] sum=quadraticroots((1,0),(b,0),(t0,0));
+  pair[] product=quadraticroots((1,0),(t0-c,0),(e,0));
+
+  if(abs(sum[0]*product[0]+sum[1]*product[1]+d) <
+     abs(sum[0]*product[1]+sum[1]*product[0]+d))
+    product=reverse(product);
+
+  for(int i=0; i < 2; ++i)
+    roots.append(quadraticroots((1,0),-sum[i],product[i]));
+
+  return roots;
+}
+
+pair[][] fft(pair[][] a, int sign=1)
+{
+  pair[][] A=new pair[a.length][];
+  int k=0;
+  for(pair[] v : a) {
+    A[k]=fft(v,sign);
+    ++k;
+  }
+  a=transpose(A);
+  k=0;
+  for(pair[] v : a) {
+    A[k]=fft(v,sign);
+    ++k;
+  }
+  return transpose(A);
+}
+
+// Given a matrix A with independent columns, return
+// the unique vector y minimizing |Ay - b|^2 (the L2 norm).
+// If the columns of A are not linearly independent,
+// throw an error (if warn == true) or return an empty array
+// (if warn == false).
+real[] leastsquares(real[][] A, real[] b, bool warn=true)
+{
+  real[] solution=solve(AtA(A),b*A,warn=false);
+  if (solution.length == 0 && warn)
+    abort("Cannot compute least-squares approximation for " +
+	  "a matrix with linearly dependent columns.");
+  return solution;
+}
+
+// Namespace
+struct rootfinder_settings {
+  static real roottolerance = 1e-4;
+}
+
+real findroot(real f(real), real a, real b,
+              real tolerance=rootfinder_settings.roottolerance,
+              real fa=f(a), real fb=f(b))
+{
+  return _findroot(f,a,b,tolerance,fa,fb);
+}
-- 
cgit v1.2.3