path: root/graphics/asymptote/base/math.asy

// 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 fuzz=-1)
{
  real[] t=times(g,x,fuzz);
  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 fuzz=-1)
{
  real[] t=times(g,z,fuzz);
  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, real fuzz=-1)
{
  return point(g,time(g,x,n,fuzz)).y;
}

// Return the nth x value of g at y=z.y.
real value(path g, explicit pair z, int n=0, real fuzz=-1)
{
  return point(g,time(g,(0,z.y),n,fuzz)).x;
}

// Return the nth slope of g at x.
real slope(path g, real x, int n=0, real fuzz=-1)
{
  pair a=dir(g,time(g,x,n,fuzz));
  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, real fuzz=-1)
{
  pair a=dir(g,time(g,(0,z.y),n,fuzz));
  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);
}