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