home *** CD-ROM | disk | FTP | other *** search
/ Crawly Crypt Collection 1 / crawlyvol1.bin / apps / math / euler / progs / util.e < prev   
Text File  |  1993-04-16  |  29KB  |  1,101 lines
  1. .. The Utilities.
  2.  
  3. "Loading utilities, Version 3.04, Apr 16, 1993."
  4.  
  5. Epsilon=epsilon();
  6. ResetEpsilon=epsilon();
  7. ResetView=view();
  8. Pi=pi();
  9. E=exp(1);
  10. I=1i;
  11.  
  12. function reset
  13. ## reset() resets espilon() and some other things
  14.     global ResetEpsilon,ResetView;
  15.     setepsilon(ResetEpsilon);
  16.     style(""); hold off;
  17.     color(1);
  18.     shrinkwindow();
  19.     shortformat();
  20.     linewidth(1);
  21.     view(ResetView);
  22.     huecolor(1);
  23.     return 0;
  24. endfunction
  25.  
  26. function wait
  27. ## waits for a key (at most 5 min)
  28.     return wait(300);
  29. endfunction
  30.  
  31. function longformat
  32. ## longformat() sets a long format for numbers
  33.     return format([20 12]);
  34. endfunction
  35.  
  36. function shortformat
  37. ## shortformat() sets a short format for numbers
  38.     return format([14 5]);
  39. endfunction;
  40.  
  41. function linspace (a,b,n)
  42. ## linspace(a,b,n) generates n+1 linear spaced points in [a,b].
  43.     if a~=b; return a; endif;
  44.     r=a:(b-a)/n:b;
  45.     r[n+1]=b;
  46.     return r;
  47. endfunction
  48.  
  49. function equispace (a,b,n)
  50. ## equispace(a,b,n) generates n+1 euqidistribution (acos) spaced values
  51. ## in [a,b].
  52.     m=(1-cos(0:pi()/n:pi()))/2;
  53.     return a+(b-a)*m;
  54. endfunction
  55.  
  56. function length (v)
  57. ## length(v) returns the length of a vector
  58.     return max(size(v));
  59. endfunction
  60.  
  61. function polydif (p)
  62. ## polydif(p) returns the polynomial p'
  63.     n=cols(p);
  64.     if (n==1); return 0; endif;
  65.     return p[2:n]*(1:n-1);
  66. endfunction
  67.  
  68. function writeform (x)
  69.     if isreal(x); return printf("%25.16e",x); endif;
  70.     if iscomplex(x);
  71.         return printf("%25.16e",re(x))|printf("+%25.16ei",im(x));
  72.     endif;
  73.     if isstring(x); return x; endif;
  74.     error("Cannot print this!");
  75. endfunction
  76.  
  77. function write (x,s)
  78.     if argn()==1; s=name(x); endif;
  79.     si=size(x);
  80.     if max(si)==1; s|" = "|writeform(x)|";", return 0; endif;
  81.     s|" = [ .."
  82.     for i=1 to si[1];
  83.         for j=1 to si[2]-1;
  84.             writeform(x[i,j])|",",
  85.         end;
  86.         if i==si[1]; writeform(x[i,si[2]])|"];",
  87.         else; writeform(x[i,si[2]])|";",
  88.         endif;
  89.     end;
  90.     return 0
  91. endfunction
  92.  
  93. .. ### plot things ###
  94.  
  95. function scalematrix(A)
  96. ## Scales a matrix A, so that its value are in the range from 0 to 1.
  97.     e=extrema(A)'; mi=min(e[1]); ma=max(e[3]);
  98.     if ma~=mi; ma=mi+1; endif;
  99.     return (A-mi)/(ma-mi);
  100. endfunction    
  101.  
  102. function select
  103. ## Returns coordinates {x,y} of mouse clicks, until the user clicked
  104. ## above the plot window.
  105.     p=plot();
  106.     x=[];
  107.     repeat
  108.         m=mouse();
  109.         if m[2]>p[4]; break; endif;
  110.         h=holding(1); mark(m[1],m[2]); holding(h);
  111.         x=x_m;
  112.     end;
  113.     x=x'; return {x[1],x[2]}
  114. endfunction
  115.  
  116. function title (text)
  117. ## title(text) plots a title to the grafik window.
  118.     ctext(text,[512 0]);
  119.     return text;
  120. endfunction
  121.  
  122. function textheight
  123. ## textheight() returns the height of a letter.
  124.     h=textsize();
  125.     return h[2];
  126. endfunction
  127.  
  128. function textwidth
  129. ## textwidth() returns the width of a letter.
  130.     h=textsize();
  131.     return h[1];
  132. endfunction
  133.  
  134. function fullwindow
  135. ## fullwindow() takes the full size (besides a title) for the
  136. ## plots.
  137.     h=textsize();
  138.     w=window();
  139.     window([12,textheight()*1.5,1011,1011]);
  140.     return w;
  141. endfunction
  142.  
  143. function shrinkwindow
  144. ## shrinkwindow() shrinks the window to allow labels.
  145.     h=textheight(); b=textwidth();
  146.     w=window();
  147.     window([8*b,1.5*h,1023-2*b,1023-3*h]);
  148.     return w;
  149. endfunction
  150.  
  151. function setplot
  152. ## setplot([xmin xmax ymin ymax]) sets the plot coordinates and holds
  153. ## the plot. Also setplot(xmin,xmax,ymin,ymax).
  154. ## setplot() resets it. 
  155.     if argn()==4; return setplot([arg1,arg2,arg3,arg4]);
  156.     else;
  157.         if argn()==0; scaling(1); 
  158.         else; error("Illegal arguments!");
  159.         endif;
  160.     endif;
  161.     return plot();
  162. endfunction
  163.  
  164. function ticks (a,b)
  165. ## returns the proper ticks to be used for intervals [a,b] and
  166. ## the factor of the ticks.
  167.     tick=10^floor(log(b-a)/log(10)-0.4);
  168.     if b-a>10*tick; tick1=tick*2; else; tick1=tick; endif;
  169.     if (tick>0.001) && (tick<1000); tick=1; endif;
  170.     return {(floor(a/tick1)+1)*tick1:tick1:(ceil(b/tick1)-1)*tick1,tick}
  171. endfunction
  172.  
  173. function xplot (x,y,grid=1,ticks=1)
  174. ## xplot(x,y) or xplot(y) works like plot, but shows axis ticks.
  175. ## xplot() shows only axis ticks and the grid.
  176.     if argn()>0;
  177.         if argn()==1;
  178.             if iscomplex(x); y=im(x); xh=re(x);
  179.             else y=x; xh=1:cols(y);
  180.             endif;
  181.         else; xh=x;
  182.         endif;
  183.         p=plotarea(xh,y);
  184.         if !holding(); clg; frame(); endif;
  185.     else
  186.         p=plot();
  187.     endif;
  188.     lw=linewidth(1);
  189.     {t,f}=ticks(p[1],p[2]);
  190.     xgrid(t,f,grid,ticks);
  191.     {t,f}=ticks(p[3],p[4]);
  192.     ygrid(t,f,grid,ticks);
  193.     linewidth(lw);
  194.     if argn()>0;
  195.         ho=holding(1); plot(xh,y); holding(ho);
  196.     endif;
  197.     return p;
  198. endfunction
  199.  
  200. function xmark (x,y,grid=1,ticks=1)
  201. ## xmark(x,y) or xmark(y) works like plot, but shows axis ticks.
  202. ## xmark() shows only axis ticks and the grid.
  203.     if !holding(); clg; frame(); endif;
  204.     if argn()==1;
  205.         if iscomplex(x); y=im(x); xh=re(x);
  206.         else; y=x; xh=1:cols(y);
  207.         endif;
  208.     else; xh=x;
  209.     endif;
  210.     p=plotarea(xh,y);
  211.     {t,f}=ticks(p[1],p[2]);
  212.     xgrid(t,f,grid,ticks);
  213.     {t,f}=ticks(p[3],p[4]);
  214.     ygrid(t,f,grid,ticks);
  215.     ho=holding(1); p=mark(xh,y); holding(ho);
  216.     return p;
  217. endfunction
  218.  
  219. function setplotm
  220. ## The user may choose the plotting coordinates with the mouse.
  221. ## Returns the plot coordinates.
  222.     h=holding(1);
  223.     k1=mouse(); mark(k1[1],k1[2]);
  224.     k2=mouse(); mark(k2[1],k2[2]);
  225.     kl=min(k1,k2); ku=max(k1,k2);
  226.     c=color(2);
  227.     plot([kl[1],kl[1],ku[1],ku[1],kl[1]],[kl[2],ku[2],ku[2],kl[2],kl[2]]);
  228.     color(c);
  229.     setplot(kl[1],ku[1],kl[2],ku[2]);
  230.     holding(h);
  231.     return plot();
  232. endfunction
  233.  
  234. function fplot (ffunction,a=0,b=0,n=200)
  235. ## fplot("f",a,b,n,...) plots the function f in [a,b].
  236. ## The arguments ... are passed to f.
  237. ## fplot("f") or fplot("f",,,n,...) plots f in the old interval.
  238.     if a~=b; s=plot(); a=s[1]; b=s[2]; endif;
  239.     t=linspace(a,b,n);
  240.     s=ffunction(t,args(4));
  241.     return xplot(t,s)
  242. endfunction
  243.  
  244. function histogram (d,n=10,grid=0,ticks=1,width=4,integer=0)
  245. ## Plots a histogram of the data d with n intervals.
  246. ## d must be a vector.
  247.     mi=min(d); ma=max(d);
  248.     if mi~=ma; ma=mi+1; endif;
  249.     if integer; n=floor(ma-mi)+1; ma=ma+0.9999; endif;
  250.     x=zeros(1,2*n+2); y=zeros(1,2*n+2);
  251.     h=(d-mi)/(ma-mi)*n;
  252.     c=count(h,n+1); c[n]=c[n]+c[n+1]; c=c[1:n];
  253.     y[2:2:2*n]=c; y[3:2:2*n+1]=c;
  254.     xx=linspace(mi,ma,n);
  255.     x[2:2:2*n]=xx[1:n]; x[3:2:2*n+1]=xx[2:n+1];
  256.     x[1]=mi; x[2*n+2]=ma;
  257.     w=linewidth(width); xplot(x,y,grid,ticks); linewidth(w);
  258.     return plot;
  259. endfunction
  260.  
  261. function printscale (x)
  262.     if (abs(x)>10000) || (abs(x)<0.00001);
  263.         return printf("%12.5e",x);
  264.     else
  265.         return printf("%10.5f",x);
  266.     endif;
  267. endfunction
  268.  
  269. function niceform (x)
  270. ## Return a string, containing a nice print of x.
  271.     y=round(x,10);
  272.     return printf("%g",y);
  273. endfunction
  274.  
  275. function xgrid(xx,f=1,grid=1,ticks=1,color=3)
  276. ## xgrid([x0,x1,...]) draws vertical grid lines on the plot window at
  277. ## x0,x1,...
  278. ## xgrid([x0,x1,...],f) additionally writes x0/f to the axis.
  279.     c=plot(); n=cols(xx); s=scaling(0); h=holding(1);
  280.     w=window();
  281.     st=linestyle("."); color(color);
  282.     ht=textheight();
  283.     if argn()<2; ticks=0; endif;
  284.     loop 1 to n;
  285.         x=xx[index()];
  286.         if (x<=c[2])&&(x>=c[1]); 
  287.             if grid; plot([x,x],[c[3],c[4]]); endif;
  288.             if ticks;
  289.                 col=w[1]+(x-c[1])/(c[2]-c[1])*(w[3]-w[1]);
  290.                 ctext(niceform(x/f),[col,w[4]+0.2*ht]);
  291.             endif;
  292.         endif;
  293.     end;
  294.     if ticks && !(f~=1);
  295.         ctext("* "|printscale(f),[(w[1]+w[3])/2,w[4]+1.5*ht]);
  296.     endif;
  297.     linestyle(st); color(1); holding(h); scaling(s);
  298.     return 0;
  299. endfunction
  300.  
  301. function ygrid(yy,f=1,grid=1,ticks=1,color=3)
  302. ## ygrid([x0,x1,...]) draws horizontal grid lines on the plot window at
  303. ## x0,x1,...
  304. ## ygrid([x0,x1,...],f) additionally writes x0/f to the axis.
  305.     c=plot(); n=cols(yy); s=scaling(0); h=holding(1);
  306.     st=linestyle("."); color(color);
  307.     w=window(); wd=textwidth(); ht=textheight();
  308.     if argn()<2; ticks=0; endif;
  309.     loop 1 to n;
  310.         y=yy[index()];
  311.         if (y>=c[3])&&(y<=c[4]);
  312.             if ticks;
  313.                 row=w[4]-(y-c[3])/(c[4]-c[3])*(w[4]-w[2]);
  314.                 rtext(niceform(y/f),[w[1]-wd/2,row-ht/2]);
  315.             endif;
  316.             if grid; plot([c[1],c[2]],[y,y]); endif;
  317.         endif;
  318.     end;
  319.     if ticks && !(f~=1);
  320.         text("* "|printscale(f),[w[1]-6*wd,0]);
  321.     endif;
  322.     linestyle(st); color(1); holding(h); scaling(s);
  323.     return 0;
  324. endfunction
  325.  
  326. function plot (x,y)
  327. ## plot(x,y) plots the values (x(i),y(i)) with the current style.
  328. ## If x is a matrix, y must be a matrix of the same size.
  329. ## The plot is then drawn for all rows of x and y.
  330. ## The plot is scaled automatically, unless hold is on.
  331. ## plot(x,y) and plot() return [x1,x2,y1,y2], where [x1,x2] is the range
  332. ## of the x-values and [y1,y2] of the y-values.
  333. ## plot(x) is the same as plot(1:cols(x),x).
  334.     if argn()==1;
  335.         if iscomplex(x); return plot(re(x),im(x));
  336.         else return plot(1:cols(x),x);
  337.         endif;
  338.     endif;
  339.     error("Illegal argument number!"),
  340. endfunction
  341.  
  342. function mark (x,y)
  343. ## mark(x,y) plots markers at (x(i),y(i)) according the the actual style.
  344. ## Works like plot.
  345.     if argn()==1;
  346.         if iscomplex(x); return mark(re(x),im(x));
  347.         else return mark(1:cols(x),x);
  348.         endif;
  349.     endif;
  350.     error("Illegal argument number!"),
  351. endfunction
  352.  
  353. function cplot (z)
  354. ## cplot(z) plots a grid of complex numbers.
  355.     plot(re(z),im(z)); s=scaling(0); h=holding(1);
  356.     w=z'; plot(re(w),im(w)); holding(h); scaling(s);
  357.     return plot();
  358. endfunction
  359.  
  360. function plotwindow
  361. ## plotwindow() sets the plot window to the screen coordinates.
  362.     w=window();
  363.     setplot(w[1],w[3],w[2],w[4]);
  364.     return plot()
  365. endfunction
  366.  
  367. function getframe (x,y,z)
  368. ## gets a box around all points in (x,y,z).
  369.     ex=extrema(x); ey=extrema(y); ez=extrema(z);
  370.     return [min(ex[:,1]'),max(ex[:,3]'),
  371.         min(ey[:,1]'),max(ey[:,3]'),min(ez[:,1]'),max(ez[:,3]')]
  372. endfunction
  373.  
  374. function framez0 (f)
  375.     wire([f[1],f[2],f[2],f[1],f[1]], ..
  376.         [f[3],f[3],f[4],f[4],f[3]],dup(f[5],5)');
  377.     return 0    
  378. endfunction
  379.  
  380. function framez1 (f)
  381.     wire([f[1],f[2],f[2],f[1],f[1]], ..
  382.         [f[3],f[3],f[4],f[4],f[3]],dup(f[6],5)');
  383.     return 0
  384. endfunction
  385.  
  386. function framexpyp (f)
  387.     wire([f[2],f[2]],[f[4],f[4]],[f[5],f[6]]);
  388.     return 0
  389. endfunction
  390.  
  391. function framexpym (f)
  392.     wire([f[2],f[2]],[f[3],f[3]],[f[5],f[6]]);
  393.     return 0
  394. endfunction
  395.  
  396. function framexmyp (f)
  397.     wire([f[1],f[1]],[f[4],f[4]],[f[5],f[6]]);
  398.     return 0
  399. endfunction
  400.  
  401. function framexmym (f)
  402.     wire([f[1],f[1]],[f[3],f[3]],[f[5],f[6]]);
  403.     return 0
  404. endfunction
  405.  
  406. function frame1 (f)
  407. ## draws the back part of the box (considering view)
  408.     v=view();
  409.     y=sin(v[4])*v[1];
  410.     if y>f[5]; framez0(f); endif;
  411.     if y<f[6]; framez1(f); endif;
  412.     x=sin(v[3])*v[1]; y=-cos(v[3])*v[1];
  413.     if x<=f[2]; framexpyp(f); framexpym(f); endif;
  414.     if x>=f[1]; framexmyp(f); framexmym(f); endif;
  415.     if y<=f[4]; framexmyp(f); framexpyp(f); endif;
  416.     if y>=f[3]; framexmym(f); framexpym(f); endif;
  417.     return 0
  418. endfunction
  419.  
  420. function frame2 (f)
  421. ## draws the front part of the box (considering view).
  422.     v=view();
  423.     y=sin(v[4])*v[1];
  424.     if y<=f[5]; framez0(f); endif;
  425.     if y>=f[6]; framez1(f); endif;
  426.     x=sin(v[3])*v[1]; y=-cos(v[3])*v[1];
  427.     if x>=f[2]; framexpyp(f); framexpym(f); endif;
  428.     if x<=f[1]; framexmyp(f); framexmym(f); endif;
  429.     if y>=f[4]; framexmyp(f); framexpyp(f); endif;
  430.     if y<=f[3]; framexmym(f); framexpym(f); endif;
  431.     return 0
  432. endfunction
  433.  
  434. function scaleframe (x,y,z,f,m)
  435.     s=max([f[2]-f[1],f[4]-f[3],f[6]-f[5]]);
  436.     xm=(f[2]+f[1])/2;
  437.     ym=(f[4]+f[3])/2;
  438.     zm=(f[6]+f[5])/2;
  439.     ff=m/s*2;
  440.     f=[f[1]-xm,f[2]-xm,f[3]-ym,f[4]-ym,f[5]-zm,f[6]-zm]*ff;
  441.     return {(x-xm)*ff,(y-ym)*ff,(z-zm)*ff}
  442. endfunction
  443.  
  444. function framedsolid (x,y,z,scale=0)
  445. ## works like solid and draws a frame around the plot
  446. ## if scale is specified, then the plot is scaled to fit into a cube of
  447. ## side length 2f centered at 0
  448.     frame=getframe(x,y,z);
  449.     if !holding(); clg; endif;
  450.     h=holding(1);
  451.     if scale; {x1,y1,z1}=scaleframe(x,y,z,frame,scale);
  452.     else; {x1,y1,z1}={x,y,z}; endif;
  453.     color(2); frame1(frame);
  454.     color(1); solid(x1,y1,z1);
  455.     color(2); frame2(frame);
  456.     color(1);
  457.     holding(h);
  458.     return frame
  459. endfunction
  460.  
  461. function framedsolidhue (x,y,z,hue,scale=0,f=1)
  462. ## works like solidhue and draws a frame around the plot
  463. ## if scale is specified, then the plot is scaled to fit into a cube of
  464. ## side length 2f centered at 0
  465.     frame=getframe(x,y,z);
  466.     if !holding(); clg; endif;
  467.     h=holding(1);
  468.     if scale; {x1,y1,z1}=scaleframe(x,y,z,frame,scale);
  469.     else; {x1,y1,z1}={x,y,z}; endif;
  470.     color(2); frame1(frame);
  471.     color(1);
  472.     if f; solidhue(x1,y1,z1,hue,f);
  473.     else; solidhue(x1,y1,z1,hue);
  474.     endif;
  475.     color(2); frame2(frame);
  476.     color(1);
  477.     holding(h);
  478.     return frame
  479. endfunction
  480.  
  481. function framedwire (x,y,z,scale=0)
  482. ## works like framedsolid
  483.     frame=getframe(x,y,z);
  484.     if !holding(); clg; endif;
  485.     h=holding(1);
  486.     if scale; {x1,y1,z1}=scaleframe(x,y,z,frame,scale);
  487.     else; {x1,y1,z1}={x,y,z}; endif;
  488.     color(2); frame1(frame);
  489.     color(1); wire(x1,y1,z1);
  490.     color(2); frame2(frame);
  491.     color(1);
  492.     holding(h);
  493.     return frame
  494. endfunction
  495.  
  496. function density (x,f)
  497. ## density(x,1) makes density plot of the values in the matrix x
  498. ## scaled to fit into [0,f].
  499.     ma=max(max(x)'); mi=min(min(x)'); h=ma-mi;
  500.     if h~=0; h=1; endif;
  501.     xx=(x-mi)/h*f;
  502.     density(xx);
  503.     return xx;
  504. endfunction
  505.  
  506. function solidhue (x,y,z,h,f)
  507. ## solidhue(x,y,z,h) makes a shaded solid 3D-plot of x,y,z.
  508. ## h is the shading and should run between 0 and 1.
  509. ## f determines, if h is scaled to fit in between 0 and f.
  510.     ma=max(max(h)'); mi=min(min(h)'); delta=ma-mi;
  511.     if delta~=0; delta=1; endif;
  512.     hh=(h-mi)/delta*f*0.9999;
  513.     return solidhue(x,y,z,hh);
  514. endfunction
  515.  
  516. .. ### linear algebra things ###
  517.  
  518. function id (n)
  519. ## id(n) creates a nxn identity matrix.
  520.     return diag([n n],0,1);
  521. endfunction
  522.  
  523. function inv (a)
  524. ## inv(a) produces the inverse of a matrix.
  525.     return a\id(cols(a));
  526. endfunction
  527.  
  528. function fit (a,b)
  529. ## fit(a,b) computes x such that ||a.x-b||_2 is minimal.
  530.     A=conj(a');
  531.     return symmult(A,a)\(A.b')
  532. endfunction
  533.  
  534. function kernel (a)
  535. ## kernel(a) computes the kernel of the quadratic matrix a.
  536.     {aa,r,c,d}=lu(a);
  537.     n=size(a); nr=size(r);
  538.     if nr[2]==n[2]; return zeros([1,n[2]])'; endif;
  539.     if nr[2]==0; return id(n[2]); endif;
  540.     c1=nonzeros(c); c2=nonzeros(!c);
  541.     b=lusolve(aa[r,c1],a[r,c2]);
  542.     m=size(b);
  543.     x=zeros([n[2] m[2]]);
  544.     x[c1,:]=-b;
  545.     x[c2,:]=id(cols(c2));
  546.     return x
  547. endfunction
  548.  
  549. function image (a)
  550. ## image(a) computes the image of the quadratic matrix a.
  551.     {aa,r,c,d}=lu(a);
  552.     return a[:,nonzeros(c));
  553. endfunction
  554.  
  555. function det (a)
  556. ## det(A) returns the determinant of A
  557.     {aa,r,c,d}=lu(a);
  558.     return d;
  559. endfunction
  560.  
  561. function eigenvalues (a)
  562. ## eigenvalues(A) returns the eigenvalues of A.
  563.     return polysolve(charpoly(a));
  564. endfunction
  565.  
  566. function eigenspace (a,l)
  567. ## eigenspace(A,l) returns the eigenspace of A to the eigenvaue l.
  568.     k=kernel(a-l*id(cols(a)));
  569.     if k==0; error("No eigenvalue found!"); endif;
  570.     si=size(k);
  571.     loop 1 to si[2];
  572.         x=k[:,index()];
  573.         k[:,index()]=x/sqrt(x'.x);
  574.     end;
  575.     return k;
  576. endfunction
  577.  
  578. function eigen (A)
  579. ## eigen(A) returns the eigenvectors and a basis of eigenvectors of A.
  580. ## {l,x,ll}=eigen(A), where l is a vector of eigenvalues,
  581. ## x is a basis of eigenvectors,
  582. ## and ll is a vector of distinct eigenvalues.
  583.     l=eigenvalues(A);
  584.     s=eigenspace(A,l[1]);
  585.     si=size(s); v=dup(l[1],si[2])'; vv=l[1];
  586.     l=eigenremove(l,si[2]);
  587.     repeat;
  588.         if min(size(l))==0; break; endif;
  589.         ll=l[1]; sp=eigenspace(A,ll);
  590.         si=size(sp); l=eigenremove(l,si[2]);
  591.         s=s|sp; v=v|dup(ll,si[2])'; vv=vv|ll;
  592.     end;
  593.     return {v,s,vv}
  594. endfunction
  595.  
  596.  
  597. function eigenspace1
  598. ## eigenspace1(A,l) returns the eigenspace of A to the eigenvalue l.
  599. ## returns {x,l1}, where l1 should be an improvement over l, and
  600. ## x contains the eigenvectors as columns.
  601.     eps=epsilon();
  602.     repeat;
  603.         k=kernel(arg1-arg2*id(cols(arg1)));
  604.         if k==0; else; break; endif;
  605.         if epsilon()>1 break; endif;
  606.         setepsilon(100*epsilon());
  607.     end;
  608.     if k==0; error("No eigenvalue found!"); endif;
  609.     setepsilon(eps);
  610.     {dummy,l}=eigennewton(arg1,k[:,1],arg2);
  611.     eps=epsilon();
  612.     repeat;
  613.         k=kernel(arg1-l*id(cols(arg1)));
  614.         if k==0; else; break; endif;
  615.         if epsilon()>1 break; endif;
  616.         setepsilon(100*epsilon());
  617.     end;
  618.     if k==0; error("No eigenvalue found!"); endif;
  619.     setepsilon(eps);
  620.     si=size(k);
  621.     loop 1 to si[2];
  622.         x=k[:,index()];
  623.         k[:,index()]=x/sqrt(x'.x);
  624.     end;
  625.     return {k,l};
  626. endfunction
  627.  
  628. function eigenremove(l)
  629. ## helping function.
  630.     return l(nonzeros(!(l[1]~=l)))
  631. endfunction
  632.  
  633. function eigen1
  634. ## eigen(A) returns the eigenvectors and a basis of eigenvectors of A.
  635. ## {l,x,ll}=eigen(A), where l is a vector of eigenvalues,
  636. ## x is a basis of eigenvectors,
  637. ## and ll is a vector of distinct eigenvalues.
  638. ## Improved version of eigen.
  639.     l=eigenvalues(A);
  640.     {s,li}=eigenspace1(A,l[1]);
  641.     si=size(s); v=dup(li,si[2])'; vv=li;
  642.     l=eigenremove(l,si[2]);
  643.     repeat;
  644.         if min(size(l))==0; break; endif;
  645.         {sp,li}=eigenspace1(A,l[1]);
  646.         si=size(sp); l=eigenremove(l,si[2]);
  647.         s=s|sp; v=v|dup(li,si[2])'; vv=vv|li;
  648.     end;
  649.     return {v,s,vv}
  650. endfunction
  651.  
  652. function eigennewton
  653. ## eigennewton(a,x,l) does a newton step to improve the eigenvalue l
  654. ## of a and the eigenvector x.
  655. ## returns {x1,l1}.
  656.     a=arg1; x=arg2; x=x/sqrt(x'.x); n=cols(a);
  657.     d=((a-arg3*id(n))|-x)_(2*x'|0);
  658.     b=d\((a.x-arg3*x)_0);
  659.     return {x-b[1:n],arg3-b[n+1]};
  660. endfunction
  661.  
  662. function hilb (n)
  663. ## hilb(n) produces the nxn-Hilbert matrix.
  664.     {i,j}=field(1:n,1:n);
  665.     return 1/(i+j-1);
  666. endfunction
  667.  
  668. .. ### polynomial fit ##
  669.  
  670. function polyfit (xx,yy,n)
  671. ## fit(x,y,degree) fits a polynomial in L_2-norm to (x,y).
  672.     A=ones(size(xx))_dup(xx,n); A=cumprod(A');
  673.     return fit(A,yy)';
  674. endfunction
  675.  
  676. function field (x,y)
  677. ## field(x,y) returns {X,Y} such that the X+i*Y is a rectangular
  678. ## grid in C containing the points x(n)+i*y(m). x and y must be
  679. ## 1xN and 1xM vectors.
  680.     return {dup(x,cols(y)),dup(y',cols(x))};
  681. endfunction
  682.  
  683. function totalsum (A)
  684. ## totalsum(a) computes the sum of all elements of a
  685.     return sum(sum(A)');
  686. endfunction
  687.  
  688. function totalmin (A)
  689. ## Returns the total minimum of the elements of a
  690.     return min(min(A)')'
  691. endfunction
  692.  
  693. function totalmax (A)
  694. ## Returns the total maximum of the elements of a
  695.     return max(max(A)')'
  696. endfunction
  697.  
  698. function sinh
  699. ## sinh(x) = (exp(x)-exp(-x))/2
  700.     h=exp(arg1);
  701.     return (h-1/h)/2;
  702. endfunction
  703.  
  704. function cosh
  705. ## cosh(x) = (exp(x)+exp(-x))/2
  706.     h=exp(arg1);
  707.     return (h+1/h)/2;
  708. endfunction
  709.  
  710. function arsinh
  711. ## arsinh(z) = log(z+sqrt(z^2+1))
  712.     return log(arg1+sqrt(arg1*arg1+1))
  713. endfunction
  714.  
  715. function arcosh
  716. ## arcosh(z) = log(z+(z^2-1))
  717.     return log(arg1+sqrt(arg1*arg1-1))
  718. endfunction
  719.  
  720. function heun (ffunction,t,y0)
  721. ## y=heun("f",t,y0,...) solves the differential equation y'=f(t,y).
  722. ## f(t,y,...) must be a function. 
  723. ## y0 is the starting value.
  724.     n=cols(t);
  725.     y=dup(y0,n);
  726.     loop 1 to n-1;
  727.         h=t[#+1]-t[#];
  728.         yh=y[#]; xh=t[#];
  729.         k1=ffunction(xh,yh,args(4));
  730.         k2=ffunction(xh+h/2,yh+h/2*k1,args(4));
  731.         k3=ffunction(xh+h,yh+2*h*k2-h*k1,args(4));
  732.         y[#+1]=yh+h/6*(k1+4*k2+k3);
  733.     end;
  734.     return y';
  735. endfunction
  736.  
  737. solveeps=epsilon();
  738.  
  739. function bisect (ffunction,a,b)
  740. ## bisect("f",a,b,...) uses the bisection method to find a root of 
  741. ## f in [a,b].
  742.     global solveeps;
  743.     if ffunction(b,args(4))<0;
  744.         b1=a; a1=b;
  745.     else;
  746.         a1=a; b1=b;
  747.     endif;
  748.     if ffunction(b1,args(4))<0 error("No zero in interval"); endif;
  749.     repeat
  750.         m=(a1+b1)/2;
  751.         if abs(a1-b1)<solveeps; return m; endif;
  752.         if ffunction(m,args(4))>0; b1=m; else a1=m; endif;
  753.     end;
  754. endfunction
  755.  
  756. function secant (ffunction,a,b)
  757. ## secant("f",a,b,...) uses the secant method to find a root of f in [a,b]
  758.     global solveeps;
  759.     x0=a; x1=b; y0=ffunction(x0,args(4)); y1=ffunction(x1,args(4));
  760.     repeat
  761.         x2=x1-y1*(x1-x0)/(y1-y0);
  762.         if abs(x2-x1)<solveeps; break; endif;
  763.         x0=x1; y0=y1; x1=x2; y1=ffunction(x2,args(4));
  764.     end;
  765.     return x2
  766. endfunction
  767.  
  768. function simpson (ffunction,a,b,n=50)
  769. ## simpson("f",a,b) or simpson("f",a,b,n,...) integrates f in [a,b] with
  770. ## the simpson method. f must be able to evaluate a vector. ... is passed
  771. ## to f.
  772.     t=linspace(a,b,2*n);
  773.     s=ffunction(t,args(5));
  774.     ff=4-mod(1:2*n+1,2)*2; ff[1]=1; ff[2*n+1]=1;
  775.     return sum(ff*s)/3*(t[2]-t[1]);
  776. endfunction
  777.  
  778. rombergeps=epsilon();
  779.  
  780. function romberg(ffunction,a,b,m=10)
  781. ## romberg(f,a,b) computes the Romberg integral of f in [a,b].
  782. ## romberg(f,a,b,m,...) specifies h=(b-a)/m/2^k for k=1,...
  783. ## ... is passed to f(x,...)
  784.     global rombergeps;
  785.     y=ffunction(linspace(a,b,m),args(5)); h=(b-a)/m;
  786.     y[1]=y[1]/2; y[m+1]=y[m+1]/2; I=sum(y);
  787.     S=I*h; H=h^2; Intalt=S;
  788.     repeat;
  789.         I=I+sum(ffunction(a+h/2:h:b,args(5))); h=h/2;
  790.         S=S|I*h;
  791.         H=H|h^2;
  792.         Int=interpval(H,interp(H,S),0);
  793.         if abs(Int-Intalt)<rombergeps; break; endif;
  794.         Intalt=Int;
  795.     end;
  796.     return Intalt
  797. endfunction
  798.  
  799. gaussz = [ ..
  800. -9.7390652851717172e-0001,
  801. -8.6506336668898451e-0001,
  802. -6.7940956829902441e-0001,
  803. -4.3339539412924719e-0001,
  804. -1.4887433898163121e-0001,
  805.  1.4887433898163121e-0001,
  806.  4.3339539412924719e-0001,
  807.  6.7940956829902440e-0001,
  808.  8.6506336668898451e-0001,
  809.  9.7390652851717172e-0001];
  810. gaussa = [ ..
  811.  6.6671344308688139e-0002,
  812.  1.4945134915058059e-0001,
  813.  2.1908636251598205e-0001,
  814.  2.6926671930999635e-0001,
  815.  2.9552422471475288e-0001,
  816.  2.9552422471475287e-0001,
  817.  2.6926671930999635e-0001,
  818.  2.1908636251598205e-0001,
  819.  1.4945134915058059e-0001,
  820.  6.6671344308688137e-0002];
  821.  
  822. function gauss10 (ffunction,a,b)
  823. ## gauss10("f",a,b,...)
  824. ## evaluates the gauss integration with 10 knots an [a,b]
  825.     global gaussa,gaussz;
  826.     z=a+(gaussz+1)*(b-a)/2;
  827.     return sum(gaussa*ffunction(z,args(4)))*(b-a)/2;
  828. endfunction
  829.  
  830. function gauss (ffunction,a,b,n=1)
  831. ## gauss("f",a,b) gauss integration with 10 knots
  832. ## gauss("f",a,b,n,...) subdivision into n subintervals.
  833. ## a and b may be vectors.
  834.     si=size(a,b); R=zeros(si);
  835.     if n==1;
  836.         loop 1 to prod(si);
  837.             sum=0;
  838.             sum=sum+gauss10(ffunction,a{#},b{#});
  839.             R{#}=sum;
  840.         end;
  841.     else;
  842.         loop 1 to prod(si);
  843.             h=linspace(a{#},b{#},n); sum=0;
  844.             loop 1 to n;
  845.                 sum=sum+gauss10(ffunction,h{#},h{#+1},args(5));
  846.             end;
  847.             R{#}=sum;
  848.         end;
  849.     endif;
  850.     return R
  851. endfunction
  852.  
  853. .. ### Broyden ###
  854.  
  855. broydenmax=50;
  856.  
  857. function broyden (ffunction,xstart,A=0)
  858. ## broyden("f",x) or broyden("f",x,a,...) finds a zero of f.
  859. ## x is the starting value and a is a approximation of the jacobian.
  860. ## if a is 0, it is neglected.
  861. ## ... is passed to f(x,...)
  862.     global solveeps;
  863.     x=xstart; n=cols(x);
  864.     global broydenmax;
  865.     delta=sqrt(solveeps);
  866.     if A==0;
  867.         A=zeros([n n]);
  868.         y=ffunction(x,args(4));
  869.         loop 1 to n;
  870.             x0=x; x0[#]=x0[#]+delta;
  871.             A[:,#]=(ffunction(x0,args(4))-y)'/delta;
  872.         end;
  873.     endif;
  874.     itercount=0;
  875.     y=ffunction(x,args(4));
  876.     repeat
  877.         d=-A\y'; x=x+d';
  878.         y1=ffunction(x,args(4)); q=y1-y;
  879.         A=A+((q'-A.d).d')/(d'.d);
  880.         if (max(abs(d))<solveeps || itercount>broydenmax) break; endif;
  881.         y=y1;
  882.         itercount=itercount+1;
  883.     end;
  884.     if (itercount>broydenmax) error("Iteration failed"); endif;
  885.     return x;
  886. endfunction
  887.  
  888. function map (ffunction,t)
  889. ## map("f",t,...) evaluates the function f at all points of the vector t.
  890. ## usually this is the same as f(t,...), but sometimes functions do not
  891. ## accept vectors as input.
  892.     si=size(t); s=zeros(si);
  893.     loop 1 to prod(si); s{#}=ffunction(t{#},args(3)); end;
  894.     return s;
  895. endfunction
  896.  
  897. itereps=epsilon();
  898.  
  899. function iterate (ffunction,x,n=0)
  900. ## iterate("f",x,n,...) iterate the function f(x,...) n times,
  901. ## starting with the point x.
  902. ## if n is missing or n is 0, then the iteration stops at a fixed point.
  903. ## Returns the value after n iterations, and its history.
  904.     global itereps;
  905.     if n==0;
  906.         R=x; Rold=x;
  907.         repeat
  908.             Rnew=ffunction(Rold,args(4));
  909.             R=R_Rnew;
  910.             if max(abs(Rold-Rnew))<itereps; return Rnew; endif;
  911.             Rold=Rnew;
  912.         end;
  913.     else;
  914.         R=zeros(n,cols(x));
  915.         R[1]=x;
  916.         loop 2 to n; R[#]=ffunction(R[#-1],args(4)); end;
  917.         return {R[n],R'}
  918.     endif;
  919. endfunction    
  920.  
  921. .. ### Remez algorithm ###
  922.  
  923. function remezhelp (x,y,n)
  924. ## {t,d,h}=remezhelp(x,y,deg) is the case, where x are exactly deg+2
  925. ## points.
  926.     d1=interp(x,y); d2=interp(x,mod(1:n+2,2)*2-1);
  927.     h=d1[n+2]/d2[n+2];
  928.     d=d1-d2*h;
  929.     return {x[1:n+1],d[1:n+1],h};
  930. endfunction
  931.  
  932. remezeps=0.00001;
  933.  
  934. function remez(x,y,deg,trace=0)
  935. ## {t,d,h,r}=remez(x,y,deg) computes the divided difference form
  936. ## of the polynomial best approximation to y at x of degree deg.
  937. ## remez(x,y,deg,1) does the same thing with tracing.
  938. ## h is the discrete error (with sign), and x the rightmost point,
  939. ## which is missing in t.
  940.     global remezeps;
  941.     n=cols(x);
  942.     ind=linspace(1,n,deg+1); h=0;
  943.     repeat
  944.         {t,d,hnew}=remezhelp(x[ind],y[ind],deg);
  945.         r=y-interpval(t,d,x); hh=hnew;
  946.         if trace;
  947.             xind=x[ind];
  948.             plot(x,r); xgrid(xind); ygrid([-h,0,h]);
  949.             title("Weiter mit Taste"); wait(180);
  950.         endif;
  951.         indnew=ind; rr=ind;
  952.         rr[1]=r[ind[1]];
  953.         rr[2]=r[ind[deg+2]];
  954.         loop 2 to deg+1;
  955.             e=extrema(r[ind[index()-1]:ind[index()+1]]);
  956.             if (hh>0);
  957.                 indnew[index()]=e[2]+ind[index()-1]; rr[index()]=e[1];
  958.             else
  959.                 indnew[index()]=e[4]+ind[index()-1]; rr[index()]=e[3];
  960.             endif;
  961.             hh=-hh;
  962.         end;
  963.         ind=indnew;
  964.         if (abs(hnew)<(1+remezeps)*abs(h)) break; endif;
  965.         h=hnew;
  966.     end;
  967.     {t,d,h}=remezhelp(x[ind],y[ind],deg);
  968.     if trace;
  969.         xind=x[ind];
  970.         plot(x,y-interpval(t,d,x)); xgrid(xind); ygrid([-h,0,h]);
  971.         title("Weiter mit Taste"); wait(180);
  972.     endif;
  973.     return {t,d,h,x[ind[deg+2]]};
  974. endfunction
  975.  
  976. .. ### Natural spline ###
  977.  
  978. function spline (x,y)
  979. ## spline(x,y) defines the natural Spline at points x(i) with
  980. ## values y(i). It returns the second derivatives at these points.
  981.     n=cols(x);
  982.     h=x[2:n]-x[1:n-1];
  983.     s=h[2:n-1]+h[1:n-2];
  984.     l=h[2:n-1]/s;
  985.     A=diag([n,n],0,2);
  986.     A=setdiag(A,1,0|l);
  987.     A=setdiag(A,-1,1-l|0);
  988.     b=6/s*((y[3:n]-y[2:n-1])/h[2:n-1]-(y[2:n-1]-y[1:n-2])/h[1:n-2]);
  989.     return (A\(0|b|0)')';
  990. endfunction
  991.  
  992. function splineval (x,y,h,t)
  993. ## splineval(x,y,h,t) evaluates the interpolation spline for
  994. ## (x(i),y(i)) with second derivatives h(i) at t.
  995.     p1=find(x,t); p2=p1+1;
  996.     y1=y[p1]; y2=y[p2];
  997.     x1=x[p1]; x2=x[p2]; f=(x2-x1)^2;
  998.     h1=h[p1]*f; h2=h[p2]*f;
  999.     b=h1/2; c=(h2-h1)/6;
  1000.     a=y2-y1-b-c;
  1001.     d=(t-x1)/(x2-x1);
  1002.     return y1+d*(a+d*(b+c*d));
  1003. endfunction
  1004.  
  1005. .. ### use zeros,... the usual way ###
  1006.  
  1007. function ctext
  1008. ## ctext(s,[c,r]) plots the centered string s at the coordinates (c,r).
  1009. ## Also ctext(s,c,r).
  1010.     if argn()==3; return ctext(arg1,[arg2,arg3]); endif;
  1011.     error("Illegal argument number!"),
  1012. endfunction
  1013.  
  1014. function text
  1015. ## text(s,[c,r]) plots the centered string s at the coordinates (c,r).
  1016. ## Also ctext(s,c,r).
  1017.     if argn()==3; return text(arg1,[arg2,arg3]); endif;
  1018.     error("Illegal argument number!"),
  1019. endfunction
  1020.  
  1021. function diag
  1022. ## diag([n,m],k,v) returns a nxm matrix A, containing v on its k-th
  1023. ## diagonal. v may be a vector or a real number. Also diag(n,m,k,v);
  1024. ## diag(A,k) returns the k-th diagonal of A.
  1025.     if argn()==4; return diag([arg1,arg2],arg3,arg4); endif;
  1026.     error("Illegal argument number!"),
  1027. endfunction
  1028.  
  1029. function format
  1030. ## format([n,m]) sets the number output format to m digits and a total
  1031. ## width of n. Also format(n,m);
  1032.     if argn()==2; return format([arg1,arg2]); endif;
  1033.     error("Illegal argument number!"),
  1034. endfunction
  1035.  
  1036. function normal
  1037. ## normal([n,m]) returns a nxm matrix of unit normal distributed random
  1038. ## values. Also normal(n,m).
  1039.     if argn()==2; return normal([arg1,arg2]); endif;
  1040.     error("Illegal argument number!"),
  1041. endfunction
  1042.  
  1043. function random
  1044. ## random([n,m]) returns a nxm matrix of uniformly distributed random 
  1045. ## values in [0,1]. Also random(n,m).
  1046.     if argn()==2; return random([arg1,arg2]); endif;
  1047.     error("Illegal argument number!"),
  1048. endfunction
  1049.  
  1050. function ones
  1051. ## ones([n,m]) returns a nxm matrix with all elements set to 1.
  1052. ## Also ones(n,m).
  1053.     if argn()==2; return ones([arg1,arg2]); endif;
  1054.     error("Illegal argument number!"),
  1055. endfunction
  1056.  
  1057. function zeros
  1058. ## zeros([n,m]) returns a nxm matrix with all elements set to 0.
  1059. ## Also zeros(n,m).
  1060.     if argn()==2; return zeros([arg1,arg2]); endif;
  1061.     error("Illegal argument number!"),
  1062. endfunction
  1063.  
  1064. function matrix
  1065. ## matrix([n,m],x) returns a nxm matrix with all elements set to x.
  1066. ## Also matrix(n,m,x).
  1067.     if argn()==3; return matrix([arg1,arg2],arg3); endif;
  1068.     error("Illegal argument number!"),
  1069. endfunction
  1070.  
  1071. function view
  1072. ## view([distance, tele, angle1, angle2]) sets the perspective for
  1073. ## solid and view. distance is the eye distance, tele a zooming factor.
  1074. ## angle1 is the angle from the negativ y-axis to the positive x-axis.
  1075. ## angle2 is the angle to the positive z-axis (the height of the eye).
  1076. ## Also view(d,t,a1,a2).
  1077. ## view() returns the values of view.
  1078.     if argn()==4; return view([arg1,arg2,arg3,arg4]); endif;
  1079.     error("Illegal argument number!"),
  1080. endfunction
  1081.  
  1082. function window
  1083. ## window([c1,r1,c2,r2]) sets a plotting window. The cooridnates must
  1084. ## be screen coordimates. Also window(c1,r1,c2,r2).
  1085. ## window() returns the active window coordinates.
  1086.     if argn()==4; return window([arg1,arg2,arg3,arg4]); endif;
  1087.     error("Illegal argument number!"),
  1088. endfunction
  1089.  
  1090. .. *** directory ***
  1091.  
  1092. function dir(pattern="*.*")
  1093. ## dir(pattern) displays a directory.
  1094. ## dir() is the same as dir("*.*").
  1095.     s=searchfile(pattern), if stringcompare(s,"")==0; return 0; endif;
  1096.     repeat
  1097.         s=searchfile(), if stringcompare(s,"")==0; return 0; endif;
  1098.     end;
  1099. endfunction
  1100.  
  1101. ə