home *** CD-ROM | disk | FTP | other *** search
/ Crawly Crypt Collection 1 / crawlyvol1.bin / apps / math / euler / source / express.c < prev    next >
C/C++ Source or Header  |  1993-03-25  |  34KB  |  1,438 lines
  1. #include <math.h>
  2. #include <stdio.h>
  3. #include <ctype.h>
  4. #include <string.h>
  5. #include <limits.h>
  6.  
  7. #include "header.h"
  8. #include "sysdep.h"
  9. #include "funcs.h"
  10.  
  11. void scan_summe (void);
  12.  
  13. extern int nosubmref,trace;
  14. int udfon=0,actargn=0;
  15. header *running;
  16.  
  17. /**************** builtin operators ****************/
  18.  
  19. header *getvalue (header *hd)
  20. /***** getvalue
  21.     get an actual value of a reference.
  22. *****/
  23. {    header *old=hd,*mhd,*result;
  24.     dims *d;
  25.     double *m,*mr,*m1,*m2,*m3;
  26.     int r,c,*rind,*cind,*cind1,i,j;
  27.     while (hd && hd->type==s_reference)
  28.         hd=referenceof(hd);
  29.     if (!hd)
  30.     {    mhd=(header *)newram;
  31.         if (exec_builtin(old->name,0,mhd)) return mhd;
  32.         hd=searchudf(old->name);
  33.         if (hd)
  34.         {    interpret_udf(hd,mhd,0);
  35.             return mhd;
  36.         }
  37.         output1("Variable %s not defined!\n",old->name);
  38.         error=10; return new_string("Fehler",6,"");
  39.     }
  40.     if (hd->type==s_submatrix)
  41.     {    mhd=submrefof(hd); d=submdimsof(hd);
  42.         rind=rowsof(hd); cind=colsof(hd);
  43.         getmatrix(mhd,&r,&c,&m);
  44.         if (d->r==1 && d->c==1)
  45.             return new_real(*mat(m,c,*rind,*cind),"");
  46.         result=new_matrix(d->r,d->c,"");
  47.         mr=matrixof(result);
  48.         for (i=0; i<d->r; i++)
  49.         {    cind1=cind;
  50.             m1=mat(mr,d->c,i,0);
  51.             m2=mat(m,c,*rind,0);
  52.             for (j=0; j<d->c; j++)
  53.             {    m1[j]=m2[*cind1];
  54.                 cind1++;
  55.             }
  56.             rind++;
  57.         }
  58.         return result;
  59.     }
  60.     if (hd->type==s_csubmatrix)
  61.     {    mhd=submrefof(hd); d=submdimsof(hd);
  62.         rind=rowsof(hd); cind=colsof(hd);
  63.         getmatrix(mhd,&r,&c,&m);
  64.         if (d->r==1 && d->c==1)
  65.         {    m=cmat(m,c,*rind,*cind);
  66.             return new_complex(*m,*(m+1),"");
  67.         }
  68.         result=new_cmatrix(d->r,d->c,"");
  69.         mr=matrixof(result);
  70.         for (i=0; i<d->r; i++)
  71.         {    cind1=cind;
  72.             m1=cmat(mr,d->c,i,0);
  73.             m2=cmat(m,c,*rind,0);
  74.             for (j=0; j<d->c; j++)
  75.             {   m3=m2+(LONG)2*(*cind1);
  76.                 *m1++=*m3++; *m1++=*m3;
  77.                 cind1++;
  78.             }
  79.             rind++;
  80.         }
  81.         return result;
  82.     }
  83.     if (hd->type==s_matrix && dimsof(hd)->c==1 && dimsof(hd)->r==1)
  84.     {    return new_real(*matrixof(hd),"");
  85.     }
  86.     if (hd->type==s_cmatrix && dimsof(hd)->c==1 && dimsof(hd)->r==1)
  87.     {    return new_complex(*matrixof(hd),*(matrixof(hd)+1),"");
  88.     }
  89.     return hd;
  90. }
  91.  
  92. header *getvariable (header *hd)
  93. /***** getvariable
  94.     get an actual variable of a reference.
  95. *****/
  96. {    header *old=hd;
  97.     while (hd && hd->type==s_reference)
  98.         hd=referenceof(hd);
  99.     if (!hd)
  100.     {    output1("Variable %s not defined!\n",old->name);
  101.         error=10; return new_string("Fehler",6,"");
  102.     }
  103.     return hd;
  104. }
  105.  
  106. void moveresult (header *stack, header *result)
  107. /***** moveresult
  108.     move the result to the start of stack.
  109. *****/
  110. {    if (stack==result) return;
  111.     memmove((char *)stack,(char *)result,result->size);
  112.     newram=(char *)stack+stack->size;
  113. }
  114.  
  115. void moveresult1 (header *stack, header *result)
  116. /***** moveresult
  117.     move several results to the start of stack.
  118. *****/
  119. {    size_t size;
  120.     if (stack==result) return;
  121.     size=newram-(char *)result;
  122.     memmove((char *)stack,(char *)result,size);
  123.     newram=(char *)stack+size;
  124. }
  125.  
  126. header *mapt2 (void f(double *,double *,double *),
  127.     void fc(double *, double *, double *, double *, double *, double *),
  128.     header *hd, header *hd1)
  129. /***** map
  130.     do the function elementwise to the two values.
  131.     store the result hd2. the values may also be complex,
  132.     then fc is used.
  133. ******/
  134. {    dims *d,*d1;
  135.     double x,y,*m,*m1,*m2,null=0.0;
  136.     header *hd2;
  137.     LONG i,n;
  138.     if (hd->type==s_real)
  139.     {    if (hd1->type==s_matrix)
  140.         {    d1=dimsof(hd1);
  141.             m1=matrixof(hd1);
  142.             x=*(realof(hd));
  143.             hd2=new_matrix(d1->r,d1->c,"");
  144.             if (error) return new_string("Fehler",6,"");
  145.             m2=matrixof(hd2);
  146.             n=(d1->c)*(d1->r);
  147.             for (i=0; i<n; i++)
  148.                 {    f(m1,&x,m2); m1++; m2++;
  149.                 }
  150.             return hd2;
  151.         }
  152.         if (fc && hd1->type==s_cmatrix)
  153.         {    d1=dimsof(hd1);
  154.             m1=matrixof(hd1);
  155.             x=*(realof(hd));
  156.             hd2=new_cmatrix(d1->r,d1->c,"");
  157.             if (error) return new_string("Fehler",6,"");
  158.             m2=matrixof(hd2);
  159.             n=(d1->c)*(d1->r);
  160.             for (i=0; i<n; i++)
  161.                 {    fc(m1,m1+1,&x,&null,m2,m2+1); m1+=2; m2+=2;
  162.                 }
  163.             return hd2;
  164.         }
  165.     }
  166.     else if (hd->type==s_complex)
  167.     {    if (hd1->type==s_matrix)
  168.         {    d1=dimsof(hd1);
  169.             m1=matrixof(hd1);
  170.             x=*(realof(hd));
  171.             y=*(imagof(hd));
  172.             hd2=new_cmatrix(d1->r,d1->c,"");
  173.             if (error) return new_string("Fehler",6,"");
  174.             m2=matrixof(hd2);
  175.             n=(d1->c)*(d1->r);
  176.             for (i=0; i<n; i++)
  177.                 {    fc(m1,&null,&x,&y,m2,m2+1); m1++; m2+=2;
  178.                 }
  179.             return hd2;
  180.         }
  181.         if (fc && hd1->type==s_cmatrix)
  182.         {    d1=dimsof(hd1);
  183.             m1=matrixof(hd1);
  184.             x=*(realof(hd));
  185.             y=*(imagof(hd));
  186.             hd2=new_cmatrix(d1->r,d1->c,"");
  187.             if (error) return new_string("Fehler",6,"");
  188.             m2=matrixof(hd2);
  189.             n=(d1->c)*(d1->r);
  190.             for (i=0; i<n; i++)
  191.                 {    fc(m1,m1+1,&x,&y,m2,m2+1); m1+=2; m2+=2;
  192.                 }
  193.             return hd2;
  194.         }
  195.     }
  196.     else if (hd->type==s_matrix)
  197.     {    if (fc && hd1->type==s_cmatrix)
  198.         {    d=dimsof(hd);
  199.             m=matrixof(hd);
  200.             d1=dimsof(hd1);
  201.             m1=matrixof(hd1);
  202.             if (d->c!=d1->c || d->r!= d1->r)
  203.             {    output("Matrix dimensions must agree!\n");
  204.                 error=4; return new_string("Fehler",6,"");
  205.             }
  206.             hd2=new_cmatrix(d->r,d->c,"");
  207.             if (error) return new_string("Fehler",6,"");
  208.             m2=matrixof(hd2);
  209.             n=(d->c)*(d->r);
  210.             for (i=0; i<n; i++)
  211.                 {    fc(m1,m1+1,m,&null,m2,m2+1); m++; m1+=2; m2+=2;
  212.                 }
  213.             return hd2;
  214.         }
  215.     }
  216.     output("Illegal operation\n"); error=3;
  217.     return new_string("Fehler",6,"");
  218. }
  219.  
  220. header *map2 (void f(double *,double *,double *),
  221.     void fc(double *, double *, double *, double *, double *, double *),
  222.     header *hd, header *hd1)
  223. /***** map
  224.     do the function elementwise to the two values.
  225.     store the result hd2. the values may also be complex,
  226.     then fc is used.
  227. ******/
  228. {    dims *d,*d1;
  229.     double x,y,*m,*m1,*m2,null=0.0;
  230.     header *hd2;
  231.     LONG i,n;
  232.     if (hd->type==s_real)
  233.     {    if (hd1->type==s_real)
  234.         {    f(realof(hd),realof(hd1),&x);
  235.             return new_real(x,"");
  236.         }
  237.         if (fc && hd1->type==s_complex)
  238.         {    fc(realof(hd),&null,realof(hd1),imagof(hd1),&x,&y);
  239.             return new_complex(x,y,"");
  240.         }
  241.         if (hd1->type==s_matrix)
  242.         {    d1=dimsof(hd1);
  243.             m1=matrixof(hd1);
  244.             x=*(realof(hd));
  245.             hd2=new_matrix(d1->r,d1->c,"");
  246.             if (error) return new_string("Fehler",6,"");
  247.             m2=matrixof(hd2);
  248.             n=(d1->c)*(d1->r);
  249.             for (i=0; i<n; i++)
  250.                 {    f(&x,m1,m2); m1++; m2++;
  251.                 }
  252.             return hd2;
  253.         }
  254.         if (fc && hd1->type==s_cmatrix)
  255.         {    d1=dimsof(hd1);
  256.             m1=matrixof(hd1);
  257.             x=*(realof(hd));
  258.             hd2=new_cmatrix(d1->r,d1->c,"");
  259.             if (error) return new_string("Fehler",6,"");
  260.             m2=matrixof(hd2);
  261.             n=(d1->c)*(d1->r);
  262.             for (i=0; i<n; i++)
  263.                 {    fc(&x,&null,m1,m1+1,m2,m2+1); m1+=2; m2+=2;
  264.                 }
  265.             return hd2;
  266.         }
  267.     }
  268.     else if (fc && hd->type==s_complex)
  269.     {    if (hd1->type==s_real)
  270.         {    fc(realof(hd),imagof(hd),realof(hd1),&null,&x,&y);
  271.             return new_complex(x,y,"");
  272.         }
  273.         if (hd1->type==s_complex)
  274.         {    fc(realof(hd),imagof(hd),realof(hd1),imagof(hd1),&x,&y);
  275.             return new_complex(x,y,"");
  276.         }
  277.         if (hd1->type==s_matrix)
  278.         {    d1=dimsof(hd1);
  279.             m1=matrixof(hd1);
  280.             x=*(realof(hd));
  281.             y=*(imagof(hd));
  282.             hd2=new_cmatrix(d1->r,d1->c,"");
  283.             if (error) return new_string("Fehler",6,"");
  284.             m2=matrixof(hd2);
  285.             n=(d1->c)*(d1->r);
  286.             for (i=0; i<n; i++)
  287.                 {    fc(&x,&y,m1,&null,m2,m2+1); m1++; m2+=2;
  288.                 }
  289.             return hd2;
  290.         }
  291.         if (hd1->type==s_cmatrix)
  292.         {    d1=dimsof(hd1);
  293.             m1=matrixof(hd1);
  294.             x=*(realof(hd));
  295.             y=*(imagof(hd));
  296.             hd2=new_cmatrix(d1->r,d1->c,"");
  297.             if (error) return new_string("Fehler",6,"");
  298.             m2=matrixof(hd2);
  299.             n=(d1->c)*(d1->r);
  300.             for (i=0; i<n; i++)
  301.                 {    fc(&x,&y,m1,m1+1,m2,m2+1); m1+=2; m2+=2;
  302.                 }
  303.             return hd2;
  304.         }
  305.     }
  306.     else if (hd->type==s_matrix)
  307.     {    if (hd1->type==s_matrix)
  308.         {    d=dimsof(hd);
  309.             m=matrixof(hd);
  310.             d1=dimsof(hd1);
  311.             m1=matrixof(hd1);
  312.             if (d->c!=d1->c || d->r!= d1->r)
  313.             {    output("Matrix dimensions must agree!\n");
  314.                 error=4;
  315.                 return new_string("Fehler",6,"");
  316.             }
  317.             hd2=new_matrix(d->r,d->c,"");
  318.             if (error) return new_string("Fehler",6,"");
  319.             m2=matrixof(hd2);
  320.             n=(d->c)*(d->r);
  321.             for (i=0; i<n; i++)
  322.                 {    f(m,m1,m2); m++; m1++; m2++;
  323.                 }
  324.             return hd2;
  325.         }
  326.         if (fc && hd1->type==s_cmatrix)
  327.         {    d=dimsof(hd);
  328.             m=matrixof(hd);
  329.             d1=dimsof(hd1);
  330.             m1=matrixof(hd1);
  331.             if (d->c!=d1->c || d->r!= d1->r)
  332.             {    output("Matrix dimensions must agree!\n");
  333.                 error=4; return new_string("Fehler",6,"");
  334.             }
  335.             hd2=new_cmatrix(d->r,d->c,"");
  336.             if (error) return new_string("Fehler",6,"");
  337.             m2=matrixof(hd2);
  338.             n=(d->c)*(d->r);
  339.             for (i=0; i<n; i++)
  340.                 {    fc(m,&null,m1,m1+1,m2,m2+1); m++; m1+=2; m2+=2;
  341.                 }
  342.             return hd2;
  343.         }
  344.         return mapt2(f,fc,hd1,hd);
  345.     }
  346.     else if (fc && hd->type==s_cmatrix)
  347.     {    if (hd1->type==s_cmatrix)
  348.         {    d=dimsof(hd);
  349.             m=matrixof(hd);
  350.             d1=dimsof(hd1);
  351.             m1=matrixof(hd1);
  352.             if (d->c!=d1->c || d->r!= d1->r)
  353.             {    output("Matrix dimensions must agree!\n");
  354.                 error=4; return new_string("Fehler",6,"");
  355.             }
  356.             hd2=new_cmatrix(d->r,d->c,"");
  357.             if (error) return new_string("Fehler",6,"");
  358.             m2=matrixof(hd2);
  359.             n=(d->c)*(d->r);
  360.             for (i=0; i<n; i++)
  361.                 {    fc(m,m+1,m1,m1+1,m2,m2+1); m+=2; m1+=2; m2+=2;
  362.                 }
  363.             return hd2;
  364.         }
  365.         return mapt2(f,fc,hd1,hd);
  366.     }
  367.     output("Illegal operation\n"); error=3;
  368.     return new_string("Fehler",6,"");
  369. }
  370.  
  371. header *map1 (void f(double *, double *), 
  372.     void fc(double *, double *, double *, double *),
  373.     header *hd)
  374. /***** map
  375.     do the function elementwise to the value.
  376.     te value may be real or complex!
  377. ******/
  378. {    double x,y;
  379.     dims *d;
  380.     header *hd1;
  381.     double *m,*m1;
  382.     int i,n;
  383.     if (hd->type==s_real)
  384.     {    f(realof(hd),&x);
  385.         return new_real(x,"");
  386.     }
  387.     else if (hd->type==s_matrix)
  388.     {    d=dimsof(hd);
  389.         hd1=new_matrix(d->r,d->c,"");
  390.         if (error) return new_string("Fehler",6,"");
  391.         m=matrixof(hd);
  392.         m1=matrixof(hd1);
  393.         n=d->c*d->r;
  394.         for (i=0; i<n; i++)
  395.             {    f(m,m1); m++; m1++;
  396.             }
  397.         return hd1;
  398.     }
  399.     else if (fc && hd->type==s_complex)
  400.     {    fc(realof(hd),imagof(hd),&x,&y);
  401.         return new_complex(x,y,"");
  402.     }
  403.     else if (fc && hd->type==s_cmatrix)
  404.     {    d=dimsof(hd);
  405.         hd1=new_cmatrix(d->r,d->c,"");
  406.         if (error) return new_string("Fehler",6,"");
  407.         m=matrixof(hd);
  408.         m1=matrixof(hd1);
  409.         n=d->c*d->r;
  410.         for (i=0; i<n; i++)
  411.             {    fc(m,m+1,m1,m1+1); m+=2; m1+=2;
  412.             }
  413.         return hd1;
  414.     }
  415.     output("Illegal operation\n"); error=3;
  416.     return new_string("Fehler",6,"");
  417. }
  418.  
  419. header *map1r (void f(double *, double *), 
  420.     void fc(double *, double *, double *),
  421.     header *hd)
  422. /***** map
  423.     do the function elementwise to the value.
  424.     te value may be real or complex! the result is always real.
  425. ******/
  426. {    double x;
  427.     dims *d;
  428.     header *hd1;
  429.     double *m,*m1;
  430.     int i,n;
  431.     if (hd->type==s_real)
  432.     {    f(realof(hd),&x);
  433.         return new_real(x,"");
  434.     }
  435.     else if (hd->type==s_matrix)
  436.     {    d=dimsof(hd);
  437.         hd1=new_matrix(d->r,d->c,"");
  438.         if (error) return new_string("Fehler",6,"");
  439.         m=matrixof(hd);
  440.         m1=matrixof(hd1);
  441.         n=d->c*d->r;
  442.         for (i=0; i<n; i++)
  443.             {    f(m,m1); m++; m1++;
  444.             }
  445.         return hd1;
  446.     }
  447.     else if (fc && hd->type==s_complex)
  448.     {    fc(realof(hd),imagof(hd),&x);
  449.         return new_real(x,"");
  450.     }
  451.     else if (fc && hd->type==s_cmatrix)
  452.     {    d=dimsof(hd);
  453.         hd1=new_matrix(d->r,d->c,"");
  454.         if (error) return new_string("Fehler",6,"");
  455.         m=matrixof(hd);
  456.         m1=matrixof(hd1);
  457.         n=d->c*d->r;
  458.         for (i=0; i<n; i++)
  459.             {    fc(m,m+1,m1); m+=2; m1++;
  460.             }
  461.         return hd1;
  462.     }
  463.     output("Illegal operation\n"); error=3;
  464.     return new_string("Fehler",6,"");
  465. }
  466.  
  467. void real_add (double *x, double *y, double *z)
  468. {    *z=*x+*y;
  469. }
  470.  
  471. void complex_add (double *x, double *xi, double *y, double *yi,
  472.     double *z, double *zi)
  473. {    *z=*x+*y;
  474.     *zi=*xi+*yi;
  475. }
  476.  
  477. void add (header *hd, header *hd1)
  478. /***** add
  479.     add the values.
  480. *****/
  481. {    header *result,*st=hd;
  482.     hd=getvalue(hd); hd1=getvalue(hd1);
  483.     if (error) return;
  484.     result=map2(real_add,complex_add,hd,hd1);
  485.     if (!error) moveresult(st,result);
  486. }
  487.  
  488. void real_subtract (double *x, double *y, double *z)
  489. {    *z=*x-*y;
  490. }
  491.  
  492. void complex_subtract (double *x, double *xi, double *y, double *yi,
  493.     double *z, double *zi)
  494. {    *z=*x-*y;
  495.     *zi=*xi-*yi;
  496. }
  497.  
  498. void subtract (header *hd, header *hd1)
  499. /***** add
  500.     add the values.
  501. *****/
  502. {    header *result,*st=hd;
  503.     hd=getvalue(hd); hd1=getvalue(hd1);
  504.     if (error) return;
  505.     result=map2(real_subtract,complex_subtract,hd,hd1);
  506.     if (!error) moveresult(st,result);
  507. }
  508.  
  509. void real_multiply (double *x, double *y, double *z)
  510. {    *z=*x*(*y);
  511. }
  512.  
  513. void complex_multiply (double *x, double *xi, double *y, double *yi,
  514.     double *z, double *zi)
  515. {    *z=*x * *y - *xi * *yi;
  516.     *zi=*x * *yi + *xi * *y;
  517. }
  518.  
  519. void dotmultiply (header *hd, header *hd1)
  520. /***** add
  521.     add the values.
  522. *****/
  523. {    header *result,*st=hd;
  524.     hd=getvalue(hd); hd1=getvalue(hd1);
  525.     if (error) return;
  526.     result=map2(real_multiply,complex_multiply,hd,hd1);
  527.     if (!error) moveresult(st,result);
  528. }
  529.  
  530. void real_divide (double *x, double *y, double *z)
  531. {    
  532. #ifdef FLOAT_TEST
  533.     if (*y==0.0) { *y=1e-10; error=1; }
  534. #endif
  535.     *z=*x/(*y);
  536. }
  537.  
  538. void complex_divide (double *x, double *xi, double *y, double *yi,
  539.     double *z, double *zi)
  540. {    double r;
  541.     r=*y * *y + *yi * *yi;
  542. #ifdef FLOAT_TEST
  543.     if (r==0) { r=1e-10; error=1; }
  544. #endif
  545.     *z=(*x * *y + *xi * *yi)/r;
  546.     *zi=(*xi * *y - *x * *yi)/r;
  547. }
  548.  
  549. void dotdivide (header *hd, header *hd1)
  550. /***** add
  551.     add the values.
  552. *****/
  553. {    header *result,*st=hd;
  554.     hd=getvalue(hd); hd1=getvalue(hd1);
  555.     if (error) return;
  556.     result=map2(real_divide,complex_divide,hd,hd1);
  557.     if (!error) moveresult(st,result);
  558. }
  559.  
  560. void cscalp (double *s, double *si, double *x, double *xi,
  561.     double *y, double *yi)
  562. {    *s += *x * *y - *xi * *yi;
  563.     *si += *x * *yi + *xi * *y;
  564. }
  565.  
  566. void ccopy (double *y, double *x, double *xi)
  567. {    *y++=*x; *y=*xi;
  568. }
  569.  
  570. void multiply (header *hd, header *hd1)
  571. /***** multiply
  572.     matrix multiplication.
  573. *****/
  574. {    header *result,*st=hd;
  575.     dims *d,*d1;
  576.     double *m,*m1,*m2,*mm1,*mm2,x,y,null=0.0;
  577.     int i,j,c,r,k;
  578.     hd=getvalue(hd); hd1=getvalue(hd1);
  579.     if (error) return;
  580.     if (hd->type==s_matrix && hd1->type==s_matrix)
  581.     {    d=dimsof(hd);
  582.         d1=dimsof(hd1);
  583.         if (d->c != d1->r)
  584.         {    error=8; output("Cannot multiply these!\n");
  585.             return;
  586.         }
  587.         r=d->r; c=d1->c;
  588.         result=new_matrix(r,c,"");
  589.         if (error) return;
  590.         m=matrixof(result);
  591.         m1=matrixof(hd);
  592.         m2=matrixof(hd1);
  593.         for (i=0; i<r; i++)
  594.             for (j=0; j<c; j++)
  595.             {    mm1=mat(m1,d->c,i,0); mm2=m2+j;
  596.                 x=0.0;
  597.                 for (k=0; k<d->c; k++)
  598.                 {    x+=(*mm1)*(*mm2);
  599.                     mm1++; mm2+=d1->c;
  600.                 }
  601.                 *mat(m,c,i,j)=x;
  602.             }
  603.         moveresult(st,result);
  604.         return;
  605.     }
  606.     if (hd->type==s_matrix && hd1->type==s_cmatrix)
  607.     {    d=dimsof(hd);
  608.         d1=dimsof(hd1);
  609.         if (d->c != d1->r)
  610.         {    error=8; output("Cannot multiply these!\n");
  611.             return;
  612.         }
  613.         r=d->r; c=d1->c;
  614.         result=new_cmatrix(r,c,"");
  615.         if (error) return;
  616.         m=matrixof(result);
  617.         m1=matrixof(hd);
  618.         m2=matrixof(hd1);
  619.         for (i=0; i<r; i++)
  620.             for (j=0; j<c; j++)
  621.             {   mm1=mat(m1,d->c,i,0); mm2=m2+(LONG)2*j;
  622.                 x=0.0; y=0.0;
  623.                 for (k=0; k<d->c; k++)
  624.                 {    cscalp(&x,&y,mm1,&null,mm2,mm2+1);
  625.                     mm1++; mm2+=2*d1->c;
  626.                 }
  627.                 ccopy(cmat(m,c,i,j),&x,&y);
  628.             }
  629.         moveresult(st,result);
  630.         return;
  631.     }
  632.     if (hd->type==s_cmatrix && hd1->type==s_matrix)
  633.     {    d=dimsof(hd);
  634.         d1=dimsof(hd1);
  635.         if (d->c != d1->r)
  636.         {    error=8; output("Cannot multiply these!\n");
  637.             return;
  638.         }
  639.         r=d->r; c=d1->c;
  640.         result=new_cmatrix(r,c,"");
  641.         if (error) return;
  642.         m=matrixof(result);
  643.         m1=matrixof(hd);
  644.         m2=matrixof(hd1);
  645.         for (i=0; i<r; i++)
  646.             for (j=0; j<c; j++)
  647.             {    mm1=cmat(m1,d->c,i,0); mm2=m2+j;
  648.                 x=0.0; y=0.0;
  649.                 for (k=0; k<d->c; k++)
  650.                 {    cscalp(&x,&y,mm1,mm1+1,mm2,&null);
  651.                     mm1+=2; mm2+=d1->c;
  652.                 }
  653.                 ccopy(cmat(m,c,i,j),&x,&y);
  654.             }
  655.         moveresult(st,result);
  656.         return;
  657.     }
  658.     if (hd->type==s_cmatrix && hd1->type==s_cmatrix)
  659.     {    d=dimsof(hd);
  660.         d1=dimsof(hd1);
  661.         if (d->c != d1->r)
  662.         {    error=8; output("Cannot multiply these!\n");
  663.             return;
  664.         }
  665.         r=d->r; c=d1->c;
  666.         result=new_cmatrix(r,c,"");
  667.         if (error) return;
  668.         m=matrixof(result);
  669.         m1=matrixof(hd);
  670.         m2=matrixof(hd1);
  671.         for (i=0; i<r; i++)
  672.             for (j=0; j<c; j++)
  673.             {   mm1=cmat(m1,d->c,i,0); mm2=m2+(LONG)2*j;
  674.                 x=0.0; y=0.0;
  675.                 for (k=0; k<d->c; k++)
  676.                 {    cscalp(&x,&y,mm1,mm1+1,mm2,mm2+1);
  677.                     mm1+=2; mm2+=2*d1->c;
  678.                 }
  679.                 ccopy(cmat(m,c,i,j),&x,&y);
  680.             }
  681.         moveresult(st,result);
  682.         return;
  683.     }
  684.     else dotmultiply(st,nextof(st));
  685. }
  686.  
  687. void divide (header *hd, header *hd1)
  688. {    dotdivide(hd,hd1);
  689. }
  690.  
  691. void real_invert (double *x, double *y)
  692. {    *y= -*x;
  693. }
  694.  
  695. void complex_invert (double *x, double *xi, double *y, double *yi)
  696. {    *y= -*x;
  697.     *yi= -*xi;
  698. }
  699.  
  700. void invert (header *hd)
  701. /***** invert
  702.     compute -matrix.
  703. *****/
  704. {    header *result,*st=hd;
  705.     hd=getvalue(hd);
  706.     if (error) return;
  707.     result=map1(real_invert,complex_invert,hd);
  708.     if (!error) moveresult(st,result);
  709. }
  710.  
  711. void get_element (int nargs, header *var, header *hd)
  712. /***** get_elements
  713.     get the element of the matrix.
  714. *****/
  715. {    header *st=hd,*result,*hd1;
  716.     var=getvalue(var); if (error) return;
  717.     if (var->type==s_string) /* interpret the string as a function */
  718.     {    if (exec_builtin(stringof(var),nargs,hd));
  719.         else
  720.         {    hd1=searchudf(stringof(var));
  721.             if (hd1) interpret_udf(hd1,hd,nargs);
  722.             else
  723.             {    output1("%s is not function name!\n",stringof(var));
  724.                 error=2020; return;
  725.             }
  726.         }
  727.         return;
  728.     }
  729.     hd=getvalue(hd); if (error) return;
  730.     if (nargs<1 || nargs>2) 
  731.     {     error=30; output("Illegal matrix reference!\n"); return; }
  732.     if (nargs==2) 
  733.     {    hd1=getvalue(next_param(st)); if (error) return;
  734.     }
  735.     else 
  736.     {    if (dimsof(var)->r==1) { hd1=hd; hd=new_real(1.0,""); }
  737.         else hd1=new_command(c_allv);
  738.         if (error) return;
  739.     }
  740.     if (var->type==s_matrix || var->type==s_real)
  741.     {    result=new_submatrix(var,hd,hd1,"");
  742.     }
  743.     else if (var->type==s_cmatrix || var->type==s_complex)
  744.     {    result=new_csubmatrix(var,hd,hd1,"");
  745.     }
  746.     else
  747.     {    error=31; output1("Not a matrix %s!\n",var->name); return;
  748.     }
  749.     if (error) return;
  750.     moveresult(st,result);
  751. }
  752.  
  753. void get_element1 (char *name, header *hd)
  754. /* get an element of a matrix, referenced by *realof(hd),
  755.    where the matrix is dentified with a vector of same length
  756. */
  757. {    header *st=hd,*result,*var;
  758.     LONG n,l;
  759.     int r,c;
  760.     double *m;
  761.     hd=getvalue(hd);
  762.     var=searchvar(name);
  763.     if (!var)
  764.     {    output1("%s not a variable!\n",name);
  765.         error=1012; return;
  766.     }
  767.     var=getvalue(var); if (error) return;
  768.     if (hd->type!=s_real)
  769.     {    output("Index must be a number!\n");
  770.         error=1013; return;
  771.     }
  772.     if (error) return;
  773.     if (var->type==s_real)
  774.     {    result=new_reference(var,"");
  775.     }
  776.     else if (var->type==s_complex)
  777.     {    result=new_reference(var,"");
  778.     }
  779.     else if (var->type==s_matrix)
  780.     {    getmatrix(var,&r,&c,&m);
  781.         l=(LONG)(*realof(hd));
  782.         n=(LONG)r*c;
  783.         if (l>n) l=n;
  784.         if (l<1) l=1;
  785.         if (nosubmref) result=new_real(*(m+l-1),"");
  786.         else result=new_subm(var,l,"");
  787.     }
  788.     else if (var->type==s_cmatrix)
  789.     {    getmatrix(var,&r,&c,&m);
  790.         l=(LONG)(*realof(hd));
  791.         n=(LONG)r*c;
  792.         if (n==0)
  793.         {    output("Matrix is empty!\n"); error=1030; return;
  794.         }
  795.         if (l>n) l=n;
  796.         if (l<1) l=1;
  797.         if (nosubmref)
  798.         {   m+=(LONG)2*(l-1);
  799.             result=new_complex(*m,*(m+1),"");
  800.         }
  801.         else result=new_csubm(var,l,"");
  802.     }
  803.     else
  804.     {    output1("%s not a variable of proper type for {}!\n");
  805.         error=1011; return;
  806.     }
  807.     moveresult(st,result);
  808. }
  809.  
  810. /****************** udf ************************/
  811.  
  812. void interpret_udf (header *var, header *args, int argn)
  813. /**** interpret_udf
  814.     interpret a user defined function.
  815. ****/
  816. {    int udfold,nargu,i,oldargn,defaults,oldtrace;
  817.     char *oldnext=next,*oldstartlocal,*oldendlocal,*udflineold,*p;
  818.     header *result,*st=args,*hd=args,*hd1,*oldrunning;
  819.     p=helpof(var);
  820.     nargu=*((int *)p); p+=sizeof(int);
  821.     for (i=0; i<argn; i++)
  822.     {    if (hd->type==s_reference && !referenceof(hd))
  823.         {    if (i<nargu && hd->name[0]==0 && *(int *)p)
  824.             {    p+=16+2*sizeof(int);
  825.                 moveresult((header *)newram,(header *)p);
  826.                 p=(char *)nextof((header *)p);
  827.                 hd=nextof(hd);
  828.                 continue;
  829.             }
  830.             else
  831.             {    hd1=getvalue(hd); if (error) return;
  832.             }
  833.         }
  834.         else hd1=hd;
  835.         if (i<nargu) 
  836.         {    defaults=*(int *)p; p+=sizeof(int);
  837.             strcpy(hd1->name,p); hd1->xor=*((int *)(p+16));
  838.             p+=16+sizeof(int);
  839.             if (defaults) p=(char *)nextof((header *)p);
  840.         }
  841.         else
  842.         {    strcpy(hd1->name,argname[i]);
  843.             hd1->xor=xors[i];
  844.         }
  845.         hd=nextof(hd);
  846.     }
  847.     for (i=argn; i<nargu; i++)
  848.     {    defaults=*(int *)p;
  849.         p+=16+2*sizeof(int);
  850.         if (defaults)
  851.         {    moveresult((header *)newram,(header *)p);
  852.             p=(char *)nextof((header *)p);
  853.         }
  854.     }
  855.     udflineold=udfline;
  856.     oldargn=actargn;
  857.     actargn=argn;
  858.     udfline=next=udfof(var); udfold=udfon; udfon=1;
  859.     oldstartlocal=startlocal; oldendlocal=endlocal;
  860.     startlocal=(char *)args; endlocal=newram;
  861.     oldrunning=running; running=var;
  862.     if ((oldtrace=trace)>0)
  863.     {    if (trace==2) trace=0;
  864.         if (trace>0) trace_udfline(next);
  865.     }
  866.     else if (var->flags&1)
  867.     {    trace=1;
  868.         if (trace>0) trace_udfline(next);
  869.     }
  870.     while (!error && udfon==1)
  871.     {    command();
  872.         if (udfon==2)
  873.         {    result=scan_value(); 
  874.             if (error) 
  875.             {    output1("Error in function %s\n",var->name);
  876.                 print_error(udfline);
  877.                 break;
  878.             }
  879.             moveresult1(st,result);
  880.             break;
  881.         }
  882.         if (test_key()==escape) 
  883.         {    output("User interrupted!\n"); error=58; break; 
  884.         }
  885.     }
  886.     endlocal=oldendlocal; startlocal=oldstartlocal;
  887.     running=oldrunning;
  888.     if (trace>=0) trace=oldtrace;
  889.     if (error) output1("Error in function %s\n",var->name);
  890.     if (udfon==0)
  891.     {    output1("Return missing in %s!\n",var->name); error=55; }
  892.     udfon=udfold; next=oldnext; udfline=udflineold;
  893.     actargn=oldargn;
  894. }
  895.  
  896. /***************** scanning ***************************/
  897.  
  898. void copy_complex (double *x, double *y)
  899. {    *x++=*y++;
  900.     *x=*y;
  901. }
  902.  
  903. int scan_arguments (void)
  904. /* look ahead for arguments */
  905. {    int count=0,olds=nosubmref,nocount=0;
  906.     header *hd,*hdold;
  907.     nosubmref=1;
  908.     while (1)
  909.     {    scan_space();
  910.         if (*next==')' || *next==']') break;
  911.         if (*next==',')
  912.             hd=new_reference(0,"");
  913.         else
  914.         {    hd=scan(); scan_space();
  915.         }
  916.         if (*next=='=')
  917.         {    next++;
  918.             hdold=(header *)newram;
  919.             scan_value();
  920.             if (!error)
  921.             {    strcpy(hdold->name,hd->name);
  922.                 hdold->xor=hd->xor;
  923.                 moveresult(hd,hdold);
  924.                 nocount=1;
  925.             }
  926.         }
  927.         else if (nocount)
  928.         {    output("Illegal parameter after named parameter!\n");
  929.             error=2700;
  930.         }
  931.         if (error)
  932.         {    nosubmref=olds;
  933.             return 0;
  934.         }
  935.         while (hd<(header *)newram)
  936.         {    if (!nocount) count++;
  937.             hd=nextof(hd);
  938.         }
  939.         if (count>=10)
  940.         {    output("To many arguments!\n"); error=56; return 0; }
  941.         if (*next!=',') break;
  942.         next++;
  943.     }
  944.     if (*next!=')' && *next!=']')
  945.     {    output("Error in arguments!\n"); error=19; return 0; }
  946.     next++;
  947.     nosubmref=olds;
  948.     return count;
  949. }
  950.  
  951. void scan_matrix (void)
  952. /***** scan_matrix
  953.     scan a matrix from input.
  954.     form: [x y z ... ; v w u ...],
  955.     where x,y,z,u,v,w are sums.
  956. *****/
  957. {    header *hd,*result;
  958.     dims *d;
  959.     int c,r,count,i,j,complex=0;
  960.     LONG ic;
  961.     size_t allcount;
  962.     double *m,*ms,cnull[2]={0,0};
  963.     hd=new_matrix(0,0,""); /* don't know how big it will be! */
  964.     if (error) return;
  965.     count=0;
  966.     getmatrix(hd,&r,&c,&m); ms=m;
  967.     r=0; c=0;
  968.     scan_space();
  969.     while (1)
  970.     {    scan_space();
  971.         if (*next==0) { next_line(); scan_space(); }
  972.             /* matrix is allowed to pass line ends */
  973.         if (*next==';' || *next==']') 
  974.             /* new column starts */
  975.         {    if (*next==';') next++;
  976.             if ((char *)(ms+count*(r+1))>=ramend)
  977.             {    output("Memory overflow!\n"); error=18; return;
  978.             }
  979.             if (count>c) /* this column is to long */
  980.             {    if (r>0) /* expand matrix */
  981.                 {    for (j=count-1; j>=0; j--)
  982.                     {    if (complex) copy_complex(cmat(ms,count,r,j),
  983.                             cmat(ms,c,r,j));
  984.                         else *mat(ms,count,r,j)=*mat(ms,c,r,j);
  985.                     }
  986.                     for (i=r-1; i>=0; i--)
  987.                     {    if (i>0)
  988.                         for (j=c-1; j>=0; j--)
  989.                         {    if (complex) 
  990.                                 copy_complex(cmat(ms,count,i,j),
  991.                                     cmat(ms,c,i,j));
  992.                             else *mat(ms,count,i,j)=*mat(ms,c,i,j);
  993.                         }
  994.                         for (j=c; j<count; j++)
  995.                         if (complex) copy_complex(cmat(ms,count,i,j),
  996.                             cnull);
  997.                         else *mat(ms,count,i,j)=0.0;
  998.                     }
  999.                 }
  1000.                 c=count;
  1001.             }
  1002.             else if (count<c)
  1003.             {    for (j=count; j<c; j++)
  1004.                     if (complex) copy_complex(cmat(ms,c,r,j),
  1005.                             cnull);
  1006.                     else *mat(ms,c,r,j)=0.0;
  1007.             }
  1008.             r++; newram=(char *)(ms+(complex?2l:1l)*(LONG)c*r);
  1009.             m=(double *)newram;
  1010.             count=0;
  1011.         }
  1012.         if (*next==']') break;
  1013.         if (*next==',') next++;
  1014.         if (*next==0) next_line();
  1015.         result=scan_value(); if (error) return;
  1016.         newram=(char *)result;
  1017.         if (!complex && result->type==s_complex)
  1018.         {    complex=1;
  1019.             /* make matrix complex up to now (oh boy!) */
  1020.             allcount=((char *)m-(char *)ms)/sizeof(double);
  1021.             if (newram+allcount*sizeof(double)+result->size>ramend)
  1022.             {    output("Memory overflow!\n"); error=16; return;
  1023.             }
  1024.             if (allcount)
  1025.             {    memmove(newram+allcount*sizeof(double),newram,result->size);
  1026.                 result=(header *)((char *)result+allcount*sizeof(double));
  1027.                 for (ic=allcount-1; ic>0; ic--)
  1028.                 {   *(ms+(LONG)2*ic)=*(ms+ic); 
  1029.                     *(ms+(LONG)2*ic+1)=0.0;
  1030.                 }
  1031.                 *(ms+1)=0.0;
  1032.             newram=(char *)(ms+(LONG)2*allcount); m=(double *)newram;
  1033.             }
  1034.             hd->type=s_cmatrix;
  1035.         }
  1036.         else if (result->type==s_real);
  1037.         else if (result->type==s_complex && complex);
  1038.         else
  1039.         {    error=-1; output("Illegal vector!\n"); return;
  1040.         }
  1041.         *m++=*realof(result); count++;
  1042.         if (complex)
  1043.         {    if (result->type==s_complex) *m++=*imagof(result);
  1044.             else *m++=0.0;
  1045.         }
  1046.         if (count>=INT_MAX)
  1047.         {    output1("Matrix has more than %d columns!\n",INT_MAX);
  1048.             error=17; return;
  1049.         }
  1050.         newram=(char *)m;
  1051.         if (newram>=ramend) 
  1052.         {    output("Memory overflow!\n"); error=16; return; 
  1053.         }
  1054.     }
  1055.     next++;
  1056.     d=(dims *)(hd+1);
  1057.     if (c==0) r=0;
  1058.     d->c=c; d->r=r;
  1059.     if (r>=INT_MAX)
  1060.     {    output1("Matrix has more than %d rows!\n",INT_MAX);
  1061.         error=18; return;
  1062.     }
  1063.     hd->size=complex?cmatrixsize(c,r):matrixsize(c,r);
  1064.     newram=(char *)hd+hd->size;
  1065. }
  1066.  
  1067. void scan_elementary (void)
  1068. /***** scan_elemtary
  1069.     scan an elementary expression, like a value or variable.
  1070.     scans also (...).
  1071. *****/
  1072. {    double x;
  1073.     int n,nargs=0,hadargs=0;
  1074.     header *hd=(header *)newram,*var;
  1075.     char name[16],*s;
  1076.     scan_space();
  1077.     if ((*next=='.' && isdigit(*(next+1))) || isdigit(*next))
  1078.     {    sscanf(next,"%lf%n",&x,&n);
  1079.         next+=n;
  1080.         if (*next=='i') /* complex number! */
  1081.         {    next++;
  1082.             new_complex(0,x,"");
  1083.         }
  1084.         else new_real(x,"");
  1085.     }
  1086.     else if (*next==2) /* a double stored in binary form */
  1087.     {    next++;
  1088. #ifdef SPECIAL_ALIGNMENT
  1089.         memmove((char *)(&x),next,sizeof(double));
  1090. #else
  1091.         x=*((double *)next);
  1092. #endif
  1093.         next+=sizeof(double);
  1094.         if (*next=='i') /* complex number! */
  1095.         {    next++;
  1096.             new_complex(0,x,"");
  1097.         }
  1098.         else new_real(x,"");
  1099.     }
  1100.     else if (*next==3)
  1101.     {    output("Command name used as variable!\n");
  1102.         error=4000; return;
  1103.     }
  1104.     else if (isalpha(*next))
  1105.     {    scan_name(name); if (error) return;
  1106.         scan_space(); nargs=0;
  1107.         if (*next=='{')
  1108.         {    next++; scan(); if (error) return; scan_space();
  1109.             if (*next!='}')
  1110.             {    output("} missing!\n"); error=1010; return;
  1111.             }
  1112.             next++;
  1113.             get_element1(name,hd);
  1114.             goto after;
  1115.         }
  1116.         if (*next=='(' || *next=='[') /* arguments or indices */
  1117.         {    hadargs=(*next=='[')?2:1;
  1118.             next++; nargs=scan_arguments();
  1119.             if (error) return;
  1120.         }
  1121.         if (hadargs==1 && exec_builtin(name,nargs,hd));
  1122.         else
  1123.         {    if (hadargs==2) var=searchvar(name);
  1124.             else if (hadargs==1)
  1125.             {    var=searchudf(name);
  1126.                 if (!var) var=searchvar(name);
  1127.             }
  1128.             else var=searchvar(name);
  1129.             if (var && var->type==s_udf && hadargs==1)
  1130.             {    interpret_udf(var,hd,nargs); if (error) return;
  1131.             }
  1132.             else if (!var && hadargs)
  1133.             {    error=24;
  1134.                 if (hadargs==2)
  1135.                 output1("%s no variable!\n",name);
  1136.                 else
  1137.                 output1(
  1138.             "%s no function or variable, or wrong argument number!\n",
  1139.             name);
  1140.                 return;
  1141.             }
  1142.             else if (var && hadargs)
  1143.             {    get_element(nargs,var,hd);
  1144.             }
  1145.             else hd=new_reference(var,name);
  1146.         }
  1147.     }
  1148.     else if (*next=='#' && *(next+1)!='#')
  1149.     {    next++; mindex(hd);
  1150.     }
  1151.     else if (*next=='+')
  1152.     {    next++; scan_elementary();
  1153.     }
  1154.     else if (*next=='-')
  1155.     {    next++; scan_elementary();
  1156.         if (!error) invert(hd);
  1157.     }
  1158.     else if (*next=='(')
  1159.     {    next++;
  1160.         scan(); if (error) return;
  1161.         scan_space();
  1162.         if (*next!=')')
  1163.         {    output("Closing bracket ) missing!\n");
  1164.             error=5; return;
  1165.         }
  1166.         newram=(char *)nextof(hd);
  1167.         next++;
  1168.     }
  1169.     else if (*next=='[')
  1170.     {    next++;
  1171.         scan_matrix();
  1172.     }
  1173.     else if (*next=='\"')
  1174.     {    next++; s=next; 
  1175.         while (*next!='\"' && *next!=0) next++;
  1176.         hd=new_string(s,next-s,"");
  1177.         if (*next=='\"') next++;
  1178.     }
  1179.     else error=1;
  1180.     after: if (error) return;
  1181.     /* for things, that come after an elementary expression */
  1182.     scan_space();
  1183.     if (*next=='\'') { next++; transpose(hd); }
  1184.     else if (*next=='^' || (*next=='*' && *(next+1)=='*'))
  1185.     {    if (*next=='^') next++; 
  1186.         else next+=2;
  1187.         newram=(char *)nextof(hd);
  1188.         scan_elementary();
  1189.         if (error) return;
  1190.         mpower(hd);
  1191.     }
  1192. }
  1193.  
  1194. void scan_factor (void)
  1195. {    scan_elementary();
  1196. }
  1197.  
  1198. void scan_summand (void)
  1199. {    header *hd=(header *)newram,*hd1;
  1200.     scan_space();
  1201.     scan_factor();
  1202.     if (error) return;
  1203.     while (1)
  1204.     {    hd1=(header *)newram;
  1205.         scan_space();
  1206.         if ((*next=='*' && *(next+1)!='*') 
  1207.                 || (*next=='.' && *(next+1)=='*'))
  1208.         {    if (*next=='*') next++;
  1209.             else next+=2;
  1210.             scan_factor();
  1211.             if (!error) dotmultiply(hd,hd1);
  1212.         }
  1213.         else if (*next=='/' || (*next=='.' && *(next+1)=='/'))
  1214.         {    if (*next=='/') next++;
  1215.             else next+=2;
  1216.             scan_factor();
  1217.             if (!error) dotdivide(hd,hd1);
  1218.         }
  1219.         else if (*next=='.') 
  1220.         {    next++;
  1221.             scan_factor();
  1222.             if (!error) multiply(hd,hd1);
  1223.         }
  1224.         else if (*next=='\\')
  1225.         {    next++;
  1226.             newram=(char *)nextof(hd);
  1227.             scan_factor();
  1228.             if (!error) msolve(hd);
  1229.         }
  1230.         else break;
  1231.         if (error) break;
  1232.     }
  1233. }
  1234.  
  1235. void scan_summe (void)
  1236. {    header *hd=(header *)newram,*hd1;
  1237.     scan_space();
  1238.     scan_summand();
  1239.     if (error) return;
  1240.     while (1)
  1241.     {    hd1=(header *)newram;
  1242.         scan_space();
  1243.         if (*next=='+')
  1244.         {    next++;
  1245.             scan_summand();
  1246.             if (!error) add(hd,hd1);
  1247.         }
  1248.         else if (*next=='-')
  1249.         {    next++;
  1250.             scan_summand();
  1251.             if (!error) subtract(hd,hd1);
  1252.         }
  1253.         else break;
  1254.         if (error) break;
  1255.     }
  1256. }
  1257.  
  1258. void scan_dp (void)
  1259. {    header *hd=(header *)newram,*hd1,*hd2;
  1260.     scan_space();
  1261.     if (*next==':')
  1262.     {    next++;
  1263.         new_command(c_allv);
  1264.         return;
  1265.     }
  1266.     scan_summe();
  1267.     if (*next==':') /* a vector a:b:c or a:c */
  1268.     {    next++;
  1269.         hd1=(header *)newram; scan_summe();
  1270.         if (error) return;
  1271.         scan_space();
  1272.         if (*next==':')
  1273.         {    next++; hd2=(header *)newram; 
  1274.             scan_summe(); if (error) return;
  1275.         }
  1276.         else
  1277.         {    hd2=hd1; hd1=new_real(1.0,"");
  1278.         }
  1279.         if (error) return;
  1280.         vectorize(hd,hd1,hd2);
  1281.     }
  1282. }
  1283.  
  1284. void scan_compare (void)
  1285. {    header *hd=(header *)newram;
  1286.     scan_space();
  1287.     if (*next=='!')
  1288.     {    next++; scan_compare(); mnot(hd); return;
  1289.     }
  1290.     scan_dp(); if (error) return;
  1291.     scan_space();
  1292.     if (*next=='<')
  1293.     {    next++;
  1294.         newram=(char *)nextof(hd);
  1295.         if (*next=='=')
  1296.         {    next++; scan_dp(); if (error) return; mlesseq(hd); return;
  1297.         }
  1298.         else if (*next=='>')
  1299.         {    next++; scan_dp(); if (error) return; munequal(hd); return;
  1300.         }
  1301.         scan_dp(); if (error) return;
  1302.         mless(hd);
  1303.     }
  1304.     else if (*next=='>')
  1305.     {    next++; 
  1306.         newram=(char *)nextof(hd);
  1307.         if (*next=='=')
  1308.         {    next++; scan_dp(); if (error) return; mgreatereq(hd); return;
  1309.         }
  1310.         scan_dp(); if (error) return;
  1311.         mgreater(hd);
  1312.     }
  1313.     else if (*next=='=' && *(next+1)=='=')
  1314.     {    next+=2;
  1315.         newram=(char *)nextof(hd);
  1316.         scan_dp(); if (error) return;
  1317.         mequal(hd);
  1318.     }
  1319.     else if (*next=='~' && *(next+1)=='=')
  1320.     {    next+=2; 
  1321.         newram=(char *)nextof(hd);
  1322.         scan_dp(); if (error) return;
  1323.         maboutequal(hd);
  1324.     }
  1325.     else if (*next=='!' && *(next+1)=='=')
  1326.     {    next+=2; 
  1327.         newram=(char *)nextof(hd);
  1328.         scan_dp(); if (error) return;
  1329.         munequal(hd);
  1330.     }
  1331.     else if (*next=='_')
  1332.     {    next++;
  1333.         newram=(char *)nextof(hd);
  1334.         scan_compare(); if (error) return;
  1335.         mvconcat(hd);
  1336.     }
  1337.     else if (*next=='|' && *(next+1)!='|')
  1338.     {    next++;
  1339.         newram=(char *)nextof(hd); 
  1340.         scan_compare(); if (error) return;
  1341.         mhconcat(hd);
  1342.     }
  1343. }
  1344.  
  1345. void scan_logical (void)
  1346. {    header *hd=(header *)newram;
  1347.     scan_compare(); if (error) return;
  1348.     scan_space();
  1349.     if (*next=='|' && *(next+1)=='|')
  1350.     {    next+=2; 
  1351.         newram=(char *)nextof(hd);
  1352.         scan_compare(); if (error) return;
  1353.         mor(hd);
  1354.     }
  1355.     else if (*next=='&' && *(next+1)=='&')
  1356.     {    next+=2; 
  1357.         newram=(char *)nextof(hd);
  1358.         scan_compare(); if (error) return;
  1359.         mand(hd);
  1360.     }
  1361. }
  1362.  
  1363. header *scan (void)
  1364. {    header *result=(header *)newram;
  1365.     scan_space();
  1366.     if (*next=='{')
  1367.     {    next++; scan_logical(); if (error) return result;
  1368.         while (1)
  1369.         {    scan_space();
  1370.             if (*next=='}') { next++; return result; }
  1371.             if (*next!=',')
  1372.             {    output("Error in {}!\n"); error=104; return result;
  1373.             }
  1374.             next++; scan_logical();
  1375.             if (error) return result;
  1376.         }
  1377.     }
  1378.     else
  1379.     {    scan_logical();
  1380.     }
  1381.     return result;
  1382. }
  1383.  
  1384. header *scan_value (void)
  1385. {    header *result=(header *)newram,*hd,*hd1,*marker,*nextresult,
  1386.         *endresults;
  1387.     int oldnosubmref;
  1388.     size_t size;
  1389.     scan_space();
  1390.     if (*next=='{')
  1391.     {    next++; 
  1392.         oldnosubmref=nosubmref; nosubmref=1; 
  1393.         scan_logical(); nosubmref=oldnosubmref;
  1394.         hd1=getvalue(result);
  1395.         if (error) return result;
  1396.         moveresult(result,hd1);
  1397.         while (1)
  1398.         {    scan_space();
  1399.             if (*next=='}') { next++; return result; }
  1400.             if (*next!=',')
  1401.             {    output("Error in {}!\n"); error=104; return result;
  1402.             }
  1403.             next++; hd=(header *)newram; scan_value();
  1404.             if (error) return result;
  1405.             hd1=getvalue(hd); if (error) return result;
  1406.             moveresult(hd,hd1);
  1407.         }
  1408.     }
  1409.     else
  1410.     {    scan_logical();
  1411.         marker=result;
  1412.         endresults=(header *)newram;
  1413.         while (marker<endresults)
  1414.         {    hd1=getvalue(marker);
  1415.             if (error) return result;
  1416.             nextresult=nextof(marker);
  1417.             if (hd1!=marker)
  1418.             {    if (nextresult==endresults)
  1419.                 {    memmove((char *)marker,(char *)hd1,hd1->size);
  1420.                     newram=(char *)nextof(marker);
  1421.                     break;
  1422.                 }
  1423.                 size=hd1->size-marker->size;
  1424.                 memmove((char *)nextresult+size,(char *)nextresult,
  1425.                     newram-(char *)nextresult);
  1426.                 if (hd1>nextresult) 
  1427.                     hd1=(header *)((char *)hd1+size);
  1428.                 nextresult=(header *)((char *)nextresult+size);
  1429.                 endresults=(header *)((char *)endresults+size);
  1430.                 memmove((char *)marker,(char *)hd1,hd1->size);
  1431.                 newram=(char *)endresults;
  1432.             }
  1433.             marker=nextresult;
  1434.         }
  1435.     }
  1436.     return result;
  1437. }
  1438.