diff options
author | Karl Berry <karl@freefriends.org> | 2013-04-07 18:19:31 +0000 |
---|---|---|
committer | Karl Berry <karl@freefriends.org> | 2013-04-07 18:19:31 +0000 |
commit | 752012c605d34cd943795527a9738475a6958fcc (patch) | |
tree | 4ee06acdd8333a662c2d6f6ef716235053468f55 /Master/texmf-dist/asymptote/math.asy | |
parent | 9789d09132f18a838e84f041b4b3aff28d3426ec (diff) |
texmf -> texmf-dist: start with unique dirs from texmf
git-svn-id: svn://tug.org/texlive/trunk@29712 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf-dist/asymptote/math.asy')
-rw-r--r-- | Master/texmf-dist/asymptote/math.asy | 417 |
1 files changed, 417 insertions, 0 deletions
diff --git a/Master/texmf-dist/asymptote/math.asy b/Master/texmf-dist/asymptote/math.asy new file mode 100644 index 00000000000..77c0664895f --- /dev/null +++ b/Master/texmf-dist/asymptote/math.asy @@ -0,0 +1,417 @@ +// 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 indices of consecutive true-element segments of bool[] b. +int[][] segment(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=0; i < edge.length-1;) { + int start=stop; + stop=edge[++i]; + if(b[start]) + segment.push(sequence(start,stop-1)); + } + return segment; +} + +// 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; + } + 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; + } + 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); +} |