home *** CD-ROM | disk | FTP | other *** search
/ Crawly Crypt Collection 1 / crawlyvol1.bin / apps / math / euler / source / funcs.c < prev    next >
Text File  |  1993-04-13  |  44KB  |  1,927 lines
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <math.h>
  5. #include <float.h>
  6. #include <limits.h>
  7.  
  8. #include "header.h"
  9. #include "sysdep.h"
  10. #include "funcs.h"
  11. #include "matheh.h"
  12. #include "polynom.h"
  13. #include "helpf.h"
  14. #include "graphics.h"
  15.  
  16. #define wrong_arg() { error=26; output("Wrong argument\n"); return; }
  17.  
  18. char *argname[] =
  19.     { "arg1","arg2","arg3","arg4","arg5","arg6","arg7","arg8","arg9",
  20.         "arg10" } ;
  21. int xors[10];
  22.  
  23. double (*func) (double);
  24.  
  25. void funceval (double *x, double *y)
  26. /* evaluates the function stored in func with pointers. */
  27. {    *y=func(*x);
  28. }
  29.  
  30. void spread1 (double f (double), 
  31.     void fc (double *, double *, double *, double *),
  32.     header *hd)
  33. {    header *result,*st=hd;
  34.     hd=getvalue(hd);
  35.     if (error) return;
  36.     func=f;
  37.     result=map1(funceval,fc,hd);
  38.     if (!error) moveresult(st,result);
  39. }
  40.  
  41. void csin (double *x, double *xi, double *z, double *zi)
  42. {    *z=cosh(*xi)*sin(*x);
  43.     *zi=sinh(*xi)*cos(*x);
  44. }
  45.  
  46. void msin (header *hd) 
  47. {    spread1(sin,csin,hd);
  48. }
  49.  
  50. void ccos (double *x, double *xi, double *z, double *zi)
  51. {    *z=cosh(*xi)*cos(*x);
  52.     *zi=-sinh(*xi)*sin(*x);
  53. }
  54.  
  55. void mcos (header *hd)
  56. {    spread1(cos,ccos,hd);
  57. }
  58.  
  59. void ctan (double *x, double *xi, double *z, double *zi)
  60. {    double s,si,c,ci;
  61.     csin(x,xi,&s,&si); ccos(x,xi,&c,&ci);
  62.     complex_divide(&s,&si,&c,&ci,z,zi);
  63. }
  64.  
  65. double rtan (double x)
  66. {    if (cos(x)==0.0) return 1e10;
  67.     return tan(x);
  68. }
  69.  
  70. void mtan (header *hd)
  71. {    spread1(
  72. #ifdef FLOAT_TEST
  73.     rtan,
  74. #else
  75.     tan,
  76. #endif
  77.     ctan,hd);
  78. }
  79.  
  80. double ratan (double x)
  81. {    if (x<=-M_PI && x>=M_PI) return 1e10;
  82.     else return atan(x);
  83. }
  84.  
  85. void carg (double *x, double *xi, double *z)
  86. {    
  87. #ifdef FLOAT_TEST
  88.     if (*x==0.0 && *xi==0.0) *z=0.0;
  89. #endif
  90.     *z = atan2(*xi,*x);
  91. }
  92.  
  93. double rlog (double x)
  94. {    if (x<=0) { error=1; return 0; }
  95.     else return log(x);
  96. }
  97.  
  98. void cclog (double *x, double *xi, double *z, double *zi)
  99. {    
  100. #ifdef FLOAT_TEST
  101.     *z=rlog(sqrt(*x * *x + *xi * *xi));
  102. #else
  103.     *z=log(sqrt(*x * *x + *xi * *xi));
  104. #endif
  105.     carg(x,xi,zi);
  106. }
  107.  
  108. double rsign (double x)
  109. {    if (x<0) return -1;
  110.     else if (x<=0) return 0;
  111.     else return 1;
  112. }
  113.  
  114. void msign (header *hd)
  115. {    spread1(rsign,0,hd);
  116. }
  117.  
  118. void catan (double *x, double *xi, double *y, double *yi)
  119. {    double h,hi,g,gi,t,ti;
  120.     h=1-*xi; hi=*x; g=1+*xi; gi=-*x;
  121.     complex_divide(&h,&hi,&g,&gi,&t,&ti);
  122.     cclog(&t,&ti,&h,&hi);
  123.     *y=hi/2; *yi=-h/2;
  124. }
  125.  
  126. void matan (header *hd)
  127. {    spread1(
  128. #ifdef FLOAT_TEST
  129.     ratan,
  130. #else
  131.     atan,
  132. #endif
  133.     catan,hd);
  134. }
  135.  
  136. double rasin (double x)
  137. {    if (x<-1 || x>1) { error=1; return 0; }
  138.     else return asin(x);
  139. }
  140.  
  141. void csqrt (double *x, double *xi, double *z, double *zi)
  142. {    double a,r;
  143.     carg(x,xi,&a); a=a/2.0;
  144.     r=sqrt(sqrt(*x * *x + *xi * *xi));
  145.     *z=r*cos(a);
  146.     *zi=r*sin(a);
  147. }
  148.  
  149. void casin (double *x, double *xi, double *y, double *yi)
  150. {    double h,hi,g,gi;
  151.     complex_multiply(x,xi,x,xi,&h,&hi);
  152.     h=1-h; hi=-hi;
  153.     csqrt(&h,&hi,&g,&gi);
  154.     h=-*xi+g; hi=*x+gi;
  155.     cclog(&h,&hi,yi,y);
  156.     *yi=-*yi;
  157. }
  158.  
  159. void masin (header *hd)
  160. {    spread1(
  161. #ifdef FLOAT_TEST
  162.     rasin,
  163. #else
  164.     asin,
  165. #endif
  166.     casin,hd);
  167. }
  168.  
  169. double racos (double x)
  170. {    if (x<-1 || x>1) { error=1; return 0; }
  171.     else return acos(x);
  172. }
  173.  
  174. void cacos (double *x, double *xi, double *y, double *yi)
  175. {    double h,hi,g,gi;
  176.     complex_multiply(x,xi,x,xi,&h,&hi);
  177.     h=1-h; hi=-hi;
  178.     csqrt(&h,&hi,&g,&gi);
  179.     hi=*xi+g; h=*x-gi;
  180.     cclog(&h,&hi,yi,y);
  181.     *yi=-*yi;
  182. }
  183.  
  184. void macos (header *hd)
  185. {    spread1(
  186. #ifdef FLOAT_TEST
  187.     racos,
  188. #else
  189.     acos,
  190. #endif
  191.     cacos,hd);
  192. }
  193.  
  194. void cexp (double *x, double *xi, double *z, double *zi)
  195. {    double r=exp(*x);
  196.     *z=cos(*xi)*r;
  197.     *zi=sin(*xi)*r;
  198. }
  199.  
  200. void mexp (header *hd)
  201. {    spread1(exp,cexp,hd);
  202. }
  203.  
  204. double rarg (double x)
  205. {    if (x>=0) return 0.0;
  206.     else return M_PI;
  207. }
  208.  
  209. void mlog (header *hd)
  210. {    spread1(log,cclog,hd);
  211. }
  212.  
  213. double rsqrt (double x)
  214. {    if (x<0.0) { error=1; return 1e10; }
  215.     else return sqrt(x);
  216. }
  217.  
  218. void msqrt (header *hd)
  219. {    spread1(
  220. #ifdef FLOAT_TEST
  221.     rsqrt,
  222. #else
  223.     sqrt,
  224. #endif
  225.     csqrt,hd);
  226. }
  227.  
  228. void mceil (header *hd)
  229. {    spread1(ceil,0,hd);
  230. }
  231.  
  232. void mfloor (header *hd)
  233. {    spread1(floor,0,hd);
  234. }
  235.  
  236. void cconj (double *x, double *xi, double *z, double *zi)
  237. {    *zi=-*xi; *z=*x;
  238. }
  239.  
  240. double ident (double x)
  241. {    return x;
  242. }
  243.  
  244. void mconj (header *hd)
  245. {    spread1(ident,cconj,hd);
  246. }
  247.  
  248. void spread1r (double f (double), 
  249.     void fc (double *, double *, double *),
  250.     header *hd)
  251. {    header *result,*st=hd;
  252.     hd=getvalue(hd);
  253.     if (error) return;
  254.     func=f;
  255.     result=map1r(funceval,fc,hd);
  256.     if (!error) moveresult(st,result);
  257. }
  258.  
  259. double rnot (double x)
  260. {    if (x!=0.0) return 0.0;
  261.     else return 1.0;
  262. }
  263.  
  264. void cnot (double *x, double *xi, double *r)
  265. {    if (*x==0.0 && *xi==0.0) *r=1.0;
  266.     else *r=0.0;
  267. }
  268.  
  269. void mnot (header *hd)
  270. {    spread1r(rnot,cnot,hd);
  271. }
  272.  
  273. void crealpart (double *x, double *xi, double *z)
  274. {    *z=*x;
  275. }
  276.  
  277. void mre (header *hd)
  278. {    spread1r(ident,crealpart,hd);
  279. }
  280.  
  281. double zero (double x)
  282. {    return 0.0;
  283. }
  284.  
  285. void cimagpart (double *x, double *xi, double *z)
  286. {    *z=*xi;
  287. }
  288.  
  289. void mim (header *hd)
  290. {    spread1r(zero,cimagpart,hd);
  291. }
  292.  
  293. void marg (header *hd)
  294. {    spread1r(rarg,carg,hd);
  295. }
  296.  
  297. void cxabs (double *x, double *xi, double *z)
  298. {    *z=sqrt(*x * *x + *xi * *xi);
  299. }
  300.  
  301. void mabs (header *hd)
  302. {    spread1r(fabs,cxabs,hd);
  303. }
  304.  
  305. void mpi (header *hd)
  306. {    new_real(M_PI,"");
  307. }
  308.  
  309. void margn (header *hd)
  310. {    new_real(actargn,"");
  311. }
  312.  
  313. void mtime (header *hd)
  314. {    hd=new_real(myclock(),"");
  315. }
  316.  
  317. void mfree (header *hd)
  318. {    new_real(ramend-endlocal,"");
  319. }
  320.  
  321. void mshrink (header *hd)
  322. {    header *st=hd,*result;
  323.     size_t size;
  324.     hd=getvalue(hd); if (error) return;
  325.     if (*realof(hd)>LONG_MAX) wrong_arg();
  326.     size=(size_t)*realof(hd);
  327.     if (ramend-size<newram)
  328.     {    output("Cannot shrink that much!\n");
  329.         error=171; return;
  330.     }
  331.     if (size) 
  332.         if (shrink(size)) ramend-=size;
  333.         else
  334.         {    output("Shrink failed!\n"); error=172; return;
  335.         }
  336.     result=new_real(ramend-ramstart,"");
  337.     moveresult(st,result);
  338. }
  339.  
  340. void mepsilon (header *hd)
  341. {    new_real(epsilon,"");
  342. }
  343.  
  344. void msetepsilon (header *hd)
  345. {    header *stack=hd,*hd1,*result;
  346.     hd1=getvalue(hd); if (error) return;
  347.     if (hd1->type!=s_real) wrong_arg();
  348.     result=new_real(epsilon,"");
  349.     epsilon=*realof(hd1);
  350.     moveresult(stack,result);
  351. }    
  352.  
  353. void mindex (header *hd)
  354. {    new_real((double)loopindex,"");
  355. }
  356.  
  357. void spread2 (void f (double *, double *, double *),
  358.     void fc (double *, double *, double *, double *, double *, double *),
  359.     header *hd)
  360. {    header *result,*st=hd,*hd1;
  361.     hd=getvalue(hd); if (error) return;
  362.     hd1=next_param(st); if (!hd1 || error) return;
  363.     hd1=getvalue(hd1); if (error) return;
  364.     result=map2(f,fc,hd,hd1);
  365.     if (!error) moveresult(st,result);
  366. }
  367.  
  368. void spread2r (void f (double *, double *, double *),
  369.     void fc (double *, double *, double *, double *, double *, double *),
  370.     header *hd)
  371. {    header *result,*st=hd,*hd1;
  372.     int complex;
  373.     hd=getvalue(hd); if (error) return;
  374.     hd1=next_param(st); if (!hd1 || error) return;
  375.     hd1=getvalue(hd1); if (error) return;
  376.     complex=(hd1->type==s_complex || hd1->type==s_cmatrix ||
  377.         hd->type==s_complex || hd->type==s_cmatrix);
  378.     result=map2(f,fc,hd,hd1);
  379.     if (complex) mcomplex(result);
  380.     if (!error) moveresult(st,result);
  381. }
  382.  
  383. void rmod (double *x, double *n, double *y)
  384. {    *y=fmod(*x,*n);
  385. }
  386.  
  387. void mmod (header *hd)
  388. {    spread2(rmod,0,hd);
  389. }
  390.  
  391. void cpow (double *x, double *xi, double *y, double *yi,
  392.     double *z, double *zi)
  393. {    double l,li,w,wi;
  394.     if (fabs(*x)<epsilon && fabs(*xi)<epsilon)
  395.     {    *z=*zi=0.0; return;
  396.     }
  397.     cclog(x,xi,&l,&li);
  398.     complex_multiply(y,yi,&l,&li,&w,&wi);
  399.     cexp(&w,&wi,z,zi);
  400. }
  401.  
  402. void rpow (double *x, double *y, double *z)
  403. {    int n;
  404.     if (*x>0.0) *z=pow(*x,*y);
  405.     else if (*x==0.0) if (*y==0.0) *z=1.0; else *z=0.0;
  406.     else
  407.     {    n=(int)*y;
  408.         if (n%2) *z=-pow(-*x,n);
  409.         else *z=pow(-*x,n);
  410.     }
  411. }
  412.  
  413. void mpower (header *hd)
  414. {    spread2(rpow,cpow,hd);
  415. }
  416.  
  417. void rgreater (double *x, double *y, double *z)
  418. {    if (*x>*y) *z=1.0;
  419.     else *z=0.0;
  420. }
  421.  
  422. void mgreater (header *hd)
  423. {    spread2(rgreater,0,hd);
  424. }
  425.  
  426. void rless (double *x, double *y, double *z)
  427. {    if (*x<*y) *z=1.0;
  428.     else *z=0.0;
  429. }
  430.  
  431. void mless (header *hd)
  432. {    spread2(rless,0,hd);
  433. }
  434.  
  435. void rgreatereq (double *x, double *y, double *z)
  436. {    if (*x>=*y) *z=1.0;
  437.     else *z=0.0;
  438. }
  439.  
  440. void mgreatereq (header *hd)
  441. {    spread2(rgreatereq,0,hd);
  442. }
  443.  
  444. void rlesseq (double *x, double *y, double *z)
  445. {    if (*x<=*y) *z=1.0;
  446.     else *z=0.0;
  447. }
  448.  
  449. void mlesseq (header *hd)
  450. {    spread2(rlesseq,0,hd);
  451. }
  452.  
  453. void ror (double *x, double *y, double *z)
  454. {    if (*x!=0.0 || *y!=0.0) *z=1.0;
  455.     else *z=0.0;
  456. }
  457.  
  458. void mor (header *hd)
  459. {    spread2(ror,0,hd);
  460. }
  461.  
  462. void rrand (double *x, double *y, double *z)
  463. {    if (*x!=0.0 && *y!=0.0) *z=1.0;
  464.     else *z=0.0;
  465. }
  466.  
  467. void mand (header *hd)
  468. {    spread2(rrand,0,hd);
  469. }
  470.  
  471. void requal (double *x, double *y, double *z)
  472. {    if (*x==*y) *z=1.0;
  473.     else *z=0.0;
  474. }
  475.  
  476. void cequal (double *x, double *xi, double *y, double *yi, double *z,
  477.     double *zi)
  478. {    if (*x==*xi && *y==*yi) *z=1.0;
  479.     else *z=0.0;
  480.     *zi=0.0;
  481. }
  482.  
  483. void mequal (header *hd)
  484. {    spread2r(requal,cequal,hd);
  485. }
  486.  
  487. void runequal (double *x, double *y, double *z)
  488. {    if (*x!=*y) *z=1.0;
  489.     else *z=0.0;
  490. }
  491.  
  492. void cunequal (double *x, double *xi, double *y, double *yi, double *z,
  493.     double *zi)
  494. {    if (*x!=*y || *xi!=*yi) *z=1.0;
  495.     else *z=0.0;
  496.     *zi=0.0;
  497. }
  498.  
  499. void munequal (header *hd)
  500. {    spread2(runequal,cunequal,hd);
  501. }
  502.  
  503. void raboutequal (double *x, double *y, double *z)
  504. {    if (fabs(*x-*y)<epsilon) *z=1.0;
  505.     else *z=0.0;
  506. }
  507.  
  508. void caboutequal 
  509.     (double *x, double *xi, double *y, double *yi, double *z,
  510.         double *zi)
  511. {    if (fabs(*x-*y)<epsilon && fabs(*xi-*yi)<epsilon) *z=1.0;
  512.     else *z=0.0;
  513.     *zi=0.0;
  514. }
  515.  
  516. void maboutequal (header *hd)
  517. {    spread2r(raboutequal,caboutequal,hd);
  518. }
  519.  
  520. void mlusolve (header *hd)
  521. {    header *st=hd,*hd1,*result;
  522.     double *m,*m1;
  523.     int r,c,r1,c1;
  524.     hd=getvalue(hd);
  525.     hd1=next_param(st);
  526.     if (hd1) hd1=getvalue(hd1);
  527.     if (error) return;
  528.     if (hd->type==s_matrix || hd->type==s_real)
  529.     {    getmatrix(hd,&r,&c,&m);
  530.         if (hd1->type==s_cmatrix)
  531.         {    make_complex(st);
  532.             mlusolve(st); return;    
  533.         }
  534.         if (hd1->type!=s_matrix && hd1->type!=s_real) wrong_arg();
  535.         getmatrix(hd1,&r1,&c1,&m1);
  536.         if (c!=r || c<1 || r!=r1) wrong_arg();
  537.         result=new_matrix(r,c1,""); if (error) return;
  538.         lu_solve(m,r,m1,c1,matrixof(result));
  539.         if (error) return;
  540.         moveresult(st,result);
  541.     }
  542.     else if (hd->type==s_cmatrix || hd->type==s_complex)
  543.     {    getmatrix(hd,&r,&c,&m);
  544.         if (hd1->type==s_matrix || hd1->type==s_real)
  545.         {    make_complex(next_param(st));
  546.             mlusolve(st); return;
  547.         }
  548.         if (hd1->type!=s_cmatrix && hd1->type!=s_complex) wrong_arg();
  549.         getmatrix(hd1,&r1,&c1,&m1);
  550.         if (c!=r || c<1 || r!=r1) wrong_arg();
  551.         result=new_cmatrix(r,c1,""); if (error) return;
  552.         clu_solve(m,r,m1,c1,matrixof(result));
  553.         if (error) return;
  554.         moveresult(st,result);
  555.     }
  556.     else wrong_arg();
  557. }
  558.  
  559. void msolve (header *hd)
  560. {    header *st=hd,*hd1,*result;
  561.     double *m,*m1;
  562.     int r,c,r1,c1;
  563.     hd=getvalue(hd);
  564.     hd1=next_param(st);
  565.     if (hd1) hd1=getvalue(hd1);
  566.     if (error) return;
  567.     if (hd->type==s_matrix || hd->type==s_real)
  568.     {    getmatrix(hd,&r,&c,&m);
  569.         if (hd1->type==s_cmatrix)
  570.         {    make_complex(st);
  571.             msolve(st); return;    
  572.         }
  573.         if (hd1->type!=s_matrix && hd1->type!=s_real) wrong_arg();
  574.         getmatrix(hd1,&r1,&c1,&m1);
  575.         if (c!=r || c<1 || r!=r1) wrong_arg();
  576.         result=new_matrix(r,c1,""); if (error) return;
  577.         solvesim(m,r,m1,c1,matrixof(result));
  578.         if (error) return;
  579.         moveresult(st,result);
  580.     }
  581.     else if (hd->type==s_cmatrix || hd->type==s_complex)
  582.     {    getmatrix(hd,&r,&c,&m);
  583.         if (hd1->type==s_matrix || hd1->type==s_real)
  584.         {    make_complex(next_param(st));
  585.             msolve(st); return;
  586.         }
  587.         if (hd1->type!=s_cmatrix && hd1->type!=s_complex) wrong_arg();
  588.         getmatrix(hd1,&r1,&c1,&m1);
  589.         if (c!=r || c<1 || r!=r1) wrong_arg();
  590.         result=new_cmatrix(r,c1,""); if (error) return;
  591.         c_solvesim(m,r,m1,c1,matrixof(result));
  592.         if (error) return;
  593.         moveresult(st,result);
  594.     }
  595.     else wrong_arg();
  596. }
  597.  
  598. void mcomplex (header *hd)
  599. {    header *st=hd,*result;
  600.     double *m,*mr;
  601.     LONG i,n;
  602.     int c,r;
  603.     hd=getvalue(hd);
  604.     if (hd->type==s_matrix)
  605.     {    getmatrix(hd,&r,&c,&m);
  606.         result=new_cmatrix(r,c,""); if (error) return;
  607.         n=(LONG)r*c;
  608.         mr=matrixof(result)+(LONG)2*(n-1);
  609.         m+=n-1;
  610.         for (i=0; i<n; i++)
  611.         {    *mr=*m--; *(mr+1)=0.0; mr-=2;
  612.         }
  613.         moveresult(st,result);
  614.     }
  615.     else if (hd->type==s_real)
  616.     {    result=new_complex(*realof(hd),0.0,""); if (error) return;
  617.         moveresult(st,result);
  618.     }
  619. }
  620.  
  621. void msum (header *hd)
  622. {    header *st=hd,*result;
  623.     int c,r,i,j;
  624.     double *m,*mr,s,si;
  625.     hd=getvalue(hd); if (error) return;
  626.     if (hd->type==s_real || hd->type==s_matrix)
  627.     {    getmatrix(hd,&r,&c,&m);
  628.         result=new_matrix(r,1,""); if (error) return;
  629.         mr=matrixof(result);
  630.         for (i=0; i<r; i++) 
  631.         {    s=0.0;
  632.             for (j=0; j<c; j++) s+=*m++;
  633.             *mr++=s;
  634.         }
  635.     }
  636.     else if (hd->type==s_complex || hd->type==s_cmatrix)
  637.     {    getmatrix(hd,&r,&c,&m);
  638.         result=new_cmatrix(r,1,""); if (error) return;
  639.         mr=matrixof(result);
  640.         for (i=0; i<r; i++) 
  641.         {    s=0.0; si=0.0;
  642.             for (j=0; j<c; j++) { s+=*m++; si+=*m++; }
  643.             *mr++=s; *mr++=si;
  644.         }
  645.     }
  646.     else wrong_arg();
  647.     moveresult(st,result);
  648. }
  649.  
  650. void mprod (header *hd)
  651. {    header *st=hd,*result;
  652.     int c,r,i,j;
  653.     double *m,*mr,s,si,h,hi;
  654.     hd=getvalue(hd); if (error) return;
  655.     if (hd->type==s_real || hd->type==s_matrix)
  656.     {    getmatrix(hd,&r,&c,&m);
  657.         result=new_matrix(r,1,""); if (error) return;
  658.         mr=matrixof(result);
  659.         for (i=0; i<r; i++) 
  660.         {    s=1.0;
  661.             for (j=0; j<c; j++) s*=*m++;
  662.             *mr++=s;
  663.         }
  664.     }
  665.     else if (hd->type==s_complex || hd->type==s_cmatrix)
  666.     {    getmatrix(hd,&r,&c,&m);
  667.         result=new_cmatrix(r,1,""); if (error) return;
  668.         mr=matrixof(result);
  669.         for (i=0; i<r; i++) 
  670.         {    s=1.0;
  671.             for (j=0; j<c; j++) 
  672.             {    complex_multiply(&s,&si,m,m+1,&h,&hi);
  673.                 s=h; si=hi; m+=2; 
  674.             }
  675.             *mr++=s; *mr++=si;
  676.         }
  677.     }
  678.     else wrong_arg();
  679.     moveresult(st,result);
  680. }
  681.  
  682. void msize (header *hd)
  683. {    header *result,*st=hd,*hd1=hd,*end=(header *)newram;
  684.     int r,c,r0=0,c0=0;
  685.     if (!hd) wrong_arg();
  686.     result=new_matrix(1,2,""); if (error) return;
  687.     while (end>hd)
  688.     {    hd1=getvariable(hd); if (error) return;
  689.         if (hd1->type==s_matrix || hd1->type==s_cmatrix)
  690.         {    r=dimsof(hd1)->r;
  691.             c=dimsof(hd1)->c;
  692.         }
  693.         else if (hd1->type==s_real || hd1->type==s_complex)
  694.         {    r=c=1;
  695.         }
  696.         else if (hd1->type==s_submatrix || hd1->type==s_csubmatrix)
  697.         {    r=submdimsof(hd1)->r;
  698.             c=submdimsof(hd1)->c;
  699.         }
  700.         else wrong_arg();
  701.         if ((r>1 || c>1) && (r0>1 || c0>1))
  702.         {    if (r0!=r && c0!=c)
  703.             {    output("Matrix dimensions must agree!\n");
  704.                 error=1021; return;
  705.             }
  706.         }
  707.         else
  708.         {    r0=(r0>r)?r0:r; c0=(c0>c)?c0:c;
  709.         }
  710.         hd=nextof(hd);
  711.     }
  712.     *matrixof(result)=r0;
  713.     *(matrixof(result)+1)=c0;
  714.     moveresult(st,result);
  715. }
  716.  
  717. void mcols (header *hd)
  718. {    header *st=hd,*res;
  719.     int n;
  720.     hd=getvalue(hd); if (error) return;
  721.     switch (hd->type)
  722.     {    case s_matrix :
  723.         case s_cmatrix : n=dimsof(hd)->c; break;
  724.         case s_submatrix :
  725.         case s_csubmatrix : n=submdimsof(hd)->c; break;
  726.         case s_real :
  727.         case s_complex : n=1; break;
  728.         case s_string : n=(int)strlen(stringof(hd)); break;
  729.         default : wrong_arg();
  730.     }
  731.     res=new_real(n,""); if (error) return;
  732.     moveresult(st,res);
  733. }
  734.  
  735. void mrows (header *hd)
  736. {    header *st=hd,*res;
  737.     int n;
  738.     hd=getvalue(hd); if (error) return;
  739.     switch (hd->type)
  740.     {    case s_matrix :
  741.         case s_cmatrix : n=dimsof(hd)->r; break;
  742.         case s_submatrix :
  743.         case s_csubmatrix : n=submdimsof(hd)->r; break;
  744.         case s_real :
  745.         case s_complex : n=1; break;
  746.         default : wrong_arg();
  747.     }
  748.     res=new_real(n,""); if (error) return;
  749.     moveresult(st,res);
  750. }
  751.  
  752. void mzerosmat (header *hd)
  753. {    header *result,*st=hd;
  754.     double rows,cols,*m;
  755.     int r,c;
  756.     LONG i,n;
  757.     hd=getvalue(hd); if (error) return;
  758.     if (hd->type!=s_matrix || dimsof(hd)->r!=1 || dimsof(hd)->c!=2)
  759.         wrong_arg();
  760.     rows=*matrixof(hd); cols=*(matrixof(hd)+1);
  761.     if (rows<0 || rows>=INT_MAX || cols<0 || cols>=INT_MAX)
  762.         wrong_arg();
  763.     r=(int)rows; c=(int)cols;
  764.     result=new_matrix(r,c,""); if (error) return;
  765.     m=matrixof(result);
  766.     n=c*r;
  767.     for (i=0; i<n; i++) *m++=0.0;
  768.     moveresult(st,result);
  769. }
  770.  
  771. void mones (header *hd)
  772. {    header *result,*st=hd;
  773.     double rows,cols,*m;
  774.     int r,c;
  775.     LONG i,n;
  776.     hd=getvalue(hd); if (error) return;
  777.     if (hd->type!=s_matrix || dimsof(hd)->r!=1 || dimsof(hd)->c!=2)
  778.         wrong_arg();
  779.     rows=*matrixof(hd); cols=*(matrixof(hd)+1);
  780.     if (rows<0 || rows>=INT_MAX || cols<0 || cols>=INT_MAX)
  781.         wrong_arg();
  782.     r=(int)rows; c=(int)cols;
  783.     result=new_matrix(r,c,""); if (error) return;
  784.     m=matrixof(result);
  785.     n=c*r;
  786.     for (i=0; i<n; i++) *m++=1.0;
  787.     moveresult(st,result);
  788. }
  789.  
  790. void mdiag (header *hd)
  791. {    header *result,*st=hd,*hd1,*hd2=0;
  792.     double rows,cols,*m,*md;
  793.     int r,c,i,ik=0,k,rd,cd;
  794.     LONG l,n;
  795.     hd=getvalue(hd); if (error) return;
  796.     if (hd->type!=s_matrix || dimsof(hd)->r!=1 || dimsof(hd)->c!=2)
  797.         wrong_arg();
  798.     rows=*matrixof(hd); cols=*(matrixof(hd)+1);
  799.     if (rows<0 || rows>=INT_MAX || cols<0 || cols>=INT_MAX)
  800.         wrong_arg();
  801.     r=(int)rows; c=(int)cols;
  802.     hd1=next_param(st); if (hd1) hd2=next_param(hd1);
  803.     if (hd1) hd1=getvalue(hd1);
  804.     if (hd2) hd2=getvalue(hd2);
  805.     if (error) return;
  806.     if    (hd1->type!=s_real) wrong_arg();
  807.     k=(int)*realof(hd1);
  808.     if (hd2->type==s_matrix || hd2->type==s_real)
  809.     {    result=new_matrix(r,c,""); if (error) return;
  810.         m=matrixof(result);
  811.         n=(LONG)c*r;
  812.         for (l=0; l<n; l++) *m++=0.0;
  813.         getmatrix(hd2,&rd,&cd,&md);
  814.         if (rd!=1 || cd<1) wrong_arg();
  815.         m=matrixof(result);
  816.         for (i=0; i<r; i++)
  817.         {    if (i+k>=0 && i+k<c) 
  818.             {    *mat(m,c,i,i+k)=*md;
  819.                 ik++; if (ik<cd) md++;
  820.             }
  821.         }
  822.     }
  823.     else if (hd2->type==s_cmatrix || hd2->type==s_complex)
  824.     {    result=new_cmatrix(r,c,""); if (error) return;
  825.         m=matrixof(result);
  826.         n=(LONG)2*(LONG)c*r;
  827.         for (l=0; l<n; l++) *m++=0.0;
  828.         getmatrix(hd2,&rd,&cd,&md);
  829.         if (rd!=1 || cd<1) wrong_arg();
  830.         m=matrixof(result);
  831.         for (i=0; i<r; i++)
  832.         {    if (i+k>=0 && i+k<c) 
  833.             {    *cmat(m,c,i,i+k)=*md;
  834.                 *(cmat(m,c,i,i+k)+1)=*(md+1);
  835.                 ik++; if (ik<cd) md+=2;
  836.             }
  837.         }
  838.     }
  839.     else wrong_arg();
  840.     moveresult(st,result);
  841. }
  842.  
  843. void msetdiag (header *hd)
  844. {    header *result,*st=hd,*hd1,*hd2=0;
  845.     double *m,*md,*mhd;
  846.     int r,c,i,ik=0,k,rd,cd;
  847.     hd=getvalue(st); if (error) return;
  848.     if (hd->type!=s_matrix && hd->type!=s_cmatrix)
  849.         wrong_arg();
  850.     getmatrix(hd,&c,&r,&mhd);
  851.     hd1=next_param(st); if (hd1) hd2=next_param(hd1);
  852.     if (hd1) hd1=getvalue(hd1);
  853.     if (hd2) hd2=getvalue(hd2);
  854.     if (error) return;
  855.     if    (hd1->type!=s_real) wrong_arg();
  856.     k=(int)*realof(hd1);
  857.     if (hd->type==s_matrix && 
  858.             (hd2->type==s_complex || hd2->type==s_cmatrix))
  859.         {    make_complex(st); msetdiag(st); return;
  860.         }
  861.     else if (hd->type==s_cmatrix &&
  862.             (hd2->type==s_real || hd2->type==s_matrix))
  863.         {    make_complex(nextof(nextof(st))); msetdiag(st); return;
  864.         }
  865.     if (hd->type==s_matrix)
  866.     {    result=new_matrix(r,c,""); if (error) return;
  867.         m=matrixof(result);
  868.         memmove((char *)m,(char *)mhd,(LONG)c*r*sizeof(double));
  869.         getmatrix(hd2,&rd,&cd,&md);
  870.         if (rd!=1 || cd<1) wrong_arg();
  871.         for (i=0; i<r; i++)
  872.         {    if (i+k>=0 && i+k<c) 
  873.             {    *mat(m,c,i,i+k)=*md;
  874.                 ik++; if (ik<cd) md++;
  875.             }
  876.         }
  877.     }
  878.     else if (hd->type==s_cmatrix)
  879.     {    result=new_cmatrix(r,c,""); if (error) return;
  880.         m=matrixof(result);
  881.         memmove((char *)m,(char *)mhd,(LONG)c*r*(LONG)2*sizeof(double));
  882.         getmatrix(hd2,&rd,&cd,&md);
  883.         if (rd!=1 || cd<1) wrong_arg();
  884.         m=matrixof(result);
  885.         for (i=0; i<r; i++)
  886.         {    if (i+k>=0 && i+k<c) 
  887.             {    *cmat(m,c,i,i+k)=*md;
  888.                 *(cmat(m,c,i,i+k)+1)=*(md+1);
  889.                 ik++; if (ik<cd) md+=2;
  890.             }
  891.         }
  892.     }
  893.     else wrong_arg();
  894.     moveresult(st,result);
  895. }
  896.  
  897. void mextrema (header *hd)
  898. {    header *result,*st=hd;
  899.     double x,*m,*mr,min,max;
  900.     int r,c,i,j,imin,imax;
  901.     hd=getvalue(hd); if (error) return;
  902.     if (hd->type==s_real || hd->type==s_matrix)
  903.     {    getmatrix(hd,&r,&c,&m);
  904.         result=new_matrix(r,4,""); if (error) return;
  905.         mr=matrixof(result);
  906.         for (i=0; i<r; i++) 
  907.         {    min=max=*m; imin=imax=0; m++;
  908.             for (j=1; j<c; j++) 
  909.             {    x=*m++;
  910.                 if (x<min) { min=x; imin=j; }
  911.                 if (x>max) { max=x; imax=j; }
  912.             }
  913.             *mr++=min; *mr++=imin+1; *mr++=max; *mr++=imax+1;
  914.         }
  915.     }
  916.     else wrong_arg();
  917.     moveresult(st,result);
  918. }
  919.  
  920. void mcumsum (header *hd)
  921. {    header *result,*st=hd;
  922.     double *m,*mr,sum=0,sumr=0,sumi=0;
  923.     int r,c,i,j;
  924.     hd=getvalue(hd); if (error) return;
  925.     if (hd->type==s_real || hd->type==s_matrix)
  926.     {    getmatrix(hd,&r,&c,&m);
  927.         if (c<1) result=new_matrix(r,1,"");
  928.         else result=new_matrix(r,c,"");
  929.         if (error) return;
  930.         mr=matrixof(result);
  931.         for (i=0; i<r; i++) 
  932.         {    if (c>=1) sum=*m++;
  933.             *mr++=sum;
  934.             for (j=1; j<c; j++) 
  935.             {    sum+=*m++;
  936.                 *mr++=sum;
  937.             }
  938.         }
  939.     }
  940.     else if (hd->type==s_complex || hd->type==s_cmatrix)
  941.     {    getmatrix(hd,&r,&c,&m);
  942.         if (c<1) result=new_cmatrix(r,1,"");
  943.         else result=new_cmatrix(r,c,"");
  944.         if (error) return;
  945.         mr=matrixof(result);
  946.         for (i=0; i<r; i++) 
  947.         {    if (c>=1) { sumr=*m++; sumi=*m++; }
  948.             *mr++=sumr; *mr++=sumi;
  949.             for (j=1; j<c; j++) 
  950.             {    sumr+=*m++; *mr++=sumr;
  951.                 sumi+=*m++; *mr++=sumi;
  952.             }
  953.         }
  954.     }
  955.     else wrong_arg();
  956.     moveresult(st,result);
  957. }
  958.  
  959. void mcumprod (header *hd)
  960. {    header *result,*st=hd;
  961.     double *m,*mr,sum=1,sumi=1,sumr=0;
  962.     int r,c,i,j;
  963.     hd=getvalue(hd); if (error) return;
  964.     if (hd->type==s_real || hd->type==s_matrix)
  965.     {    getmatrix(hd,&r,&c,&m);
  966.         if (c<1) result=new_matrix(r,1,"");
  967.         else result=new_matrix(r,c,"");
  968.         if (error) return;
  969.         mr=matrixof(result);
  970.         for (i=0; i<r; i++) 
  971.         {    if (c>=1) sum=*m++; 
  972.             *mr++=sum;
  973.             for (j=1; j<c; j++) 
  974.             {    sum*=*m++;
  975.                 *mr++=sum;
  976.             }
  977.         }
  978.     }
  979.     else if (hd->type==s_complex || hd->type==s_cmatrix)
  980.     {    getmatrix(hd,&r,&c,&m);
  981.         if (c<1) result=new_cmatrix(r,1,"");
  982.         else result=new_cmatrix(r,c,"");
  983.         if (error) return;
  984.         mr=matrixof(result);
  985.         for (i=0; i<r; i++) 
  986.         {    if (c>=1) { sumr=*m++; sumi=*m++; }
  987.             *mr++=sumr; *mr++=sumi;
  988.             for (j=1; j<c; j++)
  989.             {    sum=sumr*(*m)-sumi*(*(m+1));
  990.                 sumi=sumr*(*(m+1))+sumi*(*m);
  991.                 sumr=sum;
  992.                 m+=2;
  993.                 *mr++=sumr;
  994.                 *mr++=sumi;
  995.             }
  996.         }
  997.     }
  998.     else wrong_arg();
  999.     moveresult(st,result);
  1000. }
  1001.  
  1002. void mwait (header *hd)
  1003. {    header *st=hd,*result;
  1004.     double now;
  1005.     int h;
  1006.     hd=getvalue(hd); if (error) return;
  1007.     if (hd->type!=s_real) wrong_arg();
  1008.     now=myclock();
  1009.     sys_wait(*realof(hd),&h);
  1010.     if (h==escape) { error=1; return; }
  1011.     result=new_real(myclock()-now,"");
  1012.     if (error) return;
  1013.     moveresult(st,result);
  1014. }
  1015.  
  1016. void mkey (header *hd)
  1017. {    int scan;
  1018.     wait_key(&scan);
  1019.     new_real(scan,"");
  1020. }
  1021.  
  1022. void mformat (header *hd)
  1023. {    header *st=hd,*result;
  1024.     static int l=10,d=5;
  1025.     int oldl=l,oldd=d;
  1026.     hd=getvalue(hd); if (error) return;
  1027.     if (hd->type!=s_matrix || dimsof(hd)->r!=1 || dimsof(hd)->c!=2)
  1028.         wrong_arg();
  1029.     l=(int)*matrixof(hd); d=(int)*(matrixof(hd)+1);
  1030.     if (l<2 || l>2*DBL_DIG || d<0 || d>DBL_DIG) wrong_arg();
  1031.     if (d>l-3) d=l-3;
  1032.     sprintf(fixedformat," %c%d.%df",'%',l,d);
  1033.     sprintf(expoformat," %c%d.%de",'%',l,(l>9)?l-9:0);
  1034.     minexpo=pow(10,-d);
  1035.     maxexpo=pow(10,(l-d-3>DBL_DIG-d-1)?DBL_DIG-d-1:l-d-3)-1;
  1036.     fieldw=l+2;
  1037.     linew=linelength/fieldw;
  1038.     result=new_matrix(1,2,""); if (error) return;
  1039.     *matrixof(result)=oldl;
  1040.     *(matrixof(result)+1)=oldd;
  1041.     moveresult(st,result);
  1042. }
  1043.  
  1044. void mrandom (header *hd)
  1045. {    header *st=hd,*result;
  1046.     double *m;
  1047.     int r,c;
  1048.     LONG k,n;
  1049.     hd=getvalue(hd); if (error) return;
  1050.     if (hd->type!=s_matrix || dimsof(hd)->r!=1 || dimsof(hd)->c!=2
  1051.         || *(m=matrixof(hd))<0 || *m>=INT_MAX 
  1052.         || *(m+1)<0 || *(m+1)>INT_MAX)
  1053.         wrong_arg();
  1054.     r=(int)*m;
  1055.     c=(int)*(m+1);
  1056.     result=new_matrix(r,c,""); if (error) return;
  1057.     m=matrixof(result);
  1058.     n=(LONG)c*r;
  1059.     for (k=0; k<n; k++) *m++=(double)rand()/(double)RAND_MAX;
  1060.     moveresult(st,result);
  1061. }
  1062.  
  1063. void mnormal (header *hd)
  1064. {    header *st=hd,*result;
  1065.     double *m,r1,r2;
  1066.     int r,c;
  1067.     LONG k,n;
  1068.     hd=getvalue(hd); if (error) return;
  1069.     if (hd->type!=s_matrix || dimsof(hd)->r!=1 || dimsof(hd)->c!=2
  1070.         || *(m=matrixof(hd))<0 || *m>=INT_MAX 
  1071.         || *(m+1)<0 || *(m+1)>INT_MAX)
  1072.         wrong_arg();
  1073.     r=(int)*m;
  1074.     c=(int)*(m+1);
  1075.     result=new_matrix(r,c,""); if (error) return;
  1076.     m=matrixof(result);
  1077.     n=(LONG)c*r;
  1078.     for (k=0; k<n; k++) 
  1079.     {    r1=(double)rand()/(double)RAND_MAX;
  1080.         if (r1==0.0) *m++=0.0;
  1081.         else
  1082.         {    r2=(double)rand()/(double)RAND_MAX;
  1083.             *m++=sqrt(-2*log(r1))*cos(2*M_PI*r2);
  1084.         }
  1085.     }
  1086.     moveresult(st,result);
  1087. }
  1088.  
  1089. double gauss (double z)
  1090. {    double x,w;
  1091.     x=1/(0.2316419*fabs(z)+1);
  1092.     w=x*(0.31938153+x*(-0.356556382+x*(1.781477937+x*(
  1093.         -1.821255978+x*1.330274429))));
  1094.     w=w*exp(-z*z/2)/sqrt(2*M_PI);
  1095.     if (z<0) return w;
  1096.     else return 1-w;
  1097. }
  1098.  
  1099. void mgauss (header *hd)
  1100. {    spread1(gauss,0,hd);
  1101. }
  1102.  
  1103. double invgauss (double a)
  1104. {    double t,c,d;
  1105.     int flag=0;
  1106.     if (a<0.5) { a=1-a; flag=1; }
  1107.     t=sqrt(-2*log(fabs(1-a)));
  1108.     c=2.515517+t*(0.802853+t*0.010328);
  1109.     d=1+t*(1.432788+t*(0.189269+t*0.001308));
  1110.     if (flag) return (c/d-t);
  1111.     else return t-c/d;
  1112. }
  1113.  
  1114. void minvgauss (header *hd)
  1115. {    spread1(invgauss,0,hd);
  1116. }
  1117.  
  1118. double rfak (double x)
  1119. {    int i,n;
  1120.     double res=1;
  1121.     if (x<2 || x>INT_MAX) return 1.0;
  1122.     n=(int)x;
  1123.     for (i=2; i<=n; i++) res=res*i;
  1124.     return res;
  1125. }
  1126.  
  1127. void mfak (header *hd)
  1128. {    spread1(rfak,0,hd);
  1129. }
  1130.  
  1131. void rbin (double *x, double *y, double *z)
  1132. {   long i,n,m,k;
  1133.     double res;
  1134.     n=(long)*x; m=(long)*y;
  1135.     if (m<=0) *z=1.0;
  1136.     else
  1137.     {    res=k=(n-m+1);
  1138.         for (i=2; i<=m; i++) { k++; res=(res*k)/i; }
  1139.         *z=res;
  1140.     }
  1141. }
  1142.  
  1143. void mbin (header *hd)
  1144. {    spread2(rbin,0,hd);
  1145. }
  1146.  
  1147. void rtd (double *xa, double *yf, double *zres)
  1148. {    double t,t1,a,b,h,z,p,y,x;
  1149.     int flag=0;
  1150.     if (fabs(*xa)<epsilon) { *zres=0.5; return; }
  1151.     if (*xa<0) flag=1;
  1152.     t=*xa * *xa;
  1153.     if (t>=1) { a=1; b=*yf; t1=t; }
  1154.     else { a=*yf; b=1; t1=1/t; }
  1155.     y=2/(9*a); z=2/(9*b);
  1156.     h=pow(t1,1.0/3);
  1157.     x=fabs((1-z)*h-1+y)/sqrt(z*h*h+y);
  1158.     if (b<4) x=x*(1+0.08*x*x*x*x/(b*b*b));
  1159.     h=1+x*(0.196854+x*(0.115194+x*(0.000344+x*0.019527)));
  1160.     p=0.5/(h*h*h*h);
  1161.     if (t<0.5) *zres=p/2+0.5;
  1162.     else *zres=1-p/2;
  1163.     if (flag) *zres=1-*zres;
  1164. }
  1165.  
  1166. void mtd (header *hd)
  1167. {    spread2(rtd,0,hd);
  1168. }
  1169.  
  1170. void invrtd (double *x, double *y, double *zres)
  1171. {    double z=*x,f=*y,g1,g2,g3,g4,z2;
  1172.     int flag=0;
  1173.     if (z<0.5) { flag=1; z=1-z; }
  1174.     z=invgauss(z);
  1175.     z2=z*z;
  1176.     g1=z*(1+z2)/4.0;
  1177.     g2=z*(3+z2*(16+5*z2))/96.0;
  1178.     g3=z*(-15+z2*(17+z2*(19+z2*3)))/384.0;
  1179.     g4=z*(-945+z2*(-1920+z2*(1482+z2*(776+z2*79))))/92160.0;
  1180.     *zres=(((g4/f+g3)/f+g2)/f+g1)/f+z;
  1181.     if (flag) *zres=-*zres;
  1182. }
  1183.  
  1184. void minvtd (header *hd)
  1185. {    spread2(invrtd,0,hd);
  1186. }
  1187.  
  1188. void chi (double *xa, double *yf, double *zres)
  1189. {    double ch=*xa,x,y,s,t,g,j=1;
  1190.     long i=1,p,f;
  1191.     f=(long)*yf;
  1192.     if (ch<epsilon) { *zres=0.0; return; }
  1193.     p=(f+1)/2;
  1194.     for (i=f; i>=2; i-=2) j=j*i;
  1195.     x=pow(ch,p)*exp(-(ch/2))/j;
  1196.     if (f%2==0) y=1;
  1197.     else y=sqrt(2/(ch*M_PI));
  1198.     s=1; t=1; g=f;
  1199.     while (t>1e-9)
  1200.     {    g=g+2;
  1201.         t=t*ch/g;
  1202.         s=s+t;
  1203.     }
  1204.     *zres=x*y*s;
  1205. }
  1206.  
  1207. void mchi (header *hd)
  1208. {    spread2(chi,0,hd);
  1209. }
  1210.  
  1211. double f1,f2;
  1212.  
  1213. double rfd (double F)
  1214. {    double f0,a,b,h,z,p,y,x;
  1215.     if (F<epsilon) return 0.0;
  1216.     if (F<1) { a=f2; b=f1; f0=1/F; }
  1217.     else { a=f1; b=f2; f0=F; }
  1218.     y=2/(9*a); z=2/(9*b);
  1219.     h=pow(f0,1.0/3);
  1220.     x=fabs((1-z)*h-1+y)/sqrt(z*h*h+y);
  1221.     if (b<4) x=x*(1+0.08*x*x*x*x/(b*b*b));
  1222.     h=1+x*(0.196854+x*(0.115194+x*(0.000344+x*0.019527)));
  1223.     p=0.5/(h*h*h*h);
  1224.     if (F>=1) return 1-p;
  1225.     else return p;
  1226. }
  1227.  
  1228. void mfdis (header *hd)
  1229. {    header *st=hd,*hd1,*hd2=0;
  1230.     hd1=next_param(st);
  1231.     if (hd1)
  1232.     {    hd2=next_param(hd1);
  1233.         hd1=getvalue(hd1);
  1234.     }
  1235.     if (hd2) hd2=getvalue(hd2);
  1236.     if (error) return;
  1237.     if (hd1->type!=s_real || hd2->type!=s_real) wrong_arg();
  1238.     f1=*realof(hd1);
  1239.     f2=*realof(hd2);
  1240.     spread1(rfd,0,hd);
  1241. }
  1242.  
  1243. void rmax (double *x, double *y, double *z)
  1244. {    if (*x>*y) *z=*x;
  1245.     else *z=*y;
  1246. }
  1247.  
  1248. void mmax (header *hd)
  1249. {    spread2(rmax,0,hd);
  1250. }
  1251.  
  1252. void rmin (double *x, double *y, double *z)
  1253. {    if (*x>*y) *z=*y;
  1254.     else *z=*x;
  1255. }
  1256.  
  1257. void mmin (header *hd)
  1258. {    spread2(rmin,0,hd);
  1259. }
  1260.  
  1261. typedef struct { double val; int ind; } sorttyp;
  1262.  
  1263. int sorttyp_compare (const sorttyp *x, const sorttyp *y)
  1264. {    if (x->val>y->val) return 1;
  1265.     else if (x->val==y->val) return 0;
  1266.     else return -1;
  1267. }
  1268.  
  1269. void msort (header *hd)
  1270. {    header *st=hd,*result,*result1;
  1271.     double *m,*m1;
  1272.     sorttyp *t;
  1273.     int r,c,i;
  1274.     hd=getvalue(hd); if (error) return;
  1275.     if (hd->type!=s_real && hd->type!=s_matrix) wrong_arg();
  1276.     getmatrix(hd,&r,&c,&m);
  1277.     if (c==1 || r==1) result=new_matrix(r,c,"");
  1278.     else wrong_arg();
  1279.     if (error) return;
  1280.     result1=new_matrix(r,c,"");
  1281.     if (error) return;
  1282.     if (c==1) c=r;
  1283.     if (c==0) wrong_arg();
  1284.     if (newram+c*sizeof(sorttyp)>ramend)
  1285.     {    output("Out of memory!\n"); error=600; return; 
  1286.     }
  1287.     t=(sorttyp *)newram;
  1288.     for (i=0; i<c; i++)
  1289.     {    t->val=*m++; t->ind=i; t++;
  1290.     }
  1291.     qsort(newram,c,sizeof(sorttyp),
  1292.         (int (*) (const void *, const void *))sorttyp_compare);
  1293.     m=matrixof(result); m1=matrixof(result1);
  1294.     t=(sorttyp *)newram;
  1295.     for (i=0; i<c; i++)
  1296.     {    *m++=t->val; *m1++=t->ind+1; t++;
  1297.     }
  1298.     moveresult(st,result);
  1299.     moveresult(nextof(st),result1);
  1300. }
  1301.  
  1302. void mnonzeros (header *hd)
  1303. {    header *st=hd,*result;
  1304.     double *m,*mr;
  1305.     int r,c,i,k;
  1306.     hd=getvalue(hd); if (error) return;
  1307.     if (hd->type!=s_real && hd->type!=s_matrix) wrong_arg();
  1308.     getmatrix(hd,&r,&c,&m);
  1309.     if (r!=1 && c!=1) wrong_arg();
  1310.     if (c==1) c=r;
  1311.     result=new_matrix(1,c,""); if (error) return;
  1312.     k=0; mr=matrixof(result);
  1313.     for (i=0; i<c; i++)
  1314.     {    if (*m++!=0.0)
  1315.         {    *mr++=i+1; k++;
  1316.         }
  1317.     }
  1318.     dimsof(result)->c=k;
  1319.     result->size=matrixsize(1,k);
  1320.     moveresult(st,result);
  1321. }
  1322.  
  1323. void mstatistics (header *hd)
  1324. {    header *st=hd,*hd1,*result;
  1325.     int i,n,r,c,k;
  1326.     double *m,*mr;
  1327.     hd=getvalue(hd);
  1328.     hd1=next_param(st);
  1329.     if (hd1) hd1=getvalue(hd1); if (error) return;
  1330.     if (hd1->type!=s_real || hd->type!=s_matrix) wrong_arg();
  1331.     if (*realof(hd1)>INT_MAX || *realof(hd1)<2) wrong_arg();
  1332.     n=(int)*realof(hd1);
  1333.     getmatrix(hd,&r,&c,&m);
  1334.     if (r!=1 && c!=1) wrong_arg();
  1335.     if (c==1) c=r;
  1336.     result=new_matrix(1,n,""); if (error) return;
  1337.     mr=matrixof(result); for (i=0; i<n; i++) *mr++=0.0;
  1338.     mr=matrixof(result);
  1339.     for (i=0; i<c; i++)
  1340.     {    if (*m>=0 && *m<n)
  1341.         {    k=floor(*m);
  1342.             mr[k]+=1.0;
  1343.         }
  1344.         m++;
  1345.     }
  1346.     moveresult(st,result);
  1347. }
  1348.  
  1349. void minput (header *hd)
  1350. {    header *st=hd,*result;
  1351.     char input[1024],*oldnext;
  1352.     hd=getvalue(hd); if (error) return;
  1353.     if (hd->type!=s_string) wrong_arg();
  1354.     retry: output(stringof(hd)); output("? ");
  1355.     edit(input);
  1356.     stringon=1;
  1357.     oldnext=next; next=input; result=scan_value(); next=oldnext;
  1358.     stringon=0;
  1359.     if (error) 
  1360.     {    output("Error in input!\n"); error=0; goto retry;
  1361.     }
  1362.     moveresult(st,result);
  1363. }
  1364.  
  1365. void mlineinput (header *hd)
  1366. {    header *st=hd,*result;
  1367.     char input[1024];
  1368.     hd=getvalue(hd); if (error) return;
  1369.     if (hd->type!=s_string) wrong_arg();
  1370.     output(stringof(hd)); output("? ");
  1371.     edit(input);
  1372.     result=new_string(input,strlen(input),"");
  1373.     moveresult(st,result);
  1374. }
  1375.  
  1376. void minterpret (header *hd)
  1377. {    header *st=hd,*result;
  1378.     char *oldnext;
  1379.     hd=getvalue(hd); if (error) return;
  1380.     if (hd->type!=s_string) wrong_arg();
  1381.     stringon=1;
  1382.     oldnext=next; next=stringof(hd); result=scan(); next=oldnext;
  1383.     stringon=0;
  1384.     if (error) { result=new_string("Syntax error!",5,""); error=0; }
  1385.     moveresult(st,result);
  1386. }
  1387.  
  1388. void mmax1 (header *hd)
  1389. {    header *result,*st=hd;
  1390.     double x,*m,*mr,max;
  1391.     int r,c,i,j;
  1392.     hd=getvalue(hd); if (error) return;
  1393.     if (hd->type==s_real || hd->type==s_matrix)
  1394.     {    getmatrix(hd,&r,&c,&m);
  1395.         result=new_matrix(r,1,""); if (error) return;
  1396.         mr=matrixof(result);
  1397.         for (i=0; i<r; i++) 
  1398.         {    max=*m; m++;
  1399.             for (j=1; j<c; j++) 
  1400.             {    x=*m++;
  1401.                 if (x>max) max=x;
  1402.             }
  1403.             *mr++=max;
  1404.         }
  1405.     }
  1406.     else wrong_arg();
  1407.     moveresult(st,result);
  1408. }
  1409.  
  1410. void mmin1 (header *hd)
  1411. {    header *result,*st=hd;
  1412.     double x,*m,*mr,max;
  1413.     int r,c,i,j;
  1414.     hd=getvalue(hd); if (error) return;
  1415.     if (hd->type==s_real || hd->type==s_matrix)
  1416.     {    getmatrix(hd,&r,&c,&m);
  1417.         result=new_matrix(r,1,""); if (error) return;
  1418.         mr=matrixof(result);
  1419.         for (i=0; i<r; i++) 
  1420.         {    max=*m; m++;
  1421.             for (j=1; j<c; j++) 
  1422.             {    x=*m++;
  1423.                 if (x<max) max=x;
  1424.             }
  1425.             *mr++=max;
  1426.         }
  1427.     }
  1428.     else wrong_arg();
  1429.     moveresult(st,result);
  1430. }
  1431.  
  1432. void make_xors (void)
  1433. {    int i;
  1434.     for (i=0; i<10; i++) xors[i]=xor(argname[i]);
  1435. }
  1436.  
  1437. void mdo (header *hd)
  1438. {    header *st=hd,*hd1,*result;
  1439.     int count=0;
  1440.     size_t size;
  1441.     if (!hd) wrong_arg();
  1442.     hd=getvalue(hd);
  1443.     result=hd1=next_param(st);
  1444.     if (hd->type!=s_string) wrong_arg();
  1445.     if (error) return;
  1446.     hd=searchudf(stringof(hd));
  1447.     if (!hd || hd->type!=s_udf) wrong_arg();
  1448.     while (hd1) 
  1449.     {    strcpy(hd1->name,argname[count]);
  1450.         hd1->xor=xors[count];
  1451.         hd1=next_param(hd1); count++; 
  1452.     }
  1453.     if (result)
  1454.     {    size=(char *)result-(char *)st;
  1455.         if (size>0 && newram!=(char *)result) 
  1456.             memmove((char *)st,(char *)result,newram-(char *)result);
  1457.         newram-=size;
  1458.     }
  1459.     interpret_udf(hd,st,count);
  1460. }
  1461.  
  1462. void mlu (header *hd)
  1463. {    header *st=hd,*result,*res1,*res2,*res3;
  1464.     double *m,*mr,*m1,*m2,det,deti;
  1465.     int r,c,*rows,*cols,rank,i;
  1466.     hd=getvalue(hd); if (error) return;
  1467.     if (hd->type==s_matrix || hd->type==s_real)
  1468.     {    getmatrix(hd,&r,&c,&m);
  1469.         if (r<1) wrong_arg();
  1470.         result=new_matrix(r,c,""); if (error) return;
  1471.         mr=matrixof(result);
  1472.         memmove((char *)mr,(char *)m,(LONG)r*c*sizeof(double));
  1473.         make_lu(mr,r,c,&rows,&cols,&rank,&det); if (error) return;
  1474.         res1=new_matrix(1,rank,""); if (error) return;
  1475.         res2=new_matrix(1,c,""); if (error) return;
  1476.         res3=new_real(det,""); if (error) return;
  1477.         m1=matrixof(res1);
  1478.         for (i=0; i<rank; i++)
  1479.         {    *m1++=*rows+1;
  1480.             rows++;
  1481.         }
  1482.         m2=matrixof(res2);
  1483.         for (i=0; i<c; i++)
  1484.         {    *m2++=*cols++;
  1485.         }
  1486.         moveresult(st,getvalue(result)); st=nextof(st);
  1487.         moveresult(st,getvalue(res1)); st=nextof(st);
  1488.         moveresult(st,getvalue(res2)); st=nextof(st);
  1489.         moveresult(st,getvalue(res3));
  1490.     }
  1491.     else if (hd->type==s_cmatrix || hd->type==s_complex)
  1492.     {    getmatrix(hd,&r,&c,&m);
  1493.         if (r<1) wrong_arg();
  1494.         result=new_cmatrix(r,c,""); if (error) return;
  1495.         mr=matrixof(result);
  1496.         memmove((char *)mr,(char *)m,(LONG)r*c*(LONG)2*sizeof(double));
  1497.         cmake_lu(mr,r,c,&rows,&cols,&rank,&det,&deti); 
  1498.             if (error) return;
  1499.         res1=new_matrix(1,rank,""); if (error) return;
  1500.         res2=new_matrix(1,c,""); if (error) return;
  1501.         res3=new_complex(det,deti,""); if (error) return;
  1502.         m1=matrixof(res1);
  1503.         for (i=0; i<rank; i++)
  1504.         {    *m1++=*rows+1;
  1505.             rows++;
  1506.         }
  1507.         m2=matrixof(res2);
  1508.         for (i=0; i<c; i++)
  1509.         {    *m2++=*cols++;
  1510.         }
  1511.         moveresult(st,getvalue(result)); st=nextof(st);
  1512.         moveresult(st,getvalue(res1)); st=nextof(st);
  1513.         moveresult(st,getvalue(res2)); st=nextof(st);
  1514.         moveresult(st,getvalue(res3));
  1515.     }
  1516.     else wrong_arg();
  1517. }
  1518.  
  1519. void miscomplex (header *hd)
  1520. {    header *st=hd,*result;
  1521.     hd=getvalue(hd);
  1522.     if (hd->type==s_complex || hd->type==s_cmatrix)
  1523.         result=new_real(1.0,"");
  1524.     else result=new_real(0.0,"");
  1525.     if (error) return;
  1526.     moveresult(st,result);
  1527. }
  1528.  
  1529. void misreal (header *hd)
  1530. {    header *st=hd,*result;
  1531.     hd=getvalue(hd); if (error) return;
  1532.     if (hd->type==s_real || hd->type==s_matrix)
  1533.         result=new_real(1.0,"");
  1534.     else result=new_real(0.0,"");
  1535.     if (error) return;
  1536.     moveresult(st,result);
  1537. }
  1538.  
  1539. double rounder;
  1540.  
  1541. double rround (double x)
  1542. {    x*=rounder;
  1543.     if (x>0) x=floor(x+0.5);
  1544.     else x=-floor(-x+0.5);
  1545.     return x/rounder;
  1546. }
  1547.  
  1548. void cround (double *x, double *xi, double *z, double *zi)
  1549. {    *z=rround(*x);
  1550.     *zi=rround(*xi);
  1551. }
  1552.  
  1553. double frounder[]={1.0,10.0,100.0,1000.0,10000.0,100000.0,1000000.0,
  1554. 10000000.0,100000000.0,1000000000.0,10000000000.0};
  1555.  
  1556. void mround (header *hd)
  1557. {    header *hd1;
  1558.     int n;
  1559.     hd1=next_param(hd);
  1560.     if (hd1) hd1=getvalue(hd1); if (error) return;
  1561.     if (hd1->type!=s_real) wrong_arg();
  1562.     n=(int)(*realof(hd1));
  1563.     if (n>0 && n<11) rounder=frounder[n];
  1564.     else rounder=pow(10.0,n);
  1565.     spread1(rround,cround,hd);
  1566. }
  1567.  
  1568. void mchar (header *hd)
  1569. {    header *st=hd,*result;
  1570.     hd=getvalue(hd); if (error) return;
  1571.     if (hd->type!=s_real) wrong_arg();
  1572.     result=new_string("a",1,""); if (error) return;
  1573.     *stringof(result)=(char)*realof(hd);
  1574.     moveresult(st,result);
  1575. }
  1576.  
  1577. void merror (header *hd)
  1578. {    hd=getvalue(hd); if (error) return;
  1579.     if (hd->type!=s_string) wrong_arg();
  1580.     output1("Error : %s\n",stringof(hd));
  1581.     error=301;
  1582. }
  1583.  
  1584. void merrlevel (header *hd)
  1585. {    header *st=hd,*res;
  1586.     char *oldnext;
  1587.     int en;
  1588.     hd=getvalue(hd); if (error) return;
  1589.     if (hd->type!=s_string) wrong_arg();
  1590.     stringon=1;
  1591.     oldnext=next; next=stringof(hd); scan(); next=oldnext;
  1592.     stringon=0;
  1593.     en=error; error=0;
  1594.     res=new_real(en,""); if (error) return;
  1595.     moveresult(st,res);
  1596. }
  1597.  
  1598. void mprintf (header *hd)
  1599. {    header *st=hd,*hd1,*result;
  1600.     char string[1024];
  1601.     hd1=next_param(hd);
  1602.     hd=getvalue(hd);
  1603.     hd1=getvalue(hd1); if (error) return;
  1604.     if (hd->type!=s_string || hd1->type!=s_real)
  1605.         wrong_arg();
  1606.     sprintf(string,stringof(hd),*realof(hd1));
  1607.     result=new_string(string,strlen(string),""); if (error) return;
  1608.     moveresult(st,result);
  1609. }
  1610.  
  1611. void msetkey (header *hd)
  1612. /*****
  1613.     set a function key
  1614. *****/
  1615. {    header *st=hd,*hd1,*result;
  1616.     char *p;
  1617.     int n;
  1618.     hd=getvalue(hd); if (error) return;
  1619.     hd1=nextof(st); hd1=getvalue(hd1); if (error) return;
  1620.     if (hd->type!=s_real || hd1->type!=s_string) wrong_arg();
  1621.     n=(int)(*realof(hd))-1; p=stringof(hd1);
  1622.     if (n<0 || n>=10 || strlen(p)>63) wrong_arg();
  1623.     result=new_string(fktext[n],strlen(fktext[n]),"");
  1624.     if (error) return;
  1625.     strcpy(fktext[n],p);
  1626.     moveresult(st,result);
  1627. }
  1628.  
  1629. void many (header *hd)
  1630. {    header *st=hd,*result;
  1631.     int c,r,res=0;
  1632.     LONG i,n;
  1633.     double *m;
  1634.     hd=getvalue(hd); if (error) return;
  1635.     if (hd->type==s_real || hd->type==s_matrix)
  1636.     {    getmatrix(hd,&r,&c,&m);
  1637.         n=(LONG)(c)*r;
  1638.     }
  1639.     else if (hd->type==s_complex || hd->type==s_cmatrix)
  1640.     {    getmatrix(hd,&r,&c,&m);
  1641.         n=(LONG)2*(LONG)(c)*r;
  1642.     }
  1643.     else wrong_arg();
  1644.     for (i=0; i<n; i++)
  1645.         if (*m++!=0.0) { res=1; break; }
  1646.     result=new_real(res,""); if (error) return;
  1647.     moveresult(st,result);
  1648. }
  1649.  
  1650. void mcd (header *hd)
  1651. {    header *st=hd,*result;
  1652.     char *path;
  1653.     hd=getvalue(hd); if (error) return;
  1654.     if (hd->type!=s_string) wrong_arg();
  1655.     path=cd(stringof(hd));
  1656.     result=new_string(path,strlen(path),"");
  1657.     moveresult(st,result);
  1658. }
  1659.  
  1660. void mdir (header *hd)
  1661. {    header *st=hd,*result;
  1662.     char *name;
  1663.     hd=getvalue(hd); if (error) return;
  1664.     if (hd->type!=s_string) wrong_arg();
  1665.     name=dir(stringof(hd));
  1666.     if (name) result=new_string(name,strlen(name),"");
  1667.     else result=new_string("",0,"");
  1668.     if (error) return;
  1669.     moveresult(st,result);
  1670. }
  1671.  
  1672. void margs (header *hd)
  1673. /* return all args from realof(hd)-st argument on */
  1674. {    header *st=hd,*hd1,*result;
  1675.     int i,n;
  1676.     size_t size;
  1677.     hd=getvalue(hd);
  1678.     if (hd->type!=s_real) wrong_arg();
  1679.     n=(int)*realof(hd);
  1680.     if (n<1) wrong_arg();
  1681.     if (n>actargn)
  1682.     {    newram=(char *)st; return;
  1683.     }
  1684.     result=(header *)startlocal; i=1;
  1685.     while (i<n && result<(header *)endlocal)
  1686.     {    result=nextof(result); i++;
  1687.     }
  1688.     hd1=result;
  1689.     while (i<actargn+1 && hd1<(header *)endlocal)
  1690.     {    hd1=nextof(hd1); i++;
  1691.     }
  1692.     size=(char *)hd1-(char *)result;
  1693.     if (size<=0)
  1694.     {    output("Error in args!\n"); error=2021; return;
  1695.     }
  1696.     memmove((char *)st,(char *)result,size);
  1697.     newram=(char *)st+size;
  1698. }
  1699.  
  1700. void mname (header *hd)
  1701. {    header *st=hd,*result;
  1702.     hd=getvalue(hd); if (error) return;
  1703.     result=new_string(hd->name,strlen(hd->name),"");
  1704.     moveresult(st,result);
  1705. }
  1706.  
  1707. void mdir0 (header *hd)
  1708. {    char *name;
  1709.     name=dir(0);
  1710.     if (name) new_string(name,strlen(name),"");
  1711.     else new_string("",0,"");
  1712. }
  1713.  
  1714. void mflipx (header *hd)
  1715. {    header *st=hd,*result;
  1716.     double *m,*mr,*mr1;
  1717.     int i,j,c,r;
  1718.     hd=getvalue(hd); if (error) return;
  1719.     if (hd->type==s_real || hd->type==s_complex)
  1720.     {    moveresult(st,hd); return;
  1721.     }
  1722.     else if (hd->type==s_matrix)
  1723.     {    getmatrix(hd,&r,&c,&m);
  1724.         result=new_matrix(r,c,""); if (error) return;
  1725.         mr=matrixof(result);
  1726.         for (i=0; i<r; i++)
  1727.         {    mr1=mr+(c-1);
  1728.             for (j=0; j<c; j++) *mr1--=*m++;
  1729.             mr+=c;
  1730.         }
  1731.     }
  1732.     else if (hd->type==s_cmatrix)
  1733.     {    getmatrix(hd,&r,&c,&m);
  1734.         result=new_cmatrix(r,c,""); if (error) return;
  1735.         mr=matrixof(result);
  1736.         for (i=0; i<r; i++)
  1737.         {    mr1=mr+(2l*(c-1)+1);
  1738.             for (j=0; j<c; j++)
  1739.             {    *mr1--=*m++; *mr1--=*m++;
  1740.             }
  1741.             mr+=2l*c;
  1742.         }
  1743.     }
  1744.     else wrong_arg();
  1745.     moveresult(st,result);
  1746. }
  1747.  
  1748. void mflipy (header *hd)
  1749. {    header *st=hd,*result;
  1750.     double *m,*mr;
  1751.     int i,c,r;
  1752.     hd=getvalue(hd); if (error) return;
  1753.     if (hd->type==s_real || hd->type==s_complex)
  1754.     {    moveresult(st,hd); return;
  1755.     }
  1756.     else if (hd->type==s_matrix)
  1757.     {    getmatrix(hd,&r,&c,&m);
  1758.         result=new_matrix(r,c,""); if (error) return;
  1759.         mr=matrixof(result);
  1760.         mr+=(long)(r-1)*c;
  1761.         for (i=0; i<r; i++)
  1762.         {    memmove((char *)mr,(char *)m,c*sizeof(double));
  1763.             m+=c; mr-=c;
  1764.         }
  1765.     }
  1766.     else if (hd->type==s_cmatrix)
  1767.     {    getmatrix(hd,&r,&c,&m);
  1768.         result=new_cmatrix(r,c,""); if (error) return;
  1769.         mr=matrixof(result);
  1770.         mr+=2l*(long)(r-1)*c;
  1771.         for (i=0; i<r; i++)
  1772.         {    memmove((char *)mr,(char *)m,2l*c*sizeof(double));
  1773.             m+=2l*c; mr-=2l*c;
  1774.         }
  1775.     }
  1776.     else wrong_arg();
  1777.     moveresult(st,result);
  1778. }
  1779.  
  1780. void mmatrix (header *hd)
  1781. {    header *st=hd,*hd1,*result;
  1782.     long i,n;
  1783.     double x,xi;
  1784.     double *m,*mr;
  1785.     int c,r,c1,r1;
  1786.     hd1=nextof(hd);
  1787.     hd=getvalue(hd);
  1788.     if (error) return;
  1789.     hd1=getvalue(hd1);
  1790.     if (error) return;
  1791.     if (hd->type==s_matrix)
  1792.     {    getmatrix(hd,&r,&c,&m);
  1793.         if (*m<0 || *m>INT_MAX || *(m+1)<0 || *(m+1)>INT_MAX)
  1794.             wrong_arg();
  1795.         r1=(int)*m; c1=(int)*(m+1);
  1796.         if (hd1->type==s_real)
  1797.         {    result=new_matrix(r1,c1,"");
  1798.             mr=matrixof(result);
  1799.             x=*realof(hd1);
  1800.             n=(long)c1*r1;
  1801.             for (i=0; i<n; i++) *mr++=x;
  1802.         }
  1803.         else if (hd1->type==s_complex)
  1804.         {    result=new_cmatrix(r1,c1,"");
  1805.             mr=matrixof(result);
  1806.             x=*realof(hd1); xi=*(realof(hd1)+1);
  1807.             n=(long)c1*r1;
  1808.             for (i=0; i<n; i++) 
  1809.             {    *mr++=x; *mr++=xi;
  1810.             }
  1811.         }
  1812.         else wrong_arg();
  1813.     }
  1814.     else wrong_arg();
  1815.     moveresult(st,result);
  1816. }
  1817.  
  1818. /*************** execute builtin functions ***********/
  1819.  
  1820. int builtin_count;
  1821.  
  1822. extern builtintyp builtin_list[];
  1823.  
  1824. void print_builtin (void)
  1825. {    int linel=0,i;
  1826.     for (i=0; i<builtin_count; i++)
  1827.     {    if (linel+strlen(builtin_list[i].name)+2>linelength)
  1828.             { output("\n"); linel=0; }
  1829.         output1("%s ",builtin_list[i].name);
  1830.         linel+=(int)strlen(builtin_list[i].name)+1;
  1831.     }
  1832.     output("\n");
  1833. }
  1834.  
  1835. int builtin_compare (const builtintyp *p1, const builtintyp *p2)
  1836. {    int h;
  1837.     h=strcmp(p1->name,p2->name);
  1838.     if (h) return h;
  1839.     else 
  1840.     {    if (p1->nargs==-1 || p2->nargs==-1) return 0;
  1841.         else if (p1->nargs<p2->nargs) return -1; 
  1842.         else if (p1->nargs>p2->nargs) return 1;
  1843.         else return 0;
  1844.     }
  1845. }
  1846.  
  1847. void sort_builtin (void)
  1848. {    builtin_count=0;
  1849.     while (builtin_list[builtin_count].name) builtin_count++;
  1850.     qsort(builtin_list,builtin_count,sizeof(builtintyp),
  1851.         (int (*) (const void *, const void *))builtin_compare);
  1852. }
  1853.  
  1854. int exec_builtin (char *name, int nargs, header *hd)
  1855. {    builtintyp *b=builtin_list,h;
  1856.     h.name=name; h.nargs=nargs;
  1857.     b=bsearch(&h,builtin_list,builtin_count,sizeof(builtintyp),
  1858.         (int (*) (const void *, const void *))builtin_compare);
  1859.     if (b)
  1860.     {    if (nargs==0) hd=0;
  1861.         b->f(hd); return 1;
  1862.     }
  1863.     else return 0;
  1864. }
  1865.  
  1866. #ifndef SPLIT_MEM
  1867.  
  1868. typedef struct { size_t udfend,startlocal,endlocal,newram; }
  1869.     ptyp;
  1870.  
  1871. void mstore (header *hd)
  1872. {    FILE *file;
  1873.     ptyp p;
  1874.     hd=getvalue(hd); if (error) return;
  1875.     if (hd->type!=s_string)
  1876.     {    output("Expect file name.\n");
  1877.         error=1100; return;
  1878.     }
  1879.     p.udfend=udfend-ramstart;
  1880.     p.startlocal=startlocal-ramstart;
  1881.     p.endlocal=endlocal-ramstart;
  1882.     p.newram=newram-ramstart;
  1883.     file=fopen(stringof(hd),"wb");
  1884.     if (!file)
  1885.     {    output1("Could not open %s.\n",stringof(hd));
  1886.         error=1101; return;
  1887.     }
  1888.     fwrite(&p,sizeof(ptyp),1,file);
  1889.     fwrite(ramstart,1,newram-ramstart,file);
  1890.     if (ferror(file))
  1891.     {    output("Write error.\n");
  1892.         error=1102; return;
  1893.     }
  1894.     fclose(file);
  1895. }
  1896.     
  1897. void mrestore (header *hd)
  1898. {    FILE *file;
  1899.     ptyp p;
  1900.     hd=getvalue(hd); if (error) return;
  1901.     if (hd->type!=s_string)
  1902.     {    output("Expect file name.\n");
  1903.         error=1100; return;
  1904.     }
  1905.     file=fopen(stringof(hd),"rb");
  1906.     if (!file)
  1907.     {    output1("Could not open %s.\n",stringof(hd));
  1908.         error=1103; return;
  1909.     }
  1910.     fread(&p,sizeof(ptyp),1,file);
  1911.     if (ferror(file))
  1912.     {    output("Read error.\n");
  1913.         error=1104; return;
  1914.     }
  1915.     fread(ramstart,1,p.newram,file);
  1916.     if (ferror(file))
  1917.     {    output("Read error (fatal for EULER).\n");
  1918.         error=1104; return;
  1919.     }
  1920.     fclose(file);
  1921.     udfend=ramstart+p.udfend;
  1922.     startlocal=ramstart+p.startlocal;
  1923.     endlocal=ramstart+p.endlocal;
  1924.     newram=ramstart+p.newram;
  1925. }
  1926.     
  1927. #endifə