home *** CD-ROM | disk | FTP | other *** search
/ Magazyn Amiga 13 / MA_Cover_13.bin / source / c / apl / ae.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-11-23  |  1.4 KB  |  92 lines
  1. #include "apl.h"
  2.  
  3.  
  4. char base_com[] = {ADD, MUL};
  5.  
  6. ex_base()
  7. {
  8.     struct item *extend();
  9.     struct item *p, *q;
  10.     int i;
  11.     char *savpcp;
  12.     data d1, d2;
  13.  
  14.     p = fetch2();
  15.     q = sp[-2];
  16.     if(scalar(p)){
  17.         if(q->rank > 0) i = q->dim[0];
  18.         else i = q->size;
  19.         q = extend(DA, i, p->datap[0]);
  20.         pop();
  21.         *sp++ = p = q;
  22.         q = sp[-2];
  23.     }
  24.     d1 = p->datap[p->size-1];
  25.     p->datap[p->size-1] = 1.0;
  26.     for(i=p->size-2; i>= 0; i--){
  27.         d2 = p->datap[i];
  28.         p->datap[i] = d1;
  29.         d1 *= d2;
  30.     }
  31.     savpcp = pcp;
  32.     pcp = base_com;
  33.     ex_iprod();
  34.     pcp = savpcp;
  35. }
  36.  
  37. ex_rep()
  38. {
  39.     struct item *p, *q, *r;
  40.     double d1, d2, d3;
  41.     data *p1, *p2, *p3;
  42.  
  43.     p = fetch2();
  44.     q = sp[-2];
  45.     /*
  46.       * first map 1 element vectors to scalars:
  47.      *
  48.     if(scalar(p)) p->rank = 0;
  49.     if(scalar(q)) q->rank = 0;
  50.      */
  51.     r = newdat(DA,  p->rank+q->rank,  p->size*q->size);
  52.     copy(IN, p->dim, r->dim, p->rank);
  53.     copy(IN, q->dim, r->dim+p->rank, q->rank);
  54.     p3 = &r->datap[r->size];
  55.     for(p1 = &p->datap[p->size]; p1 > p->datap; ){
  56.         d1 = *--p1;
  57.         if(d1 == 0.0) d1 = 1.0e38;    /* all else goes here */
  58.         for(p2 = &q->datap[q->size]; p2 > q->datap; ){
  59.             d2 = *--p2;
  60.             d3 = d2 /= d1;
  61.             *p2 = d2 = floor(d2);
  62.             *--p3 = (d3 - d2)*d1;
  63.         }
  64.     }
  65.     pop();
  66.     pop();
  67.     *sp++ = r;
  68. }
  69.  
  70. /*
  71.  * scalar -- return true if arg is a scalar 
  72.  */
  73. scalar(aip)
  74. struct item *aip;
  75. {
  76.     return(aip->size == 1);
  77. }
  78.  
  79.  
  80. struct item *
  81. extend(ty, n, d)
  82. data d;
  83. {
  84.     int i;
  85.     struct item *q;
  86.  
  87.     if(ty != DA) error("extend T");
  88.     q = newdat(ty, 1, n);
  89.     for(i=0; i<n; i++) q->datap[i] = d;
  90.     return(q);
  91. }
  92.