"Fossies" - the Fresh Open Source Software Archive

Member "cfitsio-4.0.0/eval.y" (20 May 2021, 185884 Bytes) of package /linux/misc/cfitsio-4.0.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Bison source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. See also the latest Fossies "Diffs" side-by-side code changes report for "eval.y": 3.49_vs_4.0.0.

    1 %{
    2 /************************************************************************/
    3 /*                                                                      */
    4 /*                       CFITSIO Lexical Parser                         */
    5 /*                                                                      */
    6 /* This file is one of 3 files containing code which parses an          */
    7 /* arithmetic expression and evaluates it in the context of an input    */
    8 /* FITS file table extension.  The CFITSIO lexical parser is divided    */
    9 /* into the following 3 parts/files: the CFITSIO "front-end",           */
   10 /* eval_f.c, contains the interface between the user/CFITSIO and the    */
   11 /* real core of the parser; the FLEX interpreter, eval_l.c, takes the   */
   12 /* input string and parses it into tokens and identifies the FITS       */
   13 /* information required to evaluate the expression (ie, keywords and    */
   14 /* columns); and, the BISON grammar and evaluation routines, eval_y.c,  */
   15 /* receives the FLEX output and determines and performs the actual      */
   16 /* operations.  The files eval_l.c and eval_y.c are produced from       */
   17 /* running flex and bison on the files eval.l and eval.y, respectively. */
   18 /* (flex and bison are available from any GNU archive: see www.gnu.org) */
   19 /*                                                                      */
   20 /* The grammar rules, rather than evaluating the expression in situ,    */
   21 /* builds a tree, or Nodal, structure mapping out the order of          */
   22 /* operations and expression dependencies.  This "compilation" process  */
   23 /* allows for much faster processing of multiple rows.  This technique  */
   24 /* was developed by Uwe Lammers of the XMM Science Analysis System,     */
   25 /* although the CFITSIO implementation is entirely code original.       */
   26 /*                                                                      */
   27 /*                                                                      */
   28 /* Modification History:                                                */
   29 /*                                                                      */
   30 /*   Kent Blackburn      c1992  Original parser code developed for the  */
   31 /*                              FTOOLS software package, in particular, */
   32 /*                              the fselect task.                       */
   33 /*   Kent Blackburn      c1995  BIT column support added                */
   34 /*   Peter D Wilson   Feb 1998  Vector column support added             */
   35 /*   Peter D Wilson   May 1998  Ported to CFITSIO library.  User        */
   36 /*                              interface routines written, in essence  */
   37 /*                              making fselect, fcalc, and maketime     */
   38 /*                              capabilities available to all tools     */
   39 /*                              via single function calls.              */
   40 /*   Peter D Wilson   Jun 1998  Major rewrite of parser core, so as to  */
   41 /*                              create a run-time evaluation tree,      */
   42 /*                              inspired by the work of Uwe Lammers,    */
   43 /*                              resulting in a speed increase of        */
   44 /*                              10-100 times.                           */
   45 /*   Peter D Wilson   Jul 1998  gtifilter(a,b,c,d) function added       */
   46 /*   Peter D Wilson   Aug 1998  regfilter(a,b,c,d) function added       */
   47 /*   Peter D Wilson   Jul 1999  Make parser fitsfile-independent,       */
   48 /*                              allowing a purely vector-based usage    */
   49 /*  Craig B Markwardt Jun 2004  Add MEDIAN() function                   */
   50 /*  Craig B Markwardt Jun 2004  Add SUM(), and MIN/MAX() for bit arrays */
   51 /*  Craig B Markwardt Jun 2004  Allow subscripting of nX bit arrays     */
   52 /*  Craig B Markwardt Jun 2004  Implement statistical functions         */
   53 /*                              NVALID(), AVERAGE(), and STDDEV()       */
   54 /*                              for integer and floating point vectors  */
   55 /*  Craig B Markwardt Jun 2004  Use NULL values for range errors instead*/
   56 /*                              of throwing a parse error               */
   57 /*  Craig B Markwardt Oct 2004  Add ACCUM() and SEQDIFF() functions     */
   58 /*  Craig B Markwardt Feb 2005  Add ANGSEP() function                   */
   59 /*  Craig B Markwardt Aug 2005  CIRCLE, BOX, ELLIPSE, NEAR and REGFILTER*/
   60 /*                              functions now accept vector arguments   */
   61 /*  Craig B Markwardt Sum 2006  Add RANDOMN() and RANDOMP() functions   */
   62 /*  Craig B Markwardt Mar 2007  Allow arguments to RANDOM and RANDOMN to*/
   63 /*                              determine the output dimensions         */
   64 /*  Craig B Markwardt Aug 2009  Add substring STRMID() and string search*/
   65 /*                              STRSTR() functions; more overflow checks*/
   66 /*  Craig B Markwardt Dec 2019  Add bit/hex/oct literal strings and     */
   67 /*                              bitwise operatiosn between integers     */
   68 /*  Craig B Markwardt Mar 2021  Add SETNULL() function                  */
   69 /*                                                                      */
   70 /************************************************************************/
   71 
   72 #define  APPROX 1.0e-7
   73 #include "eval_defs.h"
   74 #include "region.h"
   75 #include <time.h>
   76 
   77 #include <stdlib.h>
   78 
   79 #ifndef alloca
   80 #define alloca malloc
   81 #endif
   82 
   83 /* Random number generators for various distributions */
   84 #include "simplerng.h"
   85 
   86    /*  Shrink the initial stack depth to keep local data <32K (mac limit)  */
   87    /*  yacc will allocate more space if needed, though.                    */
   88 #define  YYINITDEPTH   100
   89 
   90 /***************************************************************/
   91 /*  Replace Bison's BACKUP macro with one that fixes a bug --  */
   92 /*  must update state after popping the stack -- and allows    */
   93 /*  popping multiple terms at one time.                        */
   94 /***************************************************************/
   95 
   96 #define YYNEWBACKUP(token, value)                               \
   97    do                               \
   98      if (yychar == YYEMPTY )                    \
   99        { yychar = (token);                                      \
  100          memcpy( &yylval, &(value), sizeof(value) );            \
  101          yychar1 = YYTRANSLATE (yychar);            \
  102          while (yylen--) YYPOPSTACK;                \
  103          yystate = *yyssp;                  \
  104          goto yybackup;                     \
  105        }                            \
  106      else                           \
  107        { yyerror ("syntax error: cannot back up"); YYERROR; }   \
  108    while (0)
  109 
  110 /***************************************************************/
  111 /*  Useful macros for accessing/testing Nodes                  */
  112 /***************************************************************/
  113 
  114 #define TEST(a)        if( (a)<0 ) YYERROR
  115 #define SIZE(a)        gParse.Nodes[ a ].value.nelem
  116 #define TYPE(a)        gParse.Nodes[ a ].type
  117 #define OPER(a)        gParse.Nodes[ a ].operation
  118 #define PROMOTE(a,b)   if( TYPE(a) > TYPE(b) )                  \
  119                           b = New_Unary( TYPE(a), 0, b );       \
  120                        else if( TYPE(a) < TYPE(b) )             \
  121                       a = New_Unary( TYPE(b), 0, a );
  122 
  123 /*****  Internal functions  *****/
  124 
  125 #ifdef __cplusplus
  126 extern "C" {
  127 #endif
  128 
  129 static int  Alloc_Node    ( void );
  130 static void Free_Last_Node( void );
  131 static void Evaluate_Node ( int thisNode );
  132 
  133 static int  New_Const ( int returnType, void *value, long len );
  134 static int  New_Column( int ColNum );
  135 static int  New_Offset( int ColNum, int offset );
  136 static int  New_Unary ( int returnType, int Op, int Node1 );
  137 static int  New_BinOp ( int returnType, int Node1, int Op, int Node2 );
  138 static int  New_Func  ( int returnType, funcOp Op, int nNodes,
  139             int Node1, int Node2, int Node3, int Node4, 
  140             int Node5, int Node6, int Node7 );
  141 static int  New_FuncSize( int returnType, funcOp Op, int nNodes,
  142             int Node1, int Node2, int Node3, int Node4, 
  143               int Node5, int Node6, int Node7, int Size);
  144 static int  New_Deref ( int Var,  int nDim,
  145             int Dim1, int Dim2, int Dim3, int Dim4, int Dim5 );
  146 static int  New_GTI   ( funcOp Op, char *fname, int Node1, int Node2, char *start, char *stop );
  147 static int  New_REG   ( char *fname, int NodeX, int NodeY, char *colNames );
  148 static int  New_Vector( int subNode );
  149 static int  Close_Vec ( int vecNode );
  150 static int  Locate_Col( Node *this );
  151 static int  Test_Dims ( int Node1, int Node2 );
  152 static void Copy_Dims ( int Node1, int Node2 );
  153 
  154 static void Allocate_Ptrs( Node *this );
  155 static void Do_Unary     ( Node *this );
  156 static void Do_Offset    ( Node *this );
  157 static void Do_BinOp_bit ( Node *this );
  158 static void Do_BinOp_str ( Node *this );
  159 static void Do_BinOp_log ( Node *this );
  160 static void Do_BinOp_lng ( Node *this );
  161 static void Do_BinOp_dbl ( Node *this );
  162 static void Do_Func      ( Node *this );
  163 static void Do_Deref     ( Node *this );
  164 static void Do_GTI       ( Node *this );
  165 static void Do_GTI_Over  ( Node *this );
  166 static void Do_REG       ( Node *this );
  167 static void Do_Vector    ( Node *this );
  168 
  169 static long Search_GTI   ( double evtTime, long nGTI, double *start,
  170                double *stop, int ordered, long *nextGTI );
  171 static double GTI_Over(double evtStart, double evtStop,
  172                long nGTI, double *start, double *stop,
  173                long *gtiout);
  174 
  175 static char  saobox (double xcen, double ycen, double xwid, double ywid,
  176              double rot,  double xcol, double ycol);
  177 static char  ellipse(double xcen, double ycen, double xrad, double yrad,
  178              double rot, double xcol, double ycol);
  179 static char  circle (double xcen, double ycen, double rad,
  180              double xcol, double ycol);
  181 static char  bnear  (double x, double y, double tolerance);
  182 static char  bitcmp (char *bitstrm1, char *bitstrm2);
  183 static char  bitlgte(char *bits1, int oper, char *bits2);
  184 
  185 static void  bitand(char *result, char *bitstrm1, char *bitstrm2);
  186 static void  bitor (char *result, char *bitstrm1, char *bitstrm2);
  187 static void  bitnot(char *result, char *bits);
  188 static int cstrmid(char *dest_str, int dest_len,
  189            char *src_str,  int src_len, int pos);
  190 
  191 static void  yyerror(char *msg);
  192 
  193 #ifdef __cplusplus
  194     }
  195 #endif
  196 
  197 %}
  198 
  199 %union {
  200     int    Node;        /* Index of Node */
  201     double dbl;         /* real value    */
  202     long   lng;         /* integer value */
  203     char   log;         /* logical value */
  204     char   str[MAX_STRLEN];    /* string value  */
  205 }
  206 
  207 %token <log>   BOOLEAN        /* First 3 must be in order of        */
  208 %token <lng>   LONG           /* increasing promotion for later use */
  209 %token <dbl>   DOUBLE
  210 %token <str>   STRING
  211 %token <str>   BITSTR
  212 %token <str>   FUNCTION
  213 %token <str>   BFUNCTION      /* Bit function */
  214 %token <str>   IFUNCTION      /* Integer function */
  215 %token <str>   GTIFILTER
  216 %token <str>   GTIOVERLAP
  217 %token <str>   REGFILTER
  218 %token <lng>   COLUMN
  219 %token <lng>   BCOLUMN
  220 %token <lng>   SCOLUMN
  221 %token <lng>   BITCOL
  222 %token <lng>   ROWREF
  223 %token <lng>   NULLREF
  224 %token <lng>   SNULLREF
  225 
  226 %type <Node>  expr
  227 %type <Node>  bexpr
  228 %type <Node>  sexpr
  229 %type <Node>  bits
  230 %type <Node>  vector
  231 %type <Node>  bvector
  232 
  233 %left     ',' '=' ':' '{' '}'
  234 %right    '?'
  235 %left     OR
  236 %left     AND
  237 %left     EQ NE '~'
  238 %left     GT LT LTE GTE
  239 %left     '+' '-' '%'
  240 %left     '*' '/'
  241 %left     '|' '&' XOR
  242 %right    POWER
  243 %left     NOT
  244 %left     INTCAST FLTCAST
  245 %left     UMINUS
  246 %left     '['
  247 
  248 %right    ACCUM DIFF
  249 
  250 %%
  251 
  252 lines:   /* nothing ; was | lines line */
  253        | lines line
  254        ;
  255 
  256 line:           '\n' {}
  257        | expr   '\n'
  258                 { if( $1<0 ) {
  259              yyerror("Couldn't build node structure: out of memory?");
  260              YYERROR;  }
  261                   gParse.resultNode = $1;
  262         }
  263        | bexpr  '\n'
  264                 { if( $1<0 ) {
  265              yyerror("Couldn't build node structure: out of memory?");
  266              YYERROR;  }
  267                   gParse.resultNode = $1;
  268         }
  269        | sexpr  '\n'
  270                 { if( $1<0 ) {
  271              yyerror("Couldn't build node structure: out of memory?");
  272              YYERROR;  } 
  273                   gParse.resultNode = $1;
  274         }
  275        | bits   '\n'
  276                 { if( $1<0 ) {
  277              yyerror("Couldn't build node structure: out of memory?");
  278              YYERROR;  }
  279                   gParse.resultNode = $1;
  280         }
  281        | error  '\n' {  yyerrok;  }
  282        ;
  283 
  284 bvector: '{' bexpr
  285                 { $$ = New_Vector( $2 ); TEST($$); }
  286        | bvector ',' bexpr
  287                 {
  288                   if( gParse.Nodes[$1].nSubNodes >= MAXSUBS ) {
  289              $1 = Close_Vec( $1 ); TEST($1);
  290              $$ = New_Vector( $1 ); TEST($$);
  291                   } else {
  292                      $$ = $1;
  293                   }
  294           gParse.Nodes[$$].SubNodes[ gParse.Nodes[$$].nSubNodes++ ]
  295              = $3;
  296                 }
  297        ;
  298 
  299 vector:  '{' expr
  300                 { $$ = New_Vector( $2 ); TEST($$); }
  301        | vector ',' expr
  302                 {
  303                   if( TYPE($1) < TYPE($3) )
  304                      TYPE($1) = TYPE($3);
  305                   if( gParse.Nodes[$1].nSubNodes >= MAXSUBS ) {
  306              $1 = Close_Vec( $1 ); TEST($1);
  307              $$ = New_Vector( $1 ); TEST($$);
  308                   } else {
  309                      $$ = $1;
  310                   }
  311           gParse.Nodes[$$].SubNodes[ gParse.Nodes[$$].nSubNodes++ ]
  312              = $3;
  313                 }
  314        | vector ',' bexpr
  315                 {
  316                   if( gParse.Nodes[$1].nSubNodes >= MAXSUBS ) {
  317              $1 = Close_Vec( $1 ); TEST($1);
  318              $$ = New_Vector( $1 ); TEST($$);
  319                   } else {
  320                      $$ = $1;
  321                   }
  322           gParse.Nodes[$$].SubNodes[ gParse.Nodes[$$].nSubNodes++ ]
  323              = $3;
  324                 }
  325        | bvector ',' expr
  326                 {
  327                   TYPE($1) = TYPE($3);
  328                   if( gParse.Nodes[$1].nSubNodes >= MAXSUBS ) {
  329              $1 = Close_Vec( $1 ); TEST($1);
  330              $$ = New_Vector( $1 ); TEST($$);
  331                   } else {
  332                      $$ = $1;
  333                   }
  334           gParse.Nodes[$$].SubNodes[ gParse.Nodes[$$].nSubNodes++ ]
  335              = $3;
  336                 }
  337        ;
  338 
  339 expr:    vector '}'
  340                 { $$ = Close_Vec( $1 ); TEST($$); }
  341        ;
  342 
  343 bexpr:   bvector '}'
  344                 { $$ = Close_Vec( $1 ); TEST($$); }
  345        ;
  346 
  347 bits:    BITSTR
  348                 {
  349                   $$ = New_Const( BITSTR, $1, strlen($1)+1 ); TEST($$);
  350           SIZE($$) = strlen($1); }
  351        | BITCOL
  352                 { $$ = New_Column( $1 ); TEST($$); }
  353        | BITCOL '{' expr '}'
  354                 {
  355                   if( TYPE($3) != LONG
  356               || OPER($3) != CONST_OP ) {
  357              yyerror("Offset argument must be a constant integer");
  358              YYERROR;
  359           }
  360                   $$ = New_Offset( $1, $3 ); TEST($$);
  361                 }
  362        | bits '&' bits
  363                 { $$ = New_BinOp( BITSTR, $1, '&', $3 ); TEST($$);
  364                   SIZE($$) = ( SIZE($1)>SIZE($3) ? SIZE($1) : SIZE($3) );  }
  365        | bits '|' bits
  366                 { $$ = New_BinOp( BITSTR, $1, '|', $3 ); TEST($$);
  367                   SIZE($$) = ( SIZE($1)>SIZE($3) ? SIZE($1) : SIZE($3) );  }
  368        | bits '+' bits
  369                 { 
  370           if (SIZE($1)+SIZE($3) >= MAX_STRLEN) {
  371             yyerror("Combined bit string size exceeds " MAX_STRLEN_S " bits");
  372             YYERROR;
  373           }
  374           $$ = New_BinOp( BITSTR, $1, '+', $3 ); TEST($$);
  375                   SIZE($$) = SIZE($1) + SIZE($3); 
  376         }
  377        | bits '[' expr ']'
  378                 { $$ = New_Deref( $1, 1, $3,  0,  0,  0,   0 ); TEST($$); }
  379        | bits '[' expr ',' expr ']'
  380                 { $$ = New_Deref( $1, 2, $3, $5,  0,  0,   0 ); TEST($$); }
  381        | bits '[' expr ',' expr ',' expr ']'
  382                 { $$ = New_Deref( $1, 3, $3, $5, $7,  0,   0 ); TEST($$); }
  383        | bits '[' expr ',' expr ',' expr ',' expr ']'
  384                 { $$ = New_Deref( $1, 4, $3, $5, $7, $9,   0 ); TEST($$); }
  385        | bits '[' expr ',' expr ',' expr ',' expr ',' expr ']'
  386                 { $$ = New_Deref( $1, 5, $3, $5, $7, $9, $11 ); TEST($$); }
  387        | NOT bits
  388                 { $$ = New_Unary( BITSTR, NOT, $2 ); TEST($$);     }
  389 
  390        | '(' bits ')'
  391                 { $$ = $2; }
  392        ;
  393 
  394 expr:    LONG
  395                 { $$ = New_Const( LONG,   &($1), sizeof(long)   ); TEST($$); }
  396        | DOUBLE
  397                 { $$ = New_Const( DOUBLE, &($1), sizeof(double) ); TEST($$); }
  398        | COLUMN
  399                 { $$ = New_Column( $1 ); TEST($$); }
  400        | COLUMN '{' expr '}'
  401                 {
  402                   if( TYPE($3) != LONG
  403               || OPER($3) != CONST_OP ) {
  404              yyerror("Offset argument must be a constant integer");
  405              YYERROR;
  406           }
  407                   $$ = New_Offset( $1, $3 ); TEST($$);
  408                 }
  409        | ROWREF
  410                 { $$ = New_Func( LONG, row_fct,  0, 0, 0, 0, 0, 0, 0, 0 ); }
  411        | NULLREF
  412                 { $$ = New_Func( LONG, null_fct, 0, 0, 0, 0, 0, 0, 0, 0 ); }
  413        | expr '%' expr
  414                 { PROMOTE($1,$3); $$ = New_BinOp( TYPE($1), $1, '%', $3 );
  415           TEST($$);                                                }
  416        | expr '+' expr
  417                 { PROMOTE($1,$3); $$ = New_BinOp( TYPE($1), $1, '+', $3 );
  418           TEST($$);                                                }
  419        | expr '-' expr
  420                 { PROMOTE($1,$3); $$ = New_BinOp( TYPE($1), $1, '-', $3 ); 
  421           TEST($$);                                                }
  422        | expr '*' expr
  423                 { PROMOTE($1,$3); $$ = New_BinOp( TYPE($1), $1, '*', $3 ); 
  424           TEST($$);                                                }
  425        | expr '/' expr
  426                 { PROMOTE($1,$3); $$ = New_BinOp( TYPE($1), $1, '/', $3 ); 
  427           TEST($$);                                                }
  428        | expr '&' expr
  429                 { 
  430                    if (TYPE($1) != LONG ||
  431                TYPE($3) != LONG) {
  432                      yyerror("Bitwise operations with incompatible types; only (bit OP bit) and (int OP int) are allowed");
  433                       YYERROR;
  434                    }
  435                    $$ = New_BinOp( TYPE($1), $1, '&', $3 );
  436                 }
  437        | expr '|' expr
  438                 { 
  439                    if (TYPE($1) != LONG ||
  440                TYPE($3) != LONG) {
  441                      yyerror("Bitwise operations with incompatible types; only (bit OP bit) and (int OP int) are allowed");
  442                       YYERROR;
  443                    }
  444                    $$ = New_BinOp( TYPE($1), $1, '|', $3 );
  445                 }
  446        | expr XOR expr
  447                 { 
  448                    if (TYPE($1) != LONG ||
  449                TYPE($3) != LONG) {
  450                      yyerror("Bitwise operations with incompatible types; only (bit OP bit) and (int OP int) are allowed");
  451                       YYERROR;
  452                    }
  453                    $$ = New_BinOp( TYPE($1), $1, '^', $3 );
  454                 }
  455        | expr POWER expr
  456                 { PROMOTE($1,$3); $$ = New_BinOp( TYPE($1), $1, POWER, $3 );
  457           TEST($$);                                                }
  458        | '+' expr %prec UMINUS
  459                 { $$ = $2; }
  460        | '-' expr %prec UMINUS
  461                 { $$ = New_Unary( TYPE($2), UMINUS, $2 ); TEST($$); }
  462        |  '(' expr ')'
  463                 { $$ = $2; }
  464        | expr '*' bexpr
  465                 { $3 = New_Unary( TYPE($1), 0, $3 );
  466                   $$ = New_BinOp( TYPE($1), $1, '*', $3 ); 
  467           TEST($$);                                }
  468        | bexpr '*' expr
  469                 { $1 = New_Unary( TYPE($3), 0, $1 );
  470                   $$ = New_BinOp( TYPE($3), $1, '*', $3 );
  471                   TEST($$);                                }
  472        | bexpr '?' expr ':' expr
  473                 {
  474                   PROMOTE($3,$5);
  475                   if( ! Test_Dims($3,$5) ) {
  476                      yyerror("Incompatible dimensions in '?:' arguments");
  477              YYERROR;
  478                   }
  479                   $$ = New_Func( 0, ifthenelse_fct, 3, $3, $5, $1,
  480                                  0, 0, 0, 0 );
  481                   TEST($$);
  482                   if( SIZE($3)<SIZE($5) )  Copy_Dims($$, $5);
  483                   TYPE($1) = TYPE($3);
  484                   if( ! Test_Dims($1,$$) ) {
  485                      yyerror("Incompatible dimensions in '?:' condition");
  486              YYERROR;
  487                   }
  488                   TYPE($1) = BOOLEAN;
  489                   if( SIZE($$)<SIZE($1) )  Copy_Dims($$, $1);
  490                 }
  491        | bexpr '?' bexpr ':' expr
  492                 {
  493                   PROMOTE($3,$5);
  494                   if( ! Test_Dims($3,$5) ) {
  495                      yyerror("Incompatible dimensions in '?:' arguments");
  496              YYERROR;
  497                   }
  498                   $$ = New_Func( 0, ifthenelse_fct, 3, $3, $5, $1,
  499                                  0, 0, 0, 0 );
  500                   TEST($$);
  501                   if( SIZE($3)<SIZE($5) )  Copy_Dims($$, $5);
  502                   TYPE($1) = TYPE($3);
  503                   if( ! Test_Dims($1,$$) ) {
  504                      yyerror("Incompatible dimensions in '?:' condition");
  505              YYERROR;
  506                   }
  507                   TYPE($1) = BOOLEAN;
  508                   if( SIZE($$)<SIZE($1) )  Copy_Dims($$, $1);
  509                 }
  510        | bexpr '?' expr ':' bexpr
  511                 {
  512                   PROMOTE($3,$5);
  513                   if( ! Test_Dims($3,$5) ) {
  514                      yyerror("Incompatible dimensions in '?:' arguments");
  515              YYERROR;
  516                   }
  517                   $$ = New_Func( 0, ifthenelse_fct, 3, $3, $5, $1,
  518                                  0, 0, 0, 0 );
  519                   TEST($$);
  520                   if( SIZE($3)<SIZE($5) )  Copy_Dims($$, $5);
  521                   TYPE($1) = TYPE($3);
  522                   if( ! Test_Dims($1,$$) ) {
  523                      yyerror("Incompatible dimensions in '?:' condition");
  524              YYERROR;
  525                   }
  526                   TYPE($1) = BOOLEAN;
  527                   if( SIZE($$)<SIZE($1) )  Copy_Dims($$, $1);
  528                 }
  529        | FUNCTION ')'
  530                 { if (FSTRCMP($1,"RANDOM(") == 0) {  /* Scalar RANDOM() */
  531                      $$ = New_Func( DOUBLE, rnd_fct, 0, 0, 0, 0, 0, 0, 0, 0 );
  532           } else if (FSTRCMP($1,"RANDOMN(") == 0) {/*Scalar RANDOMN()*/
  533              $$ = New_Func( DOUBLE, gasrnd_fct, 0, 0, 0, 0, 0, 0, 0, 0 );
  534                   } else {
  535                      yyerror("Function() not supported");
  536              YYERROR;
  537           }
  538                   TEST($$); 
  539                 }
  540        | FUNCTION bexpr ')'
  541                 { if (FSTRCMP($1,"SUM(") == 0) {
  542              $$ = New_Func( LONG, sum_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
  543                   } else if (FSTRCMP($1,"NELEM(") == 0) {
  544                      $$ = New_Const( LONG, &( SIZE($2) ), sizeof(long) );
  545                   } else if (FSTRCMP($1,"ACCUM(") == 0) {
  546             long zero = 0;
  547             $$ = New_BinOp( LONG , $2, ACCUM, New_Const( LONG, &zero, sizeof(zero) ));
  548           } else {
  549                      yyerror("Function(bool) not supported");
  550              YYERROR;
  551           }
  552                   TEST($$); 
  553         }
  554        | FUNCTION sexpr ')'
  555                 { if (FSTRCMP($1,"NELEM(") == 0) {
  556                      $$ = New_Const( LONG, &( SIZE($2) ), sizeof(long) );
  557           } else if (FSTRCMP($1,"NVALID(") == 0) {
  558              $$ = New_Func( LONG, nonnull_fct, 1, $2,
  559                     0, 0, 0, 0, 0, 0 );
  560           } else {
  561                      yyerror("Function(str) not supported");
  562              YYERROR;
  563           }
  564                   TEST($$); 
  565         }
  566        | FUNCTION bits ')'
  567                 { if (FSTRCMP($1,"NELEM(") == 0) {
  568                      $$ = New_Const( LONG, &( SIZE($2) ), sizeof(long) );
  569         } else if (FSTRCMP($1,"NVALID(") == 0) { /* Bit arrays do not have NULL */
  570                      $$ = New_Const( LONG, &( SIZE($2) ), sizeof(long) );
  571         } else if (FSTRCMP($1,"SUM(") == 0) {
  572              $$ = New_Func( LONG, sum_fct, 1, $2,
  573                     0, 0, 0, 0, 0, 0 );
  574         } else if (FSTRCMP($1,"MIN(") == 0) {
  575              $$ = New_Func( TYPE($2),  /* Force 1D result */
  576                     min1_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
  577              /* Note: $2 is a vector so the result can never
  578                 be a constant.  Therefore it will never be set
  579                 inside New_Func(), and it is safe to set SIZE() */
  580              SIZE($$) = 1;
  581         } else if (FSTRCMP($1,"ACCUM(") == 0) {
  582             long zero = 0;
  583             $$ = New_BinOp( LONG , $2, ACCUM, New_Const( LONG, &zero, sizeof(zero) ));
  584         } else if (FSTRCMP($1,"MAX(") == 0) {
  585              $$ = New_Func( TYPE($2),  /* Force 1D result */
  586                     max1_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
  587              /* Note: $2 is a vector so the result can never
  588                 be a constant.  Therefore it will never be set
  589                 inside New_Func(), and it is safe to set SIZE() */
  590              SIZE($$) = 1;
  591         } else {
  592                      yyerror("Function(bits) not supported");
  593              YYERROR;
  594           }
  595                   TEST($$); 
  596         }
  597        | FUNCTION expr ')'
  598                 { if (FSTRCMP($1,"SUM(") == 0)
  599              $$ = New_Func( TYPE($2), sum_fct, 1, $2,
  600                     0, 0, 0, 0, 0, 0 );
  601           else if (FSTRCMP($1,"AVERAGE(") == 0)
  602              $$ = New_Func( DOUBLE, average_fct, 1, $2,
  603                     0, 0, 0, 0, 0, 0 );
  604           else if (FSTRCMP($1,"STDDEV(") == 0)
  605              $$ = New_Func( DOUBLE, stddev_fct, 1, $2,
  606                     0, 0, 0, 0, 0, 0 );
  607           else if (FSTRCMP($1,"MEDIAN(") == 0)
  608              $$ = New_Func( TYPE($2), median_fct, 1, $2,
  609                     0, 0, 0, 0, 0, 0 );
  610           else if (FSTRCMP($1,"NELEM(") == 0)
  611                      $$ = New_Const( LONG, &( SIZE($2) ), sizeof(long) );
  612           else if (FSTRCMP($1,"NVALID(") == 0)
  613              $$ = New_Func( LONG, nonnull_fct, 1, $2,
  614                     0, 0, 0, 0, 0, 0 );
  615           else if   ((FSTRCMP($1,"ACCUM(") == 0) && (TYPE($2) == LONG)) {
  616             long zero = 0;
  617             $$ = New_BinOp( LONG ,   $2, ACCUM, New_Const( LONG,   &zero, sizeof(zero) ));
  618           } else if ((FSTRCMP($1,"ACCUM(") == 0) && (TYPE($2) == DOUBLE)) {
  619             double zero = 0;
  620             $$ = New_BinOp( DOUBLE , $2, ACCUM, New_Const( DOUBLE, &zero, sizeof(zero) ));
  621           } else if ((FSTRCMP($1,"SEQDIFF(") == 0) && (TYPE($2) == LONG)) {
  622             long zero = 0;
  623             $$ = New_BinOp( LONG ,   $2, DIFF, New_Const( LONG,   &zero, sizeof(zero) ));
  624           } else if ((FSTRCMP($1,"SEQDIFF(") == 0) && (TYPE($2) == DOUBLE)) {
  625             double zero = 0;
  626             $$ = New_BinOp( DOUBLE , $2, DIFF, New_Const( DOUBLE, &zero, sizeof(zero) ));
  627           } else if (FSTRCMP($1,"ABS(") == 0)
  628              $$ = New_Func( 0, abs_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
  629           else if (FSTRCMP($1,"MIN(") == 0)
  630              $$ = New_Func( TYPE($2),  /* Force 1D result */
  631                     min1_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
  632           else if (FSTRCMP($1,"MAX(") == 0)
  633              $$ = New_Func( TYPE($2),  /* Force 1D result */
  634                     max1_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
  635           else if (FSTRCMP($1,"RANDOM(") == 0) { /* Vector RANDOM() */
  636                      $$ = New_Func( 0, rnd_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
  637              TEST($$);
  638              TYPE($$) = DOUBLE;
  639           } else if (FSTRCMP($1,"RANDOMN(") == 0) {
  640              $$ = New_Func( 0, gasrnd_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
  641              TEST($$);
  642              TYPE($$) = DOUBLE;
  643                   } 
  644           else {  /*  These all take DOUBLE arguments  */
  645              if( TYPE($2) != DOUBLE ) $2 = New_Unary( DOUBLE, 0, $2 );
  646                      if (FSTRCMP($1,"SIN(") == 0)
  647             $$ = New_Func( 0, sin_fct,  1, $2, 0, 0, 0, 0, 0, 0 );
  648              else if (FSTRCMP($1,"COS(") == 0)
  649             $$ = New_Func( 0, cos_fct,  1, $2, 0, 0, 0, 0, 0, 0 );
  650              else if (FSTRCMP($1,"TAN(") == 0)
  651             $$ = New_Func( 0, tan_fct,  1, $2, 0, 0, 0, 0, 0, 0 );
  652              else if (FSTRCMP($1,"ARCSIN(") == 0
  653                   || FSTRCMP($1,"ASIN(") == 0)
  654             $$ = New_Func( 0, asin_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
  655              else if (FSTRCMP($1,"ARCCOS(") == 0
  656                   || FSTRCMP($1,"ACOS(") == 0)
  657             $$ = New_Func( 0, acos_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
  658              else if (FSTRCMP($1,"ARCTAN(") == 0
  659                   || FSTRCMP($1,"ATAN(") == 0)
  660             $$ = New_Func( 0, atan_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
  661              else if (FSTRCMP($1,"SINH(") == 0)
  662             $$ = New_Func( 0, sinh_fct,  1, $2, 0, 0, 0, 0, 0, 0 );
  663              else if (FSTRCMP($1,"COSH(") == 0)
  664             $$ = New_Func( 0, cosh_fct,  1, $2, 0, 0, 0, 0, 0, 0 );
  665              else if (FSTRCMP($1,"TANH(") == 0)
  666             $$ = New_Func( 0, tanh_fct,  1, $2, 0, 0, 0, 0, 0, 0 );
  667              else if (FSTRCMP($1,"EXP(") == 0)
  668             $$ = New_Func( 0, exp_fct,  1, $2, 0, 0, 0, 0, 0, 0 );
  669              else if (FSTRCMP($1,"LOG(") == 0)
  670             $$ = New_Func( 0, log_fct,  1, $2, 0, 0, 0, 0, 0, 0 );
  671              else if (FSTRCMP($1,"LOG10(") == 0)
  672             $$ = New_Func( 0, log10_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
  673              else if (FSTRCMP($1,"SQRT(") == 0)
  674             $$ = New_Func( 0, sqrt_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
  675              else if (FSTRCMP($1,"ROUND(") == 0)
  676             $$ = New_Func( 0, round_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
  677              else if (FSTRCMP($1,"FLOOR(") == 0)
  678             $$ = New_Func( 0, floor_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
  679              else if (FSTRCMP($1,"CEIL(") == 0)
  680             $$ = New_Func( 0, ceil_fct, 1, $2, 0, 0, 0, 0, 0, 0 );
  681              else if (FSTRCMP($1,"RANDOMP(") == 0) {
  682                $$ = New_Func( 0, poirnd_fct, 1, $2, 
  683                       0, 0, 0, 0, 0, 0 );
  684                TYPE($$) = LONG;
  685              } else {
  686             yyerror("Function(expr) not supported");
  687             YYERROR;
  688              }
  689           }
  690                   TEST($$); 
  691                 }
  692        | IFUNCTION sexpr ',' sexpr ')'
  693                 { 
  694           if (FSTRCMP($1,"STRSTR(") == 0) {
  695             $$ = New_Func( LONG, strpos_fct, 2, $2, $4, 0, 
  696                    0, 0, 0, 0 );
  697             TEST($$);
  698           }
  699                 }
  700        | FUNCTION expr ',' expr ')'
  701                 { 
  702            if (FSTRCMP($1,"DEFNULL(") == 0) {
  703               if( SIZE($2)>=SIZE($4) && Test_Dims( $2, $4 ) ) {
  704              PROMOTE($2,$4);
  705              $$ = New_Func( 0, defnull_fct, 2, $2, $4, 0,
  706                     0, 0, 0, 0 );
  707              TEST($$); 
  708               } else {
  709              yyerror("Dimensions of DEFNULL arguments "
  710                  "are not compatible");
  711              YYERROR;
  712               }
  713            } else if (FSTRCMP($1,"ARCTAN2(") == 0) {
  714              if( TYPE($2) != DOUBLE ) $2 = New_Unary( DOUBLE, 0, $2 );
  715              if( TYPE($4) != DOUBLE ) $4 = New_Unary( DOUBLE, 0, $4 );
  716              if( Test_Dims( $2, $4 ) ) {
  717             $$ = New_Func( 0, atan2_fct, 2, $2, $4, 0, 0, 0, 0, 0 );
  718             TEST($$); 
  719             if( SIZE($2)<SIZE($4) ) Copy_Dims($$, $4);
  720              } else {
  721             yyerror("Dimensions of arctan2 arguments "
  722                 "are not compatible");
  723             YYERROR;
  724              }
  725            } else if (FSTRCMP($1,"MIN(") == 0) {
  726               PROMOTE( $2, $4 );
  727               if( Test_Dims( $2, $4 ) ) {
  728             $$ = New_Func( 0, min2_fct, 2, $2, $4, 0, 0, 0, 0, 0 );
  729             TEST($$);
  730             if( SIZE($2)<SIZE($4) ) Copy_Dims($$, $4);
  731               } else {
  732             yyerror("Dimensions of min(a,b) arguments "
  733                 "are not compatible");
  734             YYERROR;
  735               }
  736            } else if (FSTRCMP($1,"MAX(") == 0) {
  737               PROMOTE( $2, $4 );
  738               if( Test_Dims( $2, $4 ) ) {
  739             $$ = New_Func( 0, max2_fct, 2, $2, $4, 0, 0, 0, 0, 0 );
  740             TEST($$);
  741             if( SIZE($2)<SIZE($4) ) Copy_Dims($$, $4);
  742               } else {
  743             yyerror("Dimensions of max(a,b) arguments "
  744                 "are not compatible");
  745             YYERROR;
  746               }
  747            } else if (FSTRCMP($1,"SETNULL(") == 0) {
  748              if (OPER($2) != CONST_OP
  749              || SIZE($2) != 1) {
  750                yyerror("SETNULL first argument must be a scalar constant");
  751                YYERROR;
  752              }
  753              /* Make sure first arg is same type as second arg */
  754              if ( TYPE($2) != TYPE($4) ) $2 = New_Unary( TYPE($4), 0, $2 );
  755              $$ = New_Func( 0, setnull_fct, 2, $4, $2, 0, 0, 0, 0, 0 );
  756            } else {
  757               yyerror("Function(expr,expr) not supported");
  758               YYERROR;
  759            }
  760                 }
  761        | FUNCTION expr ',' expr ',' expr ',' expr ')'
  762                 { 
  763           if (FSTRCMP($1,"ANGSEP(") == 0) {
  764             if( TYPE($2) != DOUBLE ) $2 = New_Unary( DOUBLE, 0, $2 );
  765             if( TYPE($4) != DOUBLE ) $4 = New_Unary( DOUBLE, 0, $4 );
  766             if( TYPE($6) != DOUBLE ) $6 = New_Unary( DOUBLE, 0, $6 );
  767             if( TYPE($8) != DOUBLE ) $8 = New_Unary( DOUBLE, 0, $8 );
  768             if( Test_Dims( $2, $4 ) && Test_Dims( $4, $6 ) && 
  769             Test_Dims( $6, $8 ) ) {
  770               $$ = New_Func( 0, angsep_fct, 4, $2, $4, $6, $8,0,0,0 );
  771               TEST($$); 
  772               if( SIZE($2)<SIZE($4) ) Copy_Dims($$, $4);
  773               if( SIZE($4)<SIZE($6) ) Copy_Dims($$, $6);
  774               if( SIZE($6)<SIZE($8) ) Copy_Dims($$, $8);
  775             } else {
  776               yyerror("Dimensions of ANGSEP arguments "
  777                   "are not compatible");
  778               YYERROR;
  779             }
  780            } else {
  781               yyerror("Function(expr,expr,expr,expr) not supported");
  782               YYERROR;
  783            }
  784                 }
  785        | expr '[' expr ']'
  786                 { $$ = New_Deref( $1, 1, $3,  0,  0,  0,   0 ); TEST($$); }
  787        | expr '[' expr ',' expr ']'
  788                 { $$ = New_Deref( $1, 2, $3, $5,  0,  0,   0 ); TEST($$); }
  789        | expr '[' expr ',' expr ',' expr ']'
  790                 { $$ = New_Deref( $1, 3, $3, $5, $7,  0,   0 ); TEST($$); }
  791        | expr '[' expr ',' expr ',' expr ',' expr ']'
  792                 { $$ = New_Deref( $1, 4, $3, $5, $7, $9,   0 ); TEST($$); }
  793        | expr '[' expr ',' expr ',' expr ',' expr ',' expr ']'
  794                 { $$ = New_Deref( $1, 5, $3, $5, $7, $9, $11 ); TEST($$); }
  795        | INTCAST expr
  796         { $$ = New_Unary( LONG,   INTCAST, $2 );  TEST($$);  }
  797        | INTCAST bexpr
  798                 { $$ = New_Unary( LONG,   INTCAST, $2 );  TEST($$);  }
  799        | FLTCAST expr
  800         { $$ = New_Unary( DOUBLE, FLTCAST, $2 );  TEST($$);  }
  801        | FLTCAST bexpr
  802                 { $$ = New_Unary( DOUBLE, FLTCAST, $2 );  TEST($$);  }
  803        ;
  804 
  805 bexpr:   BOOLEAN
  806                 { $$ = New_Const( BOOLEAN, &($1), sizeof(char) ); TEST($$); }
  807        | BCOLUMN
  808                 { $$ = New_Column( $1 ); TEST($$); }
  809        | BCOLUMN '{' expr '}'
  810                 {
  811                   if( TYPE($3) != LONG
  812               || OPER($3) != CONST_OP ) {
  813              yyerror("Offset argument must be a constant integer");
  814              YYERROR;
  815           }
  816                   $$ = New_Offset( $1, $3 ); TEST($$);
  817                 }
  818        | bits EQ bits
  819                 { $$ = New_BinOp( BOOLEAN, $1, EQ,  $3 ); TEST($$);
  820           SIZE($$) = 1;                                     }
  821        | bits NE bits
  822                 { $$ = New_BinOp( BOOLEAN, $1, NE,  $3 ); TEST($$); 
  823           SIZE($$) = 1;                                     }
  824        | bits LT bits
  825                 { $$ = New_BinOp( BOOLEAN, $1, LT,  $3 ); TEST($$); 
  826           SIZE($$) = 1;                                     }
  827        | bits LTE bits
  828                 { $$ = New_BinOp( BOOLEAN, $1, LTE, $3 ); TEST($$); 
  829           SIZE($$) = 1;                                     }
  830        | bits GT bits
  831                 { $$ = New_BinOp( BOOLEAN, $1, GT,  $3 ); TEST($$); 
  832           SIZE($$) = 1;                                     }
  833        | bits GTE bits
  834                 { $$ = New_BinOp( BOOLEAN, $1, GTE, $3 ); TEST($$); 
  835           SIZE($$) = 1;                                     }
  836        | expr GT expr
  837                 { PROMOTE($1,$3); $$ = New_BinOp( BOOLEAN, $1, GT,  $3 );
  838                   TEST($$);                                               }
  839        | expr LT expr
  840                 { PROMOTE($1,$3); $$ = New_BinOp( BOOLEAN, $1, LT,  $3 );
  841                   TEST($$);                                               }
  842        | expr GTE expr
  843                 { PROMOTE($1,$3); $$ = New_BinOp( BOOLEAN, $1, GTE, $3 );
  844                   TEST($$);                                               }
  845        | expr LTE expr
  846                 { PROMOTE($1,$3); $$ = New_BinOp( BOOLEAN, $1, LTE, $3 );
  847                   TEST($$);                                               }
  848        | expr '~' expr
  849                 { PROMOTE($1,$3); $$ = New_BinOp( BOOLEAN, $1, '~', $3 );
  850                   TEST($$);                                               }
  851        | expr EQ expr
  852                 { PROMOTE($1,$3); $$ = New_BinOp( BOOLEAN, $1, EQ,  $3 );
  853                   TEST($$);                                               }
  854        | expr NE expr
  855                 { PROMOTE($1,$3); $$ = New_BinOp( BOOLEAN, $1, NE,  $3 );
  856                   TEST($$);                                               }
  857        | sexpr EQ sexpr
  858                 { $$ = New_BinOp( BOOLEAN, $1, EQ,  $3 ); TEST($$);
  859                   SIZE($$) = 1; }
  860        | sexpr NE sexpr
  861                 { $$ = New_BinOp( BOOLEAN, $1, NE,  $3 ); TEST($$);
  862                   SIZE($$) = 1; }
  863        | sexpr GT sexpr
  864                 { $$ = New_BinOp( BOOLEAN, $1, GT,  $3 ); TEST($$);
  865                   SIZE($$) = 1; }
  866        | sexpr GTE sexpr
  867                 { $$ = New_BinOp( BOOLEAN, $1, GTE, $3 ); TEST($$);
  868                   SIZE($$) = 1; }
  869        | sexpr LT sexpr
  870                 { $$ = New_BinOp( BOOLEAN, $1, LT,  $3 ); TEST($$);
  871                   SIZE($$) = 1; }
  872        | sexpr LTE sexpr
  873                 { $$ = New_BinOp( BOOLEAN, $1, LTE, $3 ); TEST($$);
  874                   SIZE($$) = 1; }
  875        | bexpr AND bexpr
  876                 { $$ = New_BinOp( BOOLEAN, $1, AND, $3 ); TEST($$); }
  877        | bexpr OR bexpr
  878                 { $$ = New_BinOp( BOOLEAN, $1, OR,  $3 ); TEST($$); }
  879        | bexpr EQ bexpr
  880                 { $$ = New_BinOp( BOOLEAN, $1, EQ,  $3 ); TEST($$); }
  881        | bexpr NE bexpr
  882                 { $$ = New_BinOp( BOOLEAN, $1, NE,  $3 ); TEST($$); }
  883 
  884        | expr '=' expr ':' expr
  885                 { PROMOTE($1,$3); PROMOTE($1,$5); PROMOTE($3,$5);
  886           $3 = New_BinOp( BOOLEAN, $3, LTE, $1 );
  887                   $5 = New_BinOp( BOOLEAN, $1, LTE, $5 );
  888                   $$ = New_BinOp( BOOLEAN, $3, AND, $5 );
  889                   TEST($$);                                         }
  890 
  891        | bexpr '?' bexpr ':' bexpr
  892                 {
  893                   if( ! Test_Dims($3,$5) ) {
  894                      yyerror("Incompatible dimensions in '?:' arguments");
  895              YYERROR;
  896                   }
  897                   $$ = New_Func( 0, ifthenelse_fct, 3, $3, $5, $1,
  898                                  0, 0, 0, 0 );
  899                   TEST($$);
  900                   if( SIZE($3)<SIZE($5) )  Copy_Dims($$, $5);
  901                   if( ! Test_Dims($1,$$) ) {
  902                      yyerror("Incompatible dimensions in '?:' condition");
  903              YYERROR;
  904                   }
  905                   if( SIZE($$)<SIZE($1) )  Copy_Dims($$, $1);
  906                 }
  907 
  908        | BFUNCTION expr ')'
  909                 {
  910            if (FSTRCMP($1,"ISNULL(") == 0) {
  911               $$ = New_Func( 0, isnull_fct, 1, $2, 0, 0,
  912                      0, 0, 0, 0 );
  913               TEST($$); 
  914                       /* Use expression's size, but return BOOLEAN */
  915               TYPE($$) = BOOLEAN;
  916            } else {
  917               yyerror("Boolean Function(expr) not supported");
  918               YYERROR;
  919            }
  920         }
  921        | BFUNCTION bexpr ')'
  922                 {
  923            if (FSTRCMP($1,"ISNULL(") == 0) {
  924               $$ = New_Func( 0, isnull_fct, 1, $2, 0, 0,
  925                      0, 0, 0, 0 );
  926               TEST($$); 
  927                       /* Use expression's size, but return BOOLEAN */
  928               TYPE($$) = BOOLEAN;
  929            } else {
  930               yyerror("Boolean Function(expr) not supported");
  931               YYERROR;
  932            }
  933         }
  934        | BFUNCTION sexpr ')'
  935                 {
  936            if (FSTRCMP($1,"ISNULL(") == 0) {
  937               $$ = New_Func( BOOLEAN, isnull_fct, 1, $2, 0, 0,
  938                      0, 0, 0, 0 );
  939               TEST($$); 
  940            } else {
  941               yyerror("Boolean Function(expr) not supported");
  942               YYERROR;
  943            }
  944         }
  945        | FUNCTION bexpr ',' bexpr ')'
  946                 {
  947            if (FSTRCMP($1,"DEFNULL(") == 0) {
  948               if( SIZE($2)>=SIZE($4) && Test_Dims( $2, $4 ) ) {
  949              $$ = New_Func( 0, defnull_fct, 2, $2, $4, 0,
  950                     0, 0, 0, 0 );
  951              TEST($$); 
  952               } else {
  953              yyerror("Dimensions of DEFNULL arguments are not compatible");
  954              YYERROR;
  955               }
  956            } else {
  957               yyerror("Boolean Function(expr,expr) not supported");
  958               YYERROR;
  959            }
  960         }
  961        | BFUNCTION expr ',' expr ',' expr ')'
  962         {
  963            if( TYPE($2) != DOUBLE ) $2 = New_Unary( DOUBLE, 0, $2 );
  964            if( TYPE($4) != DOUBLE ) $4 = New_Unary( DOUBLE, 0, $4 );
  965            if( TYPE($6) != DOUBLE ) $6 = New_Unary( DOUBLE, 0, $6 );
  966            if( ! (Test_Dims( $2, $4 ) && Test_Dims( $4, $6 ) ) ) {
  967                yyerror("Dimensions of NEAR arguments "
  968                    "are not compatible");
  969                YYERROR;
  970            } else {
  971              if (FSTRCMP($1,"NEAR(") == 0) {
  972                $$ = New_Func( BOOLEAN, near_fct, 3, $2, $4, $6,
  973                       0, 0, 0, 0 );
  974              } else {
  975                yyerror("Boolean Function not supported");
  976                YYERROR;
  977              }
  978              TEST($$); 
  979 
  980              if( SIZE($$)<SIZE($2) )  Copy_Dims($$, $2);
  981              if( SIZE($2)<SIZE($4) )  Copy_Dims($$, $4);
  982              if( SIZE($4)<SIZE($6) )  Copy_Dims($$, $6);
  983            }
  984         }
  985        | BFUNCTION expr ',' expr ',' expr ',' expr ',' expr ')'
  986             {
  987            if( TYPE($2) != DOUBLE ) $2 = New_Unary( DOUBLE, 0, $2 );
  988            if( TYPE($4) != DOUBLE ) $4 = New_Unary( DOUBLE, 0, $4 );
  989            if( TYPE($6) != DOUBLE ) $6 = New_Unary( DOUBLE, 0, $6 );
  990            if( TYPE($8) != DOUBLE ) $8 = New_Unary( DOUBLE, 0, $8 );
  991            if( TYPE($10)!= DOUBLE ) $10= New_Unary( DOUBLE, 0, $10);
  992            if( ! (Test_Dims( $2, $4 ) && Test_Dims( $4, $6 ) && 
  993               Test_Dims( $6, $8 ) && Test_Dims( $8, $10 )) ) {
  994              yyerror("Dimensions of CIRCLE arguments "
  995                  "are not compatible");
  996              YYERROR;
  997            } else {
  998              if (FSTRCMP($1,"CIRCLE(") == 0) {
  999                $$ = New_Func( BOOLEAN, circle_fct, 5, $2, $4, $6, $8,
 1000                       $10, 0, 0 );
 1001              } else {
 1002                yyerror("Boolean Function not supported");
 1003                YYERROR;
 1004              }
 1005              TEST($$); 
 1006              if( SIZE($$)<SIZE($2) )  Copy_Dims($$, $2);
 1007              if( SIZE($2)<SIZE($4) )  Copy_Dims($$, $4);
 1008              if( SIZE($4)<SIZE($6) )  Copy_Dims($$, $6);
 1009              if( SIZE($6)<SIZE($8) )  Copy_Dims($$, $8);
 1010              if( SIZE($8)<SIZE($10) ) Copy_Dims($$, $10);
 1011            }
 1012         }
 1013        | BFUNCTION expr ',' expr ',' expr ',' expr ',' expr ',' expr ',' expr ')'
 1014                 {
 1015            if( TYPE($2) != DOUBLE ) $2 = New_Unary( DOUBLE, 0, $2 );
 1016            if( TYPE($4) != DOUBLE ) $4 = New_Unary( DOUBLE, 0, $4 );
 1017            if( TYPE($6) != DOUBLE ) $6 = New_Unary( DOUBLE, 0, $6 );
 1018            if( TYPE($8) != DOUBLE ) $8 = New_Unary( DOUBLE, 0, $8 );
 1019            if( TYPE($10)!= DOUBLE ) $10= New_Unary( DOUBLE, 0, $10);
 1020            if( TYPE($12)!= DOUBLE ) $12= New_Unary( DOUBLE, 0, $12);
 1021            if( TYPE($14)!= DOUBLE ) $14= New_Unary( DOUBLE, 0, $14);
 1022            if( ! (Test_Dims( $2, $4 ) && Test_Dims( $4, $6 ) && 
 1023               Test_Dims( $6, $8 ) && Test_Dims( $8, $10 ) &&
 1024               Test_Dims($10,$12 ) && Test_Dims($12, $14 ) ) ) {
 1025              yyerror("Dimensions of BOX or ELLIPSE arguments "
 1026                  "are not compatible");
 1027              YYERROR;
 1028            } else {
 1029              if (FSTRCMP($1,"BOX(") == 0) {
 1030                $$ = New_Func( BOOLEAN, box_fct, 7, $2, $4, $6, $8,
 1031                       $10, $12, $14 );
 1032              } else if (FSTRCMP($1,"ELLIPSE(") == 0) {
 1033                $$ = New_Func( BOOLEAN, elps_fct, 7, $2, $4, $6, $8,
 1034                       $10, $12, $14 );
 1035              } else {
 1036                yyerror("SAO Image Function not supported");
 1037                YYERROR;
 1038              }
 1039              TEST($$); 
 1040              if( SIZE($$)<SIZE($2) )  Copy_Dims($$, $2);
 1041              if( SIZE($2)<SIZE($4) )  Copy_Dims($$, $4);
 1042              if( SIZE($4)<SIZE($6) )  Copy_Dims($$, $6);
 1043              if( SIZE($6)<SIZE($8) )  Copy_Dims($$, $8);
 1044              if( SIZE($8)<SIZE($10) ) Copy_Dims($$, $10);
 1045              if( SIZE($10)<SIZE($12) ) Copy_Dims($$, $12);
 1046              if( SIZE($12)<SIZE($14) ) Copy_Dims($$, $14);
 1047            }
 1048         }
 1049 
 1050        | GTIFILTER ')'
 1051                 { /* Use defaults for all elements */
 1052            $$ = New_GTI(gtifilt_fct,  "", -99, -99, "*START*", "*STOP*" );
 1053                    TEST($$);                                        }
 1054        | GTIFILTER STRING ')'
 1055                 { /* Use defaults for all except filename */
 1056           $$ = New_GTI(gtifilt_fct,  $2, -99, -99, "*START*", "*STOP*" );
 1057                    TEST($$);                                        }
 1058        | GTIFILTER STRING ',' expr ')'
 1059                 {  $$ = New_GTI(gtifilt_fct,  $2, $4, -99, "*START*", "*STOP*" );
 1060                    TEST($$);                                        }
 1061        | GTIFILTER STRING ',' expr ',' STRING ',' STRING ')'
 1062                 {  $$ = New_GTI(gtifilt_fct,  $2, $4, -99, $6, $8 );
 1063                    TEST($$);                                        }
 1064 
 1065 
 1066        | GTIOVERLAP STRING ',' expr ',' expr ')'
 1067                 {  $$ = New_GTI(gtiover_fct,  $2, $4, $6, "*START*", "*STOP*");
 1068                    TEST($$);                                        }
 1069        | GTIOVERLAP STRING ',' expr ',' expr ',' STRING ',' STRING ')'
 1070                 {  $$ = New_GTI(gtiover_fct,  $2, $4, $6, $8, $10 );
 1071                    TEST($$);                                        }
 1072 
 1073 
 1074        | REGFILTER STRING ')'
 1075                 { /* Use defaults for all except filename */
 1076                    $$ = New_REG( $2, -99, -99, "" );
 1077                    TEST($$);                                        }
 1078        | REGFILTER STRING ',' expr ',' expr ')'
 1079                 {  $$ = New_REG( $2, $4, $6, "" );
 1080                    TEST($$);                                        }
 1081        | REGFILTER STRING ',' expr ',' expr ',' STRING ')'
 1082                 {  $$ = New_REG( $2, $4, $6, $8 );
 1083                    TEST($$);                                        }
 1084 
 1085        | bexpr '[' expr ']'
 1086                 { $$ = New_Deref( $1, 1, $3,  0,  0,  0,   0 ); TEST($$); }
 1087        | bexpr '[' expr ',' expr ']'
 1088                 { $$ = New_Deref( $1, 2, $3, $5,  0,  0,   0 ); TEST($$); }
 1089        | bexpr '[' expr ',' expr ',' expr ']'
 1090                 { $$ = New_Deref( $1, 3, $3, $5, $7,  0,   0 ); TEST($$); }
 1091        | bexpr '[' expr ',' expr ',' expr ',' expr ']'
 1092                 { $$ = New_Deref( $1, 4, $3, $5, $7, $9,   0 ); TEST($$); }
 1093        | bexpr '[' expr ',' expr ',' expr ',' expr ',' expr ']'
 1094                 { $$ = New_Deref( $1, 5, $3, $5, $7, $9, $11 ); TEST($$); }
 1095        | NOT bexpr
 1096                 { $$ = New_Unary( BOOLEAN, NOT, $2 ); TEST($$); }
 1097        | '(' bexpr ')'
 1098                 { $$ = $2; }
 1099        ;
 1100 
 1101 sexpr:   STRING
 1102                 { $$ = New_Const( STRING, $1, strlen($1)+1 ); TEST($$);
 1103                   SIZE($$) = strlen($1); }
 1104        | SCOLUMN
 1105                 { $$ = New_Column( $1 ); TEST($$); }
 1106        | SCOLUMN '{' expr '}'
 1107                 {
 1108                   if( TYPE($3) != LONG
 1109               || OPER($3) != CONST_OP ) {
 1110              yyerror("Offset argument must be a constant integer");
 1111              YYERROR;
 1112           }
 1113                   $$ = New_Offset( $1, $3 ); TEST($$);
 1114                 }
 1115        | SNULLREF
 1116                 { $$ = New_Func( STRING, null_fct, 0, 0, 0, 0, 0, 0, 0, 0 ); }
 1117        | '(' sexpr ')'
 1118                 { $$ = $2; }
 1119        | sexpr '+' sexpr
 1120                 { 
 1121           if (SIZE($1)+SIZE($3) >= MAX_STRLEN) {
 1122             yyerror("Combined string size exceeds " MAX_STRLEN_S " characters");
 1123             YYERROR;
 1124           }
 1125           $$ = New_BinOp( STRING, $1, '+', $3 );  TEST($$);
 1126           SIZE($$) = SIZE($1) + SIZE($3);
 1127         }
 1128        | bexpr '?' sexpr ':' sexpr
 1129                 {
 1130           int outSize;
 1131                   if( SIZE($1)!=1 ) {
 1132                      yyerror("Cannot have a vector string column");
 1133              YYERROR;
 1134                   }
 1135           /* Since the output can be calculated now, as a constant
 1136              scalar, we must precalculate the output size, in
 1137              order to avoid an overflow. */
 1138           outSize = SIZE($3);
 1139           if (SIZE($5) > outSize) outSize = SIZE($5);
 1140                   $$ = New_FuncSize( 0, ifthenelse_fct, 3, $3, $5, $1,
 1141                      0, 0, 0, 0, outSize);
 1142           
 1143                   TEST($$);
 1144                   if( SIZE($3)<SIZE($5) )  Copy_Dims($$, $5);
 1145                 }
 1146 
 1147        | FUNCTION sexpr ',' sexpr ')'
 1148                 { 
 1149           if (FSTRCMP($1,"DEFNULL(") == 0) {
 1150              int outSize;
 1151              /* Since the output can be calculated now, as a constant
 1152             scalar, we must precalculate the output size, in
 1153             order to avoid an overflow. */
 1154              outSize = SIZE($2);
 1155              if (SIZE($4) > outSize) outSize = SIZE($4);
 1156              
 1157              $$ = New_FuncSize( 0, defnull_fct, 2, $2, $4, 0,
 1158                     0, 0, 0, 0, outSize );
 1159              TEST($$); 
 1160              if( SIZE($4)>SIZE($2) ) SIZE($$) = SIZE($4);
 1161           } else {
 1162              yyerror("Function(string,string) not supported");
 1163              YYERROR;
 1164           }
 1165         }
 1166        | FUNCTION sexpr ',' expr ',' expr ')'
 1167                 { 
 1168           if (FSTRCMP($1,"STRMID(") == 0) {
 1169             int len;
 1170             if( TYPE($4) != LONG || SIZE($4) != 1 ||
 1171             TYPE($6) != LONG || SIZE($6) != 1) {
 1172               yyerror("When using STRMID(S,P,N), P and N must be integers (and not vector columns)");
 1173               YYERROR;
 1174             }
 1175             if (OPER($6) == CONST_OP) {
 1176               /* Constant value: use that directly */
 1177               len = (gParse.Nodes[$6].value.data.lng);
 1178             } else {
 1179               /* Variable value: use the maximum possible (from $2) */
 1180               len = SIZE($2);
 1181             }
 1182             if (len <= 0 || len >= MAX_STRLEN) {
 1183               yyerror("STRMID(S,P,N), N must be 1-" MAX_STRLEN_S);
 1184               YYERROR;
 1185             }
 1186             $$ = New_FuncSize( 0, strmid_fct, 3, $2, $4,$6,0,0,0,0,len);
 1187             TEST($$);
 1188           } else {
 1189              yyerror("Function(string,expr,expr) not supported");
 1190              YYERROR;
 1191           }
 1192         }
 1193 
 1194     ;
 1195 
 1196 %%
 1197 
 1198 /*************************************************************************/
 1199 /*  Start of "New" routines which build the expression Nodal structure   */
 1200 /*************************************************************************/
 1201 
 1202 static int Alloc_Node( void )
 1203 {
 1204                       /* Use this for allocation to guarantee *Nodes */
 1205    Node *newNodePtr;  /* survives on failure, making it still valid  */
 1206                       /* while working our way out of this error     */
 1207 
 1208    if( gParse.nNodes == gParse.nNodesAlloc ) {
 1209       if( gParse.Nodes ) {
 1210      gParse.nNodesAlloc += gParse.nNodesAlloc;
 1211      newNodePtr = (Node *)realloc( gParse.Nodes,
 1212                        sizeof(Node)*gParse.nNodesAlloc );
 1213       } else {
 1214      gParse.nNodesAlloc = 100;
 1215      newNodePtr = (Node *)malloc ( sizeof(Node)*gParse.nNodesAlloc );
 1216       }  
 1217 
 1218       if( newNodePtr ) {
 1219      gParse.Nodes = newNodePtr;
 1220       } else {
 1221      gParse.status = MEMORY_ALLOCATION;
 1222      return( -1 );
 1223       }
 1224    }
 1225 
 1226    return ( gParse.nNodes++ );
 1227 }
 1228 
 1229 static void Free_Last_Node( void )
 1230 {
 1231    if( gParse.nNodes ) gParse.nNodes--;
 1232 }
 1233 
 1234 static int New_Const( int returnType, void *value, long len )
 1235 {
 1236    Node *this;
 1237    int n;
 1238 
 1239    n = Alloc_Node();
 1240    if( n>=0 ) {
 1241       this             = gParse.Nodes + n;
 1242       this->operation  = CONST_OP;             /* Flag a constant */
 1243       this->DoOp       = NULL;
 1244       this->nSubNodes  = 0;
 1245       this->type       = returnType;
 1246       memcpy( &(this->value.data), value, len );
 1247       this->value.undef = NULL;
 1248       this->value.nelem = 1;
 1249       this->value.naxis = 1;
 1250       this->value.naxes[0] = 1;
 1251    }
 1252    return(n);
 1253 }
 1254 
 1255 static int New_Column( int ColNum )
 1256 {
 1257    Node *this;
 1258    int  n, i;
 1259 
 1260    n = Alloc_Node();
 1261    if( n>=0 ) {
 1262       this              = gParse.Nodes + n;
 1263       this->operation   = -ColNum;
 1264       this->DoOp        = NULL;
 1265       this->nSubNodes   = 0;
 1266       this->type        = gParse.varData[ColNum].type;
 1267       this->value.nelem = gParse.varData[ColNum].nelem;
 1268       this->value.naxis = gParse.varData[ColNum].naxis;
 1269       for( i=0; i<gParse.varData[ColNum].naxis; i++ )
 1270      this->value.naxes[i] = gParse.varData[ColNum].naxes[i];
 1271    }
 1272    return(n);
 1273 }
 1274 
 1275 static int New_Offset( int ColNum, int offsetNode )
 1276 {
 1277    Node *this;
 1278    int  n, i, colNode;
 1279 
 1280    colNode = New_Column( ColNum );
 1281    if( colNode<0 ) return(-1);
 1282 
 1283    n = Alloc_Node();
 1284    if( n>=0 ) {
 1285       this              = gParse.Nodes + n;
 1286       this->operation   = '{';
 1287       this->DoOp        = Do_Offset;
 1288       this->nSubNodes   = 2;
 1289       this->SubNodes[0] = colNode;
 1290       this->SubNodes[1] = offsetNode;
 1291       this->type        = gParse.varData[ColNum].type;
 1292       this->value.nelem = gParse.varData[ColNum].nelem;
 1293       this->value.naxis = gParse.varData[ColNum].naxis;
 1294       for( i=0; i<gParse.varData[ColNum].naxis; i++ )
 1295      this->value.naxes[i] = gParse.varData[ColNum].naxes[i];
 1296    }
 1297    return(n);
 1298 }
 1299 
 1300 static int New_Unary( int returnType, int Op, int Node1 )
 1301 {
 1302    Node *this, *that;
 1303    int  i,n;
 1304 
 1305    if( Node1<0 ) return(-1);
 1306    that = gParse.Nodes + Node1;
 1307 
 1308    if( !Op ) Op = returnType;
 1309 
 1310    if( (Op==DOUBLE || Op==FLTCAST) && that->type==DOUBLE  ) return( Node1 );
 1311    if( (Op==LONG   || Op==INTCAST) && that->type==LONG    ) return( Node1 );
 1312    if( (Op==BOOLEAN              ) && that->type==BOOLEAN ) return( Node1 );
 1313    
 1314    n = Alloc_Node();
 1315    if( n>=0 ) {
 1316       this              = gParse.Nodes + n;
 1317       this->operation   = Op;
 1318       this->DoOp        = Do_Unary;
 1319       this->nSubNodes   = 1;
 1320       this->SubNodes[0] = Node1;
 1321       this->type        = returnType;
 1322 
 1323       that              = gParse.Nodes + Node1; /* Reset in case .Nodes mv'd */
 1324       this->value.nelem = that->value.nelem;
 1325       this->value.naxis = that->value.naxis;
 1326       for( i=0; i<that->value.naxis; i++ )
 1327      this->value.naxes[i] = that->value.naxes[i];
 1328 
 1329       if( that->operation==CONST_OP ) this->DoOp( this );
 1330    }
 1331    return( n );
 1332 }
 1333 
 1334 static int New_BinOp( int returnType, int Node1, int Op, int Node2 )
 1335 {
 1336    Node *this,*that1,*that2;
 1337    int  n,i,constant;
 1338 
 1339    if( Node1<0 || Node2<0 ) return(-1);
 1340 
 1341    n = Alloc_Node();
 1342    if( n>=0 ) {
 1343       this             = gParse.Nodes + n;
 1344       this->operation  = Op;
 1345       this->nSubNodes  = 2;
 1346       this->SubNodes[0]= Node1;
 1347       this->SubNodes[1]= Node2;
 1348       this->type       = returnType;
 1349 
 1350       that1            = gParse.Nodes + Node1;
 1351       that2            = gParse.Nodes + Node2;
 1352       constant         = (that1->operation==CONST_OP
 1353                           && that2->operation==CONST_OP);
 1354       if( that1->type!=STRING && that1->type!=BITSTR )
 1355      if( !Test_Dims( Node1, Node2 ) ) {
 1356         Free_Last_Node();
 1357         yyerror("Array sizes/dims do not match for binary operator");
 1358         return(-1);
 1359      }
 1360       if( that1->value.nelem == 1 ) that1 = that2;
 1361 
 1362       this->value.nelem = that1->value.nelem;
 1363       this->value.naxis = that1->value.naxis;
 1364       for( i=0; i<that1->value.naxis; i++ )
 1365      this->value.naxes[i] = that1->value.naxes[i];
 1366 
 1367       if ( Op == ACCUM && that1->type == BITSTR ) {
 1368     /* ACCUM is rank-reducing on bit strings */
 1369     this->value.nelem = 1;
 1370     this->value.naxis = 1;
 1371     this->value.naxes[0] = 1;
 1372       }
 1373 
 1374       /*  Both subnodes should be of same time  */
 1375       switch( that1->type ) {
 1376       case BITSTR:  this->DoOp = Do_BinOp_bit;  break;
 1377       case STRING:  this->DoOp = Do_BinOp_str;  break;
 1378       case BOOLEAN: this->DoOp = Do_BinOp_log;  break;
 1379       case LONG:    this->DoOp = Do_BinOp_lng;  break;
 1380       case DOUBLE:  this->DoOp = Do_BinOp_dbl;  break;
 1381       }
 1382       if( constant ) this->DoOp( this );
 1383    }
 1384    return( n );
 1385 }
 1386 
 1387 static int New_Func( int returnType, funcOp Op, int nNodes,
 1388              int Node1, int Node2, int Node3, int Node4, 
 1389              int Node5, int Node6, int Node7 )
 1390 {
 1391   return New_FuncSize(returnType, Op, nNodes,
 1392               Node1, Node2, Node3, Node4, 
 1393               Node5, Node6, Node7, 0);
 1394 }
 1395 
 1396 static int New_FuncSize( int returnType, funcOp Op, int nNodes,
 1397              int Node1, int Node2, int Node3, int Node4, 
 1398              int Node5, int Node6, int Node7, int Size )
 1399 /* If returnType==0 , use Node1's type and vector sizes as returnType, */
 1400 /* else return a single value of type returnType                       */
 1401 {
 1402    Node *this, *that;
 1403    int  i,n,constant;
 1404 
 1405    if( Node1<0 || Node2<0 || Node3<0 || Node4<0 || 
 1406        Node5<0 || Node6<0 || Node7<0 ) return(-1);
 1407 
 1408    n = Alloc_Node();
 1409    if( n>=0 ) {
 1410       this              = gParse.Nodes + n;
 1411       this->operation   = (int)Op;
 1412       this->DoOp        = Do_Func;
 1413       this->nSubNodes   = nNodes;
 1414       this->SubNodes[0] = Node1;
 1415       this->SubNodes[1] = Node2;
 1416       this->SubNodes[2] = Node3;
 1417       this->SubNodes[3] = Node4;
 1418       this->SubNodes[4] = Node5;
 1419       this->SubNodes[5] = Node6;
 1420       this->SubNodes[6] = Node7;
 1421       i = constant = nNodes;    /* Functions with zero params are not const */
 1422       if (Op == poirnd_fct) constant = 0; /* Nor is Poisson deviate */
 1423 
 1424       while( i-- )
 1425     constant = ( constant && OPER(this->SubNodes[i]) == CONST_OP );
 1426       
 1427       if( returnType ) {
 1428      this->type           = returnType;
 1429      this->value.nelem    = 1;
 1430      this->value.naxis    = 1;
 1431      this->value.naxes[0] = 1;
 1432       } else {
 1433      that              = gParse.Nodes + Node1;
 1434      this->type        = that->type;
 1435      this->value.nelem = that->value.nelem;
 1436      this->value.naxis = that->value.naxis;
 1437      for( i=0; i<that->value.naxis; i++ )
 1438         this->value.naxes[i] = that->value.naxes[i];
 1439       }
 1440       /* Force explicit size before evaluating */
 1441       if (Size > 0) this->value.nelem = Size;
 1442 
 1443       if( constant ) this->DoOp( this );
 1444    }
 1445    return( n );
 1446 }
 1447 
 1448 static int New_Deref( int Var,  int nDim,
 1449               int Dim1, int Dim2, int Dim3, int Dim4, int Dim5 )
 1450 {
 1451    int n, idx, constant;
 1452    long elem=0;
 1453    Node *this, *theVar, *theDim[MAXDIMS];
 1454 
 1455    if( Var<0 || Dim1<0 || Dim2<0 || Dim3<0 || Dim4<0 || Dim5<0 ) return(-1);
 1456 
 1457    theVar = gParse.Nodes + Var;
 1458    if( theVar->operation==CONST_OP || theVar->value.nelem==1 ) {
 1459       yyerror("Cannot index a scalar value");
 1460       return(-1);
 1461    }
 1462 
 1463    n = Alloc_Node();
 1464    if( n>=0 ) {
 1465       this              = gParse.Nodes + n;
 1466       this->nSubNodes   = nDim+1;
 1467       theVar            = gParse.Nodes + (this->SubNodes[0]=Var);
 1468       theDim[0]         = gParse.Nodes + (this->SubNodes[1]=Dim1);
 1469       theDim[1]         = gParse.Nodes + (this->SubNodes[2]=Dim2);
 1470       theDim[2]         = gParse.Nodes + (this->SubNodes[3]=Dim3);
 1471       theDim[3]         = gParse.Nodes + (this->SubNodes[4]=Dim4);
 1472       theDim[4]         = gParse.Nodes + (this->SubNodes[5]=Dim5);
 1473       constant          = theVar->operation==CONST_OP;
 1474       for( idx=0; idx<nDim; idx++ )
 1475      constant = (constant && theDim[idx]->operation==CONST_OP);
 1476 
 1477       for( idx=0; idx<nDim; idx++ )
 1478      if( theDim[idx]->value.nelem>1 ) {
 1479         Free_Last_Node();
 1480         yyerror("Cannot use an array as an index value");
 1481         return(-1);
 1482      } else if( theDim[idx]->type!=LONG ) {
 1483         Free_Last_Node();
 1484         yyerror("Index value must be an integer type");
 1485         return(-1);
 1486      }
 1487 
 1488       this->operation   = '[';
 1489       this->DoOp        = Do_Deref;
 1490       this->type        = theVar->type;
 1491 
 1492       if( theVar->value.naxis == nDim ) { /* All dimensions specified */
 1493      this->value.nelem    = 1;
 1494      this->value.naxis    = 1;
 1495      this->value.naxes[0] = 1;
 1496       } else if( nDim==1 ) { /* Dereference only one dimension */
 1497      elem=1;
 1498      this->value.naxis = theVar->value.naxis-1;
 1499      for( idx=0; idx<this->value.naxis; idx++ ) {
 1500         elem *= ( this->value.naxes[idx] = theVar->value.naxes[idx] );
 1501      }
 1502      this->value.nelem = elem;
 1503       } else {
 1504      Free_Last_Node();
 1505      yyerror("Must specify just one or all indices for vector");
 1506      return(-1);
 1507       }
 1508       if( constant ) this->DoOp( this );
 1509    }
 1510    return(n);
 1511 }
 1512 
 1513 extern int yyGetVariable( char *varName, YYSTYPE *varVal );
 1514 
 1515 static int New_GTI( funcOp Op, char *fname, int Node1, int Node2, char *start, char *stop )
 1516 {
 1517    fitsfile *fptr;
 1518    Node *this, *that0, *that1, *that2;
 1519    int  type,i,n, startCol, stopCol, Node0;
 1520    int  hdutype, hdunum, evthdu, samefile, extvers, movetotype, tstat;
 1521    char extname[100];
 1522    long nrows;
 1523    double timeZeroI[2], timeZeroF[2], dt, timeSpan;
 1524    char xcol[20], xexpr[20];
 1525    YYSTYPE colVal;
 1526 
 1527    if( Op == gtifilt_fct && Node1==-99 ) {
 1528       type = yyGetVariable( "TIME", &colVal );
 1529       if( type==COLUMN ) {
 1530      Node1 = New_Column( (int)colVal.lng );
 1531       } else {
 1532      yyerror("Could not build TIME column for GTIFILTER");
 1533      return(-1);
 1534       }
 1535    }
 1536 
 1537    if (Op == gtiover_fct) {
 1538      if (Node1 == -99 || Node2 == -99) {
 1539        yyerror("startExpr and stopExpr values must be defined for GTIOVERLAP");
 1540        return(-1);
 1541      }
 1542      /* Also case TIME_STOP to double precision */
 1543      Node2 = New_Unary( DOUBLE, 0, Node2 );
 1544      if (Node2 < 0) return(-1);
 1545 
 1546    }
 1547 
 1548    /* Type cast TIME to double precision */
 1549    Node1 = New_Unary( DOUBLE, 0, Node1 );
 1550    Node0 = Alloc_Node(); /* This will hold the START/STOP times */
 1551    if( Node1<0 || Node0<0 ) return(-1);
 1552 
 1553    /*  Record current HDU number in case we need to move within this file  */
 1554 
 1555    fptr = gParse.def_fptr;
 1556    ffghdn( fptr, &evthdu );
 1557 
 1558    /*  Look for TIMEZERO keywords in current extension  */
 1559 
 1560    tstat = 0;
 1561    if( ffgkyd( fptr, "TIMEZERO", timeZeroI, NULL, &tstat ) ) {
 1562       tstat = 0;
 1563       if( ffgkyd( fptr, "TIMEZERI", timeZeroI, NULL, &tstat ) ) {
 1564      timeZeroI[0] = timeZeroF[0] = 0.0;
 1565       } else if( ffgkyd( fptr, "TIMEZERF", timeZeroF, NULL, &tstat ) ) {
 1566      timeZeroF[0] = 0.0;
 1567       }
 1568    } else {
 1569       timeZeroF[0] = 0.0;
 1570    }
 1571 
 1572    /*  Resolve filename parameter  */
 1573 
 1574    switch( fname[0] ) {
 1575    case '\0':
 1576       samefile = 1;
 1577       hdunum = 1;
 1578       break;
 1579    case '[':
 1580       samefile = 1;
 1581       i = 1;
 1582       while( fname[i] != '\0' && fname[i] != ']' ) i++;
 1583       if( fname[i] ) {
 1584      fname[i] = '\0';
 1585      fname++;
 1586      ffexts( fname, &hdunum, extname, &extvers, &movetotype,
 1587          xcol, xexpr, &gParse.status );
 1588          if( *extname ) {
 1589         ffmnhd( fptr, movetotype, extname, extvers, &gParse.status );
 1590         ffghdn( fptr, &hdunum );
 1591      } else if( hdunum ) {
 1592         ffmahd( fptr, ++hdunum, &hdutype, &gParse.status );
 1593      } else if( !gParse.status ) {
 1594         yyerror("Cannot use primary array for GTI filter");
 1595         return( -1 );
 1596      }
 1597       } else {
 1598      yyerror("File extension specifier lacks closing ']'");
 1599      return( -1 );
 1600       }
 1601       break;
 1602    case '+':
 1603       samefile = 1;
 1604       hdunum = atoi( fname ) + 1;
 1605       if( hdunum>1 )
 1606      ffmahd( fptr, hdunum, &hdutype, &gParse.status );
 1607       else {
 1608      yyerror("Cannot use primary array for GTI filter");
 1609      return( -1 );
 1610       }
 1611       break;
 1612    default:
 1613       samefile = 0;
 1614       if( ! ffopen( &fptr, fname, READONLY, &gParse.status ) )
 1615      ffghdn( fptr, &hdunum );
 1616       break;
 1617    }
 1618    if( gParse.status ) return(-1);
 1619 
 1620    /*  If at primary, search for GTI extension  */
 1621 
 1622    if( hdunum==1 ) {
 1623       while( 1 ) {
 1624      hdunum++;
 1625      if( ffmahd( fptr, hdunum, &hdutype, &gParse.status ) ) break;
 1626      if( hdutype==IMAGE_HDU ) continue;
 1627      tstat = 0;
 1628      if( ffgkys( fptr, "EXTNAME", extname, NULL, &tstat ) ) continue;
 1629      ffupch( extname );
 1630      if( strstr( extname, "GTI" ) ) break;
 1631       }
 1632       if( gParse.status ) {
 1633      if( gParse.status==END_OF_FILE )
 1634         yyerror("GTI extension not found in this file");
 1635      return(-1);
 1636       }
 1637    }
 1638 
 1639    /*  Locate START/STOP Columns  */
 1640 
 1641    ffgcno( fptr, CASEINSEN, start, &startCol, &gParse.status );
 1642    ffgcno( fptr, CASEINSEN, stop,  &stopCol,  &gParse.status );
 1643    if( gParse.status ) return(-1);
 1644 
 1645    /*  Look for TIMEZERO keywords in GTI extension  */
 1646 
 1647    tstat = 0;
 1648    if( ffgkyd( fptr, "TIMEZERO", timeZeroI+1, NULL, &tstat ) ) {
 1649       tstat = 0;
 1650       if( ffgkyd( fptr, "TIMEZERI", timeZeroI+1, NULL, &tstat ) ) {
 1651      timeZeroI[1] = timeZeroF[1] = 0.0;
 1652       } else if( ffgkyd( fptr, "TIMEZERF", timeZeroF+1, NULL, &tstat ) ) {
 1653      timeZeroF[1] = 0.0;
 1654       }
 1655    } else {
 1656       timeZeroF[1] = 0.0;
 1657    }
 1658 
 1659    n = Alloc_Node();
 1660    if( n >= 0 ) {
 1661       this                 = gParse.Nodes + n;
 1662       this->SubNodes[1]    = Node1;
 1663       this->operation      = (int) Op;
 1664       if (Op == gtifilt_fct) {
 1665     this->nSubNodes      = 2;
 1666     this->DoOp           = Do_GTI;
 1667     this->type           = BOOLEAN;
 1668       } else {
 1669     this->nSubNodes      = 3;
 1670     this->DoOp           = Do_GTI_Over;
 1671     this->type           = DOUBLE;
 1672       }
 1673       that1                = gParse.Nodes + Node1;
 1674       this->value.nelem    = that1->value.nelem;
 1675       this->value.naxis    = that1->value.naxis;
 1676       for( i=0; i < that1->value.naxis; i++ )
 1677      this->value.naxes[i] = that1->value.naxes[i];
 1678       if (Op == gtiover_fct) {
 1679     this->SubNodes[2]  = Node2;
 1680     that2 = gParse.Nodes + Node2;
 1681     if (that1->value.nelem != that2->value.nelem) {
 1682       yyerror("Dimensions of TIME and TIME_STOP must match for GTIOVERLAP");
 1683       return(-1);
 1684     }
 1685       }
 1686 
 1687       /* Init START/STOP node to be treated as a "constant" */
 1688 
 1689       this->SubNodes[0]    = Node0;
 1690       that0                = gParse.Nodes + Node0;
 1691       that0->operation     = CONST_OP;
 1692       that0->DoOp          = NULL;
 1693       that0->value.data.ptr= NULL;
 1694 
 1695       /*  Read in START/STOP times  */
 1696 
 1697       if( ffgkyj( fptr, "NAXIS2", &nrows, NULL, &gParse.status ) )
 1698      return(-1);
 1699       that0->value.nelem = nrows;
 1700       if( nrows ) {
 1701 
 1702      that0->value.data.dblptr = (double*)malloc( 2*nrows*sizeof(double) );
 1703      if( !that0->value.data.dblptr ) {
 1704         gParse.status = MEMORY_ALLOCATION;
 1705         return(-1);
 1706      }
 1707      
 1708      ffgcvd( fptr, startCol, 1L, 1L, nrows, 0.0,
 1709          that0->value.data.dblptr, &i, &gParse.status );
 1710      ffgcvd( fptr, stopCol, 1L, 1L, nrows, 0.0,
 1711          that0->value.data.dblptr+nrows, &i, &gParse.status );
 1712      if( gParse.status ) {
 1713         free( that0->value.data.dblptr );
 1714         return(-1);
 1715      }
 1716 
 1717      /*  Test for fully time-ordered GTI... both START && STOP  */
 1718 
 1719      that0->type = 1; /*  Assume yes  */
 1720      i = nrows;
 1721      while( --i )
 1722         if(    that0->value.data.dblptr[i-1]
 1723                    >= that0->value.data.dblptr[i]
 1724         || that0->value.data.dblptr[i-1+nrows]
 1725            >= that0->value.data.dblptr[i+nrows] ) {
 1726            that0->type = 0;
 1727            break;
 1728         }
 1729 
 1730      /* GTIOVERLAP() requires ordered GTI */
 1731      if (that0->type != 1 && Op == gtiover_fct) {
 1732        yyerror("Input GTI must be time-ordered for GTIOVERLAP");
 1733        return(-1);
 1734      }
 1735      
 1736      /*  Handle TIMEZERO offset, if any  */
 1737      
 1738      dt = (timeZeroI[1] - timeZeroI[0]) + (timeZeroF[1] - timeZeroF[0]);
 1739      timeSpan = that0->value.data.dblptr[nrows+nrows-1]
 1740         - that0->value.data.dblptr[0];
 1741      
 1742      if( fabs( dt / timeSpan ) > 1e-12 ) {
 1743         for( i=0; i<(nrows+nrows); i++ )
 1744            that0->value.data.dblptr[i] += dt;
 1745      }
 1746       }
 1747       /* If Node1 is constant (gtifilt_fct) or
 1748      Node1 and Node2 are constant (gtiover_fct), then evaluate now */
 1749       if( OPER(Node1)==CONST_OP && (Op == gtifilt_fct || OPER(Node2)==CONST_OP)) {
 1750     this->DoOp( this );
 1751       }
 1752    }
 1753 
 1754    if( samefile )
 1755       ffmahd( fptr, evthdu, &hdutype, &gParse.status );
 1756    else
 1757       ffclos( fptr, &gParse.status );
 1758 
 1759    return( n );
 1760 }
 1761 
 1762 static int New_REG( char *fname, int NodeX, int NodeY, char *colNames )
 1763 {
 1764    Node *this, *that0;
 1765    int  type, n, Node0;
 1766    int  Xcol, Ycol, tstat;
 1767    WCSdata wcs;
 1768    SAORegion *Rgn;
 1769    char *cX, *cY;
 1770    YYSTYPE colVal;
 1771 
 1772    if( NodeX==-99 ) {
 1773       type = yyGetVariable( "X", &colVal );
 1774       if( type==COLUMN ) {
 1775      NodeX = New_Column( (int)colVal.lng );
 1776       } else {
 1777      yyerror("Could not build X column for REGFILTER");
 1778      return(-1);
 1779       }
 1780    }
 1781    if( NodeY==-99 ) {
 1782       type = yyGetVariable( "Y", &colVal );
 1783       if( type==COLUMN ) {
 1784      NodeY = New_Column( (int)colVal.lng );
 1785       } else {
 1786      yyerror("Could not build Y column for REGFILTER");
 1787      return(-1);
 1788       }
 1789    }
 1790    NodeX = New_Unary( DOUBLE, 0, NodeX );
 1791    NodeY = New_Unary( DOUBLE, 0, NodeY );
 1792    Node0 = Alloc_Node(); /* This will hold the Region Data */
 1793    if( NodeX<0 || NodeY<0 || Node0<0 ) return(-1);
 1794 
 1795    if( ! (Test_Dims( NodeX, NodeY ) ) ) {
 1796      yyerror("Dimensions of REGFILTER arguments are not compatible");
 1797      return (-1);
 1798    }
 1799 
 1800    n = Alloc_Node();
 1801    if( n >= 0 ) {
 1802       this                 = gParse.Nodes + n;
 1803       this->nSubNodes      = 3;
 1804       this->SubNodes[0]    = Node0;
 1805       this->SubNodes[1]    = NodeX;
 1806       this->SubNodes[2]    = NodeY;
 1807       this->operation      = (int)regfilt_fct;
 1808       this->DoOp           = Do_REG;
 1809       this->type           = BOOLEAN;
 1810       this->value.nelem    = 1;
 1811       this->value.naxis    = 1;
 1812       this->value.naxes[0] = 1;
 1813       
 1814       Copy_Dims(n, NodeX);
 1815       if( SIZE(NodeX)<SIZE(NodeY) )  Copy_Dims(n, NodeY);
 1816 
 1817       /* Init Region node to be treated as a "constant" */
 1818 
 1819       that0                = gParse.Nodes + Node0;
 1820       that0->operation     = CONST_OP;
 1821       that0->DoOp          = NULL;
 1822 
 1823       /*  Identify what columns to use for WCS information  */
 1824 
 1825       Xcol = Ycol = 0;
 1826       if( *colNames ) {
 1827      /*  Use the column names in this string for WCS info  */
 1828      while( *colNames==' ' ) colNames++;
 1829      cX = cY = colNames;
 1830      while( *cY && *cY!=' ' && *cY!=',' ) cY++;
 1831      if( *cY )
 1832         *(cY++) = '\0';
 1833      while( *cY==' ' ) cY++;
 1834      if( !*cY ) {
 1835         yyerror("Could not extract valid pair of column names from REGFILTER");
 1836         Free_Last_Node();
 1837         return( -1 );
 1838      }
 1839      fits_get_colnum( gParse.def_fptr, CASEINSEN, cX, &Xcol,
 1840               &gParse.status );
 1841      fits_get_colnum( gParse.def_fptr, CASEINSEN, cY, &Ycol,
 1842               &gParse.status );
 1843      if( gParse.status ) {
 1844         yyerror("Could not locate columns indicated for WCS info");
 1845         Free_Last_Node();
 1846         return( -1 );
 1847      }
 1848 
 1849       } else {
 1850      /*  Try to find columns used in X/Y expressions  */
 1851      Xcol = Locate_Col( gParse.Nodes + NodeX );
 1852      Ycol = Locate_Col( gParse.Nodes + NodeY );
 1853      if( Xcol<0 || Ycol<0 ) {
 1854         yyerror("Found multiple X/Y column references in REGFILTER");
 1855         Free_Last_Node();
 1856         return( -1 );
 1857      }
 1858       }
 1859 
 1860       /*  Now, get the WCS info, if it exists, from the indicated columns  */
 1861       wcs.exists = 0;
 1862       if( Xcol>0 && Ycol>0 ) {
 1863      tstat = 0;
 1864      ffgtcs( gParse.def_fptr, Xcol, Ycol,
 1865          &wcs.xrefval, &wcs.yrefval,
 1866          &wcs.xrefpix, &wcs.yrefpix,
 1867          &wcs.xinc,    &wcs.yinc,
 1868          &wcs.rot,      wcs.type,
 1869          &tstat );
 1870      if( tstat==NO_WCS_KEY ) {
 1871         wcs.exists = 0;
 1872      } else if( tstat ) {
 1873         gParse.status = tstat;
 1874         Free_Last_Node();
 1875         return( -1 );
 1876      } else {
 1877         wcs.exists = 1;
 1878      }
 1879       }
 1880 
 1881       /*  Read in Region file  */
 1882 
 1883       fits_read_rgnfile( fname, &wcs, &Rgn, &gParse.status );
 1884       if( gParse.status ) {
 1885      Free_Last_Node();
 1886      return( -1 );
 1887       }
 1888 
 1889       that0->value.data.ptr = Rgn;
 1890 
 1891       if( OPER(NodeX)==CONST_OP && OPER(NodeY)==CONST_OP )
 1892      this->DoOp( this );
 1893    }
 1894 
 1895    return( n );
 1896 }
 1897 
 1898 static int New_Vector( int subNode )
 1899 {
 1900    Node *this, *that;
 1901    int n;
 1902 
 1903    n = Alloc_Node();
 1904    if( n >= 0 ) {
 1905       this              = gParse.Nodes + n;
 1906       that              = gParse.Nodes + subNode;
 1907       this->type        = that->type;
 1908       this->nSubNodes   = 1;
 1909       this->SubNodes[0] = subNode;
 1910       this->operation   = '{';
 1911       this->DoOp        = Do_Vector;
 1912    }
 1913 
 1914    return( n );
 1915 }
 1916 
 1917 static int Close_Vec( int vecNode )
 1918 {
 1919    Node *this;
 1920    int n, nelem=0;
 1921 
 1922    this = gParse.Nodes + vecNode;
 1923    for( n=0; n < this->nSubNodes; n++ ) {
 1924       if( TYPE( this->SubNodes[n] ) != this->type ) {
 1925      this->SubNodes[n] = New_Unary( this->type, 0, this->SubNodes[n] );
 1926      if( this->SubNodes[n]<0 ) return(-1);
 1927       }
 1928       nelem += SIZE(this->SubNodes[n]);
 1929    }
 1930    this->value.naxis    = 1;
 1931    this->value.nelem    = nelem;
 1932    this->value.naxes[0] = nelem;
 1933 
 1934    return( vecNode );
 1935 }
 1936 
 1937 static int Locate_Col( Node *this )
 1938 /*  Locate the TABLE column number of any columns in "this" calculation.  */
 1939 /*  Return ZERO if none found, or negative if more than 1 found.          */
 1940 {
 1941    Node *that;
 1942    int  i, col=0, newCol, nfound=0;
 1943    
 1944    if( this->nSubNodes==0
 1945        && this->operation<=0 && this->operation!=CONST_OP )
 1946       return gParse.colData[ - this->operation].colnum;
 1947 
 1948    for( i=0; i<this->nSubNodes; i++ ) {
 1949       that = gParse.Nodes + this->SubNodes[i];
 1950       if( that->operation>0 ) {
 1951      newCol = Locate_Col( that );
 1952      if( newCol<=0 ) {
 1953         nfound += -newCol;
 1954      } else {
 1955         if( !nfound ) {
 1956            col = newCol;
 1957            nfound++;
 1958         } else if( col != newCol ) {
 1959            nfound++;
 1960         }
 1961      }
 1962       } else if( that->operation!=CONST_OP ) {
 1963      /*  Found a Column  */
 1964      newCol = gParse.colData[- that->operation].colnum;
 1965      if( !nfound ) {
 1966         col = newCol;
 1967         nfound++;
 1968      } else if( col != newCol ) {
 1969         nfound++;
 1970      }
 1971       }
 1972    }
 1973    if( nfound!=1 )
 1974       return( - nfound );
 1975    else
 1976       return( col );
 1977 }
 1978 
 1979 static int Test_Dims( int Node1, int Node2 )
 1980 {
 1981    Node *that1, *that2;
 1982    int valid, i;
 1983 
 1984    if( Node1<0 || Node2<0 ) return(0);
 1985 
 1986    that1 = gParse.Nodes + Node1;
 1987    that2 = gParse.Nodes + Node2;
 1988 
 1989    if( that1->value.nelem==1 || that2->value.nelem==1 )
 1990       valid = 1;
 1991    else if( that1->type==that2->type
 1992         && that1->value.nelem==that2->value.nelem
 1993         && that1->value.naxis==that2->value.naxis ) {
 1994       valid = 1;
 1995       for( i=0; i<that1->value.naxis; i++ ) {
 1996      if( that1->value.naxes[i]!=that2->value.naxes[i] )
 1997         valid = 0;
 1998       }
 1999    } else
 2000       valid = 0;
 2001    return( valid );
 2002 }   
 2003 
 2004 static void Copy_Dims( int Node1, int Node2 )
 2005 {
 2006    Node *that1, *that2;
 2007    int i;
 2008 
 2009    if( Node1<0 || Node2<0 ) return;
 2010 
 2011    that1 = gParse.Nodes + Node1;
 2012    that2 = gParse.Nodes + Node2;
 2013 
 2014    that1->value.nelem = that2->value.nelem;
 2015    that1->value.naxis = that2->value.naxis;
 2016    for( i=0; i<that2->value.naxis; i++ )
 2017       that1->value.naxes[i] = that2->value.naxes[i];
 2018 }
 2019 
 2020 /********************************************************************/
 2021 /*    Routines for actually evaluating the expression start here    */
 2022 /********************************************************************/
 2023 
 2024 void Evaluate_Parser( long firstRow, long nRows )
 2025     /***********************************************************************/
 2026     /*  Reset the parser for processing another batch of data...           */
 2027     /*    firstRow:  Row number of the first element to evaluate           */
 2028     /*    nRows:     Number of rows to be processed                        */
 2029     /*  Initialize each COLUMN node so that its UNDEF and DATA pointers    */
 2030     /*  point to the appropriate column arrays.                            */
 2031     /*  Finally, call Evaluate_Node for final node.                        */
 2032     /***********************************************************************/
 2033 {
 2034    int     i, column;
 2035    long    offset, rowOffset;
 2036    static int rand_initialized = 0;
 2037 
 2038    /* Initialize the random number generator once and only once */
 2039    if (rand_initialized == 0) {
 2040      simplerng_srand( (unsigned int) time(NULL) );
 2041      rand_initialized = 1;
 2042    }
 2043 
 2044    gParse.firstRow = firstRow;
 2045    gParse.nRows    = nRows;
 2046 
 2047    /*  Reset Column Nodes' pointers to point to right data and UNDEF arrays  */
 2048 
 2049    rowOffset = firstRow - gParse.firstDataRow;
 2050    for( i=0; i<gParse.nNodes; i++ ) {
 2051      if(    OPER(i) >  0 || OPER(i) == CONST_OP ) continue;
 2052 
 2053       column = -OPER(i);
 2054       offset = gParse.varData[column].nelem * rowOffset;
 2055 
 2056       gParse.Nodes[i].value.undef = gParse.varData[column].undef + offset;
 2057 
 2058       switch( gParse.Nodes[i].type ) {
 2059       case BITSTR:
 2060      gParse.Nodes[i].value.data.strptr =
 2061         (char**)gParse.varData[column].data + rowOffset;
 2062      gParse.Nodes[i].value.undef       = NULL;
 2063      break;
 2064       case STRING:
 2065      gParse.Nodes[i].value.data.strptr = 
 2066         (char**)gParse.varData[column].data + rowOffset;
 2067      gParse.Nodes[i].value.undef = gParse.varData[column].undef + rowOffset;
 2068      break;
 2069       case BOOLEAN:
 2070      gParse.Nodes[i].value.data.logptr = 
 2071         (char*)gParse.varData[column].data + offset;
 2072      break;
 2073       case LONG:
 2074      gParse.Nodes[i].value.data.lngptr = 
 2075         (long*)gParse.varData[column].data + offset;
 2076      break;
 2077       case DOUBLE:
 2078      gParse.Nodes[i].value.data.dblptr = 
 2079         (double*)gParse.varData[column].data + offset;
 2080      break;
 2081       }
 2082    }
 2083 
 2084    Evaluate_Node( gParse.resultNode );
 2085 }
 2086 
 2087 static void Evaluate_Node( int thisNode )
 2088     /**********************************************************************/
 2089     /*  Recursively evaluate thisNode's subNodes, then call one of the    */
 2090     /*  Do_<Action> functions pointed to by thisNode's DoOp element.      */
 2091     /**********************************************************************/
 2092 {
 2093    Node *this;
 2094    int i;
 2095    
 2096    if( gParse.status ) return;
 2097 
 2098    this = gParse.Nodes + thisNode;
 2099    if( this->operation>0 ) {  /* <=0 indicate constants and columns */
 2100       i = this->nSubNodes;
 2101       while( i-- ) {
 2102      Evaluate_Node( this->SubNodes[i] );
 2103      if( gParse.status ) return;
 2104       }
 2105       this->DoOp( this );
 2106    }
 2107 }
 2108 
 2109 static void Allocate_Ptrs( Node *this )
 2110 {
 2111    long elem, row, size;
 2112 
 2113    if( this->type==BITSTR || this->type==STRING ) {
 2114 
 2115       this->value.data.strptr = (char**)malloc( gParse.nRows
 2116                         * sizeof(char*) );
 2117       if( this->value.data.strptr ) {
 2118      this->value.data.strptr[0] = (char*)malloc( gParse.nRows
 2119                              * (this->value.nelem+2)
 2120                              * sizeof(char) );
 2121      if( this->value.data.strptr[0] ) {
 2122         row = 0;
 2123         while( (++row)<gParse.nRows ) {
 2124            this->value.data.strptr[row] =
 2125           this->value.data.strptr[row-1] + this->value.nelem+1;
 2126         }
 2127         if( this->type==STRING ) {
 2128            this->value.undef = this->value.data.strptr[row-1]
 2129                                    + this->value.nelem+1;
 2130         } else {
 2131            this->value.undef = NULL;  /* BITSTRs don't use undef array */
 2132         }
 2133      } else {
 2134         gParse.status = MEMORY_ALLOCATION;
 2135         free( this->value.data.strptr );
 2136      }
 2137       } else {
 2138      gParse.status = MEMORY_ALLOCATION;
 2139       }
 2140 
 2141    } else {
 2142 
 2143       elem = this->value.nelem * gParse.nRows;
 2144       switch( this->type ) {
 2145       case DOUBLE:  size = sizeof( double ); break;
 2146       case LONG:    size = sizeof( long   ); break;
 2147       case BOOLEAN: size = sizeof( char   ); break;
 2148       default:      size = 1;                break;
 2149       }
 2150 
 2151       this->value.data.ptr = calloc(size+1, elem);
 2152 
 2153       if( this->value.data.ptr==NULL ) {
 2154      gParse.status = MEMORY_ALLOCATION;
 2155       } else {
 2156      this->value.undef = (char *)this->value.data.ptr + elem*size;
 2157       }
 2158    }
 2159 }
 2160 
 2161 static void Do_Unary( Node *this )
 2162 {
 2163    Node *that;
 2164    long elem;
 2165 
 2166    that = gParse.Nodes + this->SubNodes[0];
 2167 
 2168    if( that->operation==CONST_OP ) {  /* Operating on a constant! */
 2169       switch( this->operation ) {
 2170       case DOUBLE:
 2171       case FLTCAST:
 2172      if( that->type==LONG )
 2173         this->value.data.dbl = (double)that->value.data.lng;
 2174      else if( that->type==BOOLEAN )
 2175         this->value.data.dbl = ( that->value.data.log ? 1.0 : 0.0 );
 2176      break;
 2177       case LONG:
 2178       case INTCAST:
 2179      if( that->type==DOUBLE )
 2180         this->value.data.lng = (long)that->value.data.dbl;
 2181      else if( that->type==BOOLEAN )
 2182         this->value.data.lng = ( that->value.data.log ? 1L : 0L );
 2183      break;
 2184       case BOOLEAN:
 2185      if( that->type==DOUBLE )
 2186         this->value.data.log = ( that->value.data.dbl != 0.0 );
 2187      else if( that->type==LONG )
 2188         this->value.data.log = ( that->value.data.lng != 0L );
 2189      break;
 2190       case UMINUS:
 2191      if( that->type==DOUBLE )
 2192         this->value.data.dbl = - that->value.data.dbl;
 2193      else if( that->type==LONG )
 2194         this->value.data.lng = - that->value.data.lng;
 2195      break;
 2196       case NOT:
 2197      if( that->type==BOOLEAN )
 2198         this->value.data.log = ( ! that->value.data.log );
 2199      else if( that->type==BITSTR )
 2200         bitnot( this->value.data.str, that->value.data.str );
 2201      break;
 2202       }
 2203       this->operation = CONST_OP;
 2204 
 2205    } else {
 2206 
 2207       Allocate_Ptrs( this );
 2208 
 2209       if( !gParse.status ) {
 2210 
 2211      if( this->type!=BITSTR ) {
 2212         elem = gParse.nRows;
 2213         if( this->type!=STRING )
 2214            elem *= this->value.nelem;
 2215         while( elem-- )
 2216            this->value.undef[elem] = that->value.undef[elem];
 2217      }
 2218 
 2219      elem = gParse.nRows * this->value.nelem;
 2220 
 2221      switch( this->operation ) {
 2222 
 2223      case BOOLEAN:
 2224         if( that->type==DOUBLE )
 2225            while( elem-- )
 2226           this->value.data.logptr[elem] =
 2227              ( that->value.data.dblptr[elem] != 0.0 );
 2228         else if( that->type==LONG )
 2229            while( elem-- )
 2230           this->value.data.logptr[elem] =
 2231              ( that->value.data.lngptr[elem] != 0L );
 2232         break;
 2233 
 2234      case DOUBLE:
 2235      case FLTCAST:
 2236         if( that->type==LONG )
 2237            while( elem-- )
 2238           this->value.data.dblptr[elem] =
 2239              (double)that->value.data.lngptr[elem];
 2240         else if( that->type==BOOLEAN )
 2241            while( elem-- )
 2242           this->value.data.dblptr[elem] =
 2243              ( that->value.data.logptr[elem] ? 1.0 : 0.0 );
 2244         break;
 2245 
 2246      case LONG:
 2247      case INTCAST:
 2248         if( that->type==DOUBLE )
 2249            while( elem-- )
 2250           this->value.data.lngptr[elem] =
 2251              (long)that->value.data.dblptr[elem];
 2252         else if( that->type==BOOLEAN )
 2253            while( elem-- )
 2254           this->value.data.lngptr[elem] =
 2255              ( that->value.data.logptr[elem] ? 1L : 0L );
 2256         break;
 2257 
 2258      case UMINUS:
 2259         if( that->type==DOUBLE ) {
 2260            while( elem-- )
 2261           this->value.data.dblptr[elem] =
 2262              - that->value.data.dblptr[elem];
 2263         } else if( that->type==LONG ) {
 2264            while( elem-- )
 2265           this->value.data.lngptr[elem] =
 2266              - that->value.data.lngptr[elem];
 2267         }
 2268         break;
 2269 
 2270      case NOT:
 2271         if( that->type==BOOLEAN ) {
 2272            while( elem-- )
 2273           this->value.data.logptr[elem] =
 2274              ( ! that->value.data.logptr[elem] );
 2275         } else if( that->type==BITSTR ) {
 2276            elem = gParse.nRows;
 2277            while( elem-- )
 2278           bitnot( this->value.data.strptr[elem],
 2279               that->value.data.strptr[elem] );
 2280         }
 2281         break;
 2282      }
 2283       }
 2284    }
 2285 
 2286    if( that->operation>0 ) {
 2287       free( that->value.data.ptr );
 2288    }
 2289 }
 2290 
 2291 static void Do_Offset( Node *this )
 2292 {
 2293    Node *col;
 2294    long fRow, nRowOverlap, nRowReload, rowOffset;
 2295    long nelem, elem, offset, nRealElem;
 2296    int status;
 2297 
 2298    col       = gParse.Nodes + this->SubNodes[0];
 2299    rowOffset = gParse.Nodes[  this->SubNodes[1] ].value.data.lng;
 2300 
 2301    Allocate_Ptrs( this );
 2302 
 2303    fRow   = gParse.firstRow + rowOffset;
 2304    if( this->type==STRING || this->type==BITSTR )
 2305       nRealElem = 1;
 2306    else
 2307       nRealElem = this->value.nelem;
 2308 
 2309    nelem = nRealElem;
 2310 
 2311    if( fRow < gParse.firstDataRow ) {
 2312 
 2313       /* Must fill in data at start of array */
 2314 
 2315       nRowReload = gParse.firstDataRow - fRow;
 2316       if( nRowReload > gParse.nRows ) nRowReload = gParse.nRows;
 2317       nRowOverlap = gParse.nRows - nRowReload;
 2318 
 2319       offset = 0;
 2320 
 2321       /*  NULLify any values falling out of bounds  */
 2322 
 2323       while( fRow<1 && nRowReload>0 ) {
 2324      if( this->type == BITSTR ) {
 2325         nelem = this->value.nelem;
 2326         this->value.data.strptr[offset][ nelem ] = '\0';
 2327         while( nelem-- ) this->value.data.strptr[offset][nelem] = '0';
 2328         offset++;
 2329      } else {
 2330         while( nelem-- )
 2331            this->value.undef[offset++] = 1;
 2332      }
 2333      nelem = nRealElem;
 2334      fRow++;
 2335      nRowReload--;
 2336       }
 2337 
 2338    } else if( fRow + gParse.nRows > gParse.firstDataRow + gParse.nDataRows ) {
 2339 
 2340       /* Must fill in data at end of array */
 2341 
 2342       nRowReload = (fRow+gParse.nRows) - (gParse.firstDataRow+gParse.nDataRows);
 2343       if( nRowReload>gParse.nRows ) {
 2344      nRowReload = gParse.nRows;
 2345       } else {
 2346      fRow = gParse.firstDataRow + gParse.nDataRows;
 2347       }
 2348       nRowOverlap = gParse.nRows - nRowReload;
 2349 
 2350       offset = nRowOverlap * nelem;
 2351 
 2352       /*  NULLify any values falling out of bounds  */
 2353 
 2354       elem = gParse.nRows * nelem;
 2355       while( fRow+nRowReload>gParse.totalRows && nRowReload>0 ) {
 2356      if( this->type == BITSTR ) {
 2357         nelem = this->value.nelem;
 2358         elem--;
 2359         this->value.data.strptr[elem][ nelem ] = '\0';
 2360         while( nelem-- ) this->value.data.strptr[elem][nelem] = '0';
 2361      } else {
 2362         while( nelem-- )
 2363            this->value.undef[--elem] = 1;
 2364      }
 2365      nelem = nRealElem;
 2366      nRowReload--;
 2367       }
 2368 
 2369    } else {
 2370 
 2371       nRowReload  = 0;
 2372       nRowOverlap = gParse.nRows;
 2373       offset      = 0;
 2374 
 2375    }
 2376 
 2377    if( nRowReload>0 ) {
 2378       switch( this->type ) {
 2379       case BITSTR:
 2380       case STRING:
 2381      status = (*gParse.loadData)( -col->operation, fRow, nRowReload,
 2382                       this->value.data.strptr+offset,
 2383                       this->value.undef+offset );
 2384      break;
 2385       case BOOLEAN:
 2386      status = (*gParse.loadData)( -col->operation, fRow, nRowReload,
 2387                       this->value.data.logptr+offset,
 2388                       this->value.undef+offset );
 2389      break;
 2390       case LONG:
 2391      status = (*gParse.loadData)( -col->operation, fRow, nRowReload,
 2392                       this->value.data.lngptr+offset,
 2393                       this->value.undef+offset );
 2394      break;
 2395       case DOUBLE:
 2396      status = (*gParse.loadData)( -col->operation, fRow, nRowReload,
 2397                       this->value.data.dblptr+offset,
 2398                       this->value.undef+offset );
 2399      break;
 2400       }
 2401    }
 2402 
 2403    /*  Now copy over the overlapping region, if any  */
 2404 
 2405    if( nRowOverlap <= 0 ) return;
 2406 
 2407    if( rowOffset>0 )
 2408       elem = nRowOverlap * nelem;
 2409    else
 2410       elem = gParse.nRows * nelem;
 2411 
 2412    offset = nelem * rowOffset;
 2413    while( nRowOverlap-- && !gParse.status ) {
 2414       while( nelem-- && !gParse.status ) {
 2415      elem--;
 2416      if( this->type != BITSTR )
 2417         this->value.undef[elem] = col->value.undef[elem+offset];
 2418      switch( this->type ) {
 2419      case BITSTR:
 2420         strcpy( this->value.data.strptr[elem       ],
 2421                      col->value.data.strptr[elem+offset] );
 2422         break;
 2423      case STRING:
 2424         strcpy( this->value.data.strptr[elem       ],
 2425                      col->value.data.strptr[elem+offset] );
 2426         break;
 2427      case BOOLEAN:
 2428         this->value.data.logptr[elem] = col->value.data.logptr[elem+offset];
 2429         break;
 2430      case LONG:
 2431         this->value.data.lngptr[elem] = col->value.data.lngptr[elem+offset];
 2432         break;
 2433      case DOUBLE:
 2434         this->value.data.dblptr[elem] = col->value.data.dblptr[elem+offset];
 2435         break;
 2436      }
 2437       }
 2438       nelem = nRealElem;
 2439    }
 2440 }
 2441 
 2442 static void Do_BinOp_bit( Node *this )
 2443 {
 2444    Node *that1, *that2;
 2445    char *sptr1=NULL, *sptr2=NULL;
 2446    int  const1, const2;
 2447    long rows;
 2448 
 2449    that1 = gParse.Nodes + this->SubNodes[0];
 2450    that2 = gParse.Nodes + this->SubNodes[1];
 2451 
 2452    const1 = ( that1->operation==CONST_OP );
 2453    const2 = ( that2->operation==CONST_OP );
 2454    sptr1  = ( const1 ? that1->value.data.str : NULL );
 2455    sptr2  = ( const2 ? that2->value.data.str : NULL );
 2456 
 2457    if( const1 && const2 ) {
 2458       switch( this->operation ) {
 2459       case NE:
 2460      this->value.data.log = !bitcmp( sptr1, sptr2 );
 2461      break;
 2462       case EQ:
 2463      this->value.data.log =  bitcmp( sptr1, sptr2 );
 2464      break;
 2465       case GT:
 2466       case LT:
 2467       case LTE:
 2468       case GTE:
 2469      this->value.data.log = bitlgte( sptr1, this->operation, sptr2 );
 2470      break;
 2471       case '|': 
 2472      bitor( this->value.data.str, sptr1, sptr2 );
 2473      break;
 2474       case '&': 
 2475      bitand( this->value.data.str, sptr1, sptr2 );
 2476      break;
 2477       case '+':
 2478      strcpy( this->value.data.str, sptr1 );
 2479      strcat( this->value.data.str, sptr2 );
 2480      break;
 2481       case ACCUM:
 2482     this->value.data.lng = 0;
 2483     while( *sptr1 ) {
 2484       if ( *sptr1 == '1' ) this->value.data.lng ++;
 2485       sptr1 ++;
 2486     }
 2487     break;
 2488     
 2489       }
 2490       this->operation = CONST_OP;
 2491 
 2492    } else {
 2493 
 2494       Allocate_Ptrs( this );
 2495 
 2496       if( !gParse.status ) {
 2497      rows  = gParse.nRows;
 2498      switch( this->operation ) {
 2499 
 2500         /*  BITSTR comparisons  */
 2501 
 2502      case NE:
 2503      case EQ:
 2504      case GT:
 2505      case LT:
 2506      case LTE:
 2507      case GTE:
 2508         while( rows-- ) {
 2509            if( !const1 )
 2510           sptr1 = that1->value.data.strptr[rows];
 2511            if( !const2 )
 2512           sptr2 = that2->value.data.strptr[rows];
 2513            switch( this->operation ) {
 2514            case NE:  this->value.data.logptr[rows] = 
 2515                                                       !bitcmp( sptr1, sptr2 );
 2516                          break;
 2517            case EQ:  this->value.data.logptr[rows] = 
 2518                                                        bitcmp( sptr1, sptr2 );
 2519                          break;
 2520            case GT:
 2521            case LT:
 2522            case LTE:
 2523            case GTE: this->value.data.logptr[rows] = 
 2524                                      bitlgte( sptr1, this->operation, sptr2 );
 2525                      break;
 2526            }
 2527            this->value.undef[rows] = 0;
 2528         }
 2529         break;
 2530      
 2531         /*  BITSTR AND/ORs ...  no UNDEFS in or out */
 2532       
 2533      case '|': 
 2534      case '&': 
 2535      case '+':
 2536         while( rows-- ) {
 2537            if( !const1 )
 2538           sptr1 = that1->value.data.strptr[rows];
 2539            if( !const2 )
 2540           sptr2 = that2->value.data.strptr[rows];
 2541            if( this->operation=='|' )
 2542           bitor(  this->value.data.strptr[rows], sptr1, sptr2 );
 2543            else if( this->operation=='&' )
 2544           bitand( this->value.data.strptr[rows], sptr1, sptr2 );
 2545            else {
 2546           strcpy( this->value.data.strptr[rows], sptr1 );
 2547           strcat( this->value.data.strptr[rows], sptr2 );
 2548            }
 2549         }
 2550         break;
 2551 
 2552         /* Accumulate 1 bits */
 2553      case ACCUM:
 2554        { 
 2555          long i, previous, curr;
 2556 
 2557          previous = that2->value.data.lng;
 2558          
 2559           /* Cumulative sum of this chunk */
 2560          for (i=0; i<rows; i++) {
 2561            sptr1 = that1->value.data.strptr[i];
 2562            for (curr = 0; *sptr1; sptr1 ++) {
 2563          if ( *sptr1 == '1' ) curr ++;
 2564            }
 2565            previous += curr;
 2566            this->value.data.lngptr[i] = previous;
 2567            this->value.undef[i] = 0;
 2568          }
 2569          
 2570           /* Store final cumulant for next pass */
 2571          that2->value.data.lng = previous;
 2572        }
 2573      }
 2574       }
 2575    }
 2576 
 2577    if( that1->operation>0 ) {
 2578       free( that1->value.data.strptr[0] );
 2579       free( that1->value.data.strptr    );
 2580    }
 2581    if( that2->operation>0 ) {
 2582       free( that2->value.data.strptr[0] );
 2583       free( that2->value.data.strptr    );
 2584    }
 2585 }
 2586 
 2587 static void Do_BinOp_str( Node *this )
 2588 {
 2589    Node *that1, *that2;
 2590    char *sptr1, *sptr2, null1=0, null2=0;
 2591    int const1, const2, val;
 2592    long rows;
 2593 
 2594    that1 = gParse.Nodes + this->SubNodes[0];
 2595    that2 = gParse.Nodes + this->SubNodes[1];
 2596 
 2597    const1 = ( that1->operation==CONST_OP );
 2598    const2 = ( that2->operation==CONST_OP );
 2599    sptr1  = ( const1 ? that1->value.data.str : NULL );
 2600    sptr2  = ( const2 ? that2->value.data.str : NULL );
 2601 
 2602    if( const1 && const2 ) {  /*  Result is a constant  */
 2603       switch( this->operation ) {
 2604 
 2605      /*  Compare Strings  */
 2606 
 2607       case NE:
 2608       case EQ:
 2609      val = ( FSTRCMP( sptr1, sptr2 ) == 0 );
 2610      this->value.data.log = ( this->operation==EQ ? val : !val );
 2611      break;
 2612       case GT:
 2613      this->value.data.log = ( FSTRCMP( sptr1, sptr2 ) > 0 );
 2614      break;
 2615       case LT:
 2616      this->value.data.log = ( FSTRCMP( sptr1, sptr2 ) < 0 );
 2617      break;
 2618       case GTE:
 2619      this->value.data.log = ( FSTRCMP( sptr1, sptr2 ) >= 0 );
 2620      break;
 2621       case LTE:
 2622      this->value.data.log = ( FSTRCMP( sptr1, sptr2 ) <= 0 );
 2623      break;
 2624 
 2625      /*  Concat Strings  */
 2626 
 2627       case '+':
 2628      strcpy( this->value.data.str, sptr1 );
 2629      strcat( this->value.data.str, sptr2 );
 2630      break;
 2631       }
 2632       this->operation = CONST_OP;
 2633 
 2634    } else {  /*  Not a constant  */
 2635 
 2636       Allocate_Ptrs( this );
 2637 
 2638       if( !gParse.status ) {
 2639 
 2640      rows = gParse.nRows;
 2641      switch( this->operation ) {
 2642 
 2643         /*  Compare Strings  */
 2644 
 2645      case NE:
 2646      case EQ:
 2647         while( rows-- ) {
 2648            if( !const1 ) null1 = that1->value.undef[rows];
 2649            if( !const2 ) null2 = that2->value.undef[rows];
 2650            this->value.undef[rows] = (null1 || null2);
 2651            if( ! this->value.undef[rows] ) {
 2652           if( !const1 ) sptr1  = that1->value.data.strptr[rows];
 2653           if( !const2 ) sptr2  = that2->value.data.strptr[rows];
 2654           val = ( FSTRCMP( sptr1, sptr2 ) == 0 );
 2655           this->value.data.logptr[rows] =
 2656              ( this->operation==EQ ? val : !val );
 2657            }
 2658         }
 2659         break;
 2660         
 2661      case GT:
 2662      case LT:
 2663         while( rows-- ) {
 2664            if( !const1 ) null1 = that1->value.undef[rows];
 2665            if( !const2 ) null2 = that2->value.undef[rows];
 2666            this->value.undef[rows] = (null1 || null2);
 2667            if( ! this->value.undef[rows] ) {
 2668           if( !const1 ) sptr1  = that1->value.data.strptr[rows];
 2669           if( !const2 ) sptr2  = that2->value.data.strptr[rows];
 2670           val = ( FSTRCMP( sptr1, sptr2 ) );
 2671           this->value.data.logptr[rows] =
 2672              ( this->operation==GT ? val>0 : val<0 );
 2673            }
 2674         }
 2675         break;
 2676 
 2677      case GTE:
 2678      case LTE:
 2679         while( rows-- ) {
 2680            if( !const1 ) null1 = that1->value.undef[rows];
 2681            if( !const2 ) null2 = that2->value.undef[rows];
 2682            this->value.undef[rows] = (null1 || null2);
 2683            if( ! this->value.undef[rows] ) {
 2684           if( !const1 ) sptr1  = that1->value.data.strptr[rows];
 2685           if( !const2 ) sptr2  = that2->value.data.strptr[rows];
 2686           val = ( FSTRCMP( sptr1, sptr2 ) );
 2687           this->value.data.logptr[rows] =
 2688              ( this->operation==GTE ? val>=0 : val<=0 );
 2689            }
 2690         }
 2691         break;
 2692 
 2693         /*  Concat Strings  */
 2694         
 2695      case '+':
 2696         while( rows-- ) {
 2697            if( !const1 ) null1 = that1->value.undef[rows];
 2698            if( !const2 ) null2 = that2->value.undef[rows];
 2699            this->value.undef[rows] = (null1 || null2);
 2700            if( ! this->value.undef[rows] ) {
 2701           if( !const1 ) sptr1  = that1->value.data.strptr[rows];
 2702           if( !const2 ) sptr2  = that2->value.data.strptr[rows];
 2703           strcpy( this->value.data.strptr[rows], sptr1 );
 2704           strcat( this->value.data.strptr[rows], sptr2 );
 2705            }
 2706         }
 2707         break;
 2708      }
 2709       }
 2710    }
 2711 
 2712    if( that1->operation>0 ) {
 2713       free( that1->value.data.strptr[0] );
 2714       free( that1->value.data.strptr );
 2715    }
 2716    if( that2->operation>0 ) {
 2717       free( that2->value.data.strptr[0] );
 2718       free( that2->value.data.strptr );
 2719    }
 2720 }
 2721 
 2722 static void Do_BinOp_log( Node *this )
 2723 {
 2724    Node *that1, *that2;
 2725    int vector1, vector2;
 2726    char val1=0, val2=0, null1=0, null2=0;
 2727    long rows, nelem, elem;
 2728 
 2729    that1 = gParse.Nodes + this->SubNodes[0];
 2730    that2 = gParse.Nodes + this->SubNodes[1];
 2731 
 2732    vector1 = ( that1->operation!=CONST_OP );
 2733    if( vector1 )
 2734       vector1 = that1->value.nelem;
 2735    else {
 2736       val1  = that1->value.data.log;
 2737    }
 2738 
 2739    vector2 = ( that2->operation!=CONST_OP );
 2740    if( vector2 )
 2741       vector2 = that2->value.nelem;
 2742    else {
 2743       val2  = that2->value.data.log;
 2744    }
 2745 
 2746    if( !vector1 && !vector2 ) {  /*  Result is a constant  */
 2747       switch( this->operation ) {
 2748       case OR:
 2749      this->value.data.log = (val1 || val2);
 2750      break;
 2751       case AND:
 2752      this->value.data.log = (val1 && val2);
 2753      break;
 2754       case EQ:
 2755      this->value.data.log = ( (val1 && val2) || (!val1 && !val2) );
 2756      break;
 2757       case NE:
 2758      this->value.data.log = ( (val1 && !val2) || (!val1 && val2) );
 2759      break;
 2760       case ACCUM:
 2761      this->value.data.lng = val1;
 2762      break;
 2763       }
 2764       this->operation=CONST_OP;
 2765    } else if (this->operation == ACCUM) {
 2766       long i, previous, curr;
 2767       rows  = gParse.nRows;
 2768       nelem = this->value.nelem;
 2769       elem  = this->value.nelem * rows;
 2770       
 2771       Allocate_Ptrs( this );
 2772       
 2773       if( !gParse.status ) {
 2774     previous = that2->value.data.lng;
 2775     
 2776     /* Cumulative sum of this chunk */
 2777     for (i=0; i<elem; i++) {
 2778       if (!that1->value.undef[i]) {
 2779         curr = that1->value.data.logptr[i];
 2780         previous += curr;
 2781       }
 2782       this->value.data.lngptr[i] = previous;
 2783       this->value.undef[i] = 0;
 2784     }
 2785     
 2786     /* Store final cumulant for next pass */
 2787     that2->value.data.lng = previous;
 2788       }
 2789       
 2790    } else {
 2791       rows  = gParse.nRows;
 2792       nelem = this->value.nelem;
 2793       elem  = this->value.nelem * rows;
 2794 
 2795       Allocate_Ptrs( this );
 2796 
 2797       if( !gParse.status ) {
 2798     
 2799      if (this->operation == ACCUM) {
 2800        long i, previous, curr;
 2801        
 2802        previous = that2->value.data.lng;
 2803        
 2804        /* Cumulative sum of this chunk */
 2805        for (i=0; i<elem; i++) {
 2806          if (!that1->value.undef[i]) {
 2807            curr = that1->value.data.logptr[i];
 2808            previous += curr;
 2809          }
 2810          this->value.data.lngptr[i] = previous;
 2811          this->value.undef[i] = 0;
 2812        }
 2813        
 2814        /* Store final cumulant for next pass */
 2815        that2->value.data.lng = previous;
 2816      }
 2817     
 2818      while( rows-- ) {
 2819         while( nelem-- ) {
 2820            elem--;
 2821 
 2822            if( vector1>1 ) {
 2823           val1  = that1->value.data.logptr[elem];
 2824           null1 = that1->value.undef[elem];
 2825            } else if( vector1 ) {
 2826           val1  = that1->value.data.logptr[rows];
 2827           null1 = that1->value.undef[rows];
 2828            }
 2829 
 2830            if( vector2>1 ) {
 2831           val2  = that2->value.data.logptr[elem];
 2832           null2 = that2->value.undef[elem];
 2833            } else if( vector2 ) {
 2834           val2  = that2->value.data.logptr[rows];
 2835           null2 = that2->value.undef[rows];
 2836            }
 2837 
 2838            this->value.undef[elem] = (null1 || null2);
 2839            switch( this->operation ) {
 2840 
 2841            case OR:
 2842           /*  This is more complicated than others to suppress UNDEFs */
 2843           /*  in those cases where the other argument is DEF && TRUE  */
 2844 
 2845           if( !null1 && !null2 ) {
 2846              this->value.data.logptr[elem] = (val1 || val2);
 2847           } else if( (null1 && !null2 && val2)
 2848                  || ( !null1 && null2 && val1 ) ) {
 2849              this->value.data.logptr[elem] = 1;
 2850              this->value.undef[elem] = 0;
 2851           }
 2852           break;
 2853 
 2854            case AND:
 2855           /*  This is more complicated than others to suppress UNDEFs */
 2856           /*  in those cases where the other argument is DEF && FALSE */
 2857 
 2858           if( !null1 && !null2 ) {
 2859              this->value.data.logptr[elem] = (val1 && val2);
 2860           } else if( (null1 && !null2 && !val2)
 2861                  || ( !null1 && null2 && !val1 ) ) {
 2862              this->value.data.logptr[elem] = 0;
 2863              this->value.undef[elem] = 0;
 2864           }
 2865           break;
 2866 
 2867            case EQ:
 2868           this->value.data.logptr[elem] = 
 2869              ( (val1 && val2) || (!val1 && !val2) );
 2870           break;
 2871 
 2872            case NE:
 2873           this->value.data.logptr[elem] =
 2874              ( (val1 && !val2) || (!val1 && val2) );
 2875           break;
 2876            }
 2877         }
 2878         nelem = this->value.nelem;
 2879      }
 2880       }
 2881    }
 2882 
 2883    if( that1->operation>0 ) {
 2884       free( that1->value.data.ptr );
 2885    }
 2886    if( that2->operation>0 ) {
 2887       free( that2->value.data.ptr );
 2888    }
 2889 }
 2890 
 2891 static void Do_BinOp_lng( Node *this )
 2892 {
 2893    Node *that1, *that2;
 2894    int  vector1, vector2;
 2895    long val1=0, val2=0;
 2896    char null1=0, null2=0;
 2897    long rows, nelem, elem;
 2898 
 2899    that1 = gParse.Nodes + this->SubNodes[0];
 2900    that2 = gParse.Nodes + this->SubNodes[1];
 2901 
 2902    vector1 = ( that1->operation!=CONST_OP );
 2903    if( vector1 )
 2904       vector1 = that1->value.nelem;
 2905    else {
 2906       val1  = that1->value.data.lng;
 2907    }
 2908 
 2909    vector2 = ( that2->operation!=CONST_OP );
 2910    if( vector2 )
 2911       vector2 = that2->value.nelem;
 2912    else {
 2913       val2  = that2->value.data.lng;
 2914    }
 2915 
 2916    if( !vector1 && !vector2 ) {  /*  Result is a constant  */
 2917 
 2918       switch( this->operation ) {
 2919       case '~':   /* Treat as == for LONGS */
 2920       case EQ:    this->value.data.log = (val1 == val2);   break;
 2921       case NE:    this->value.data.log = (val1 != val2);   break;
 2922       case GT:    this->value.data.log = (val1 >  val2);   break;
 2923       case LT:    this->value.data.log = (val1 <  val2);   break;
 2924       case LTE:   this->value.data.log = (val1 <= val2);   break;
 2925       case GTE:   this->value.data.log = (val1 >= val2);   break;
 2926 
 2927       case '+':   this->value.data.lng = (val1  + val2);   break;
 2928       case '-':   this->value.data.lng = (val1  - val2);   break;
 2929       case '*':   this->value.data.lng = (val1  * val2);   break;
 2930 
 2931       case '&':   this->value.data.lng = (val1  & val2);   break;
 2932       case '|':   this->value.data.lng = (val1  | val2);   break;
 2933       case '^':   this->value.data.lng = (val1  ^ val2);   break;
 2934 
 2935       case '%':
 2936      if( val2 ) this->value.data.lng = (val1 % val2);
 2937      else       yyerror("Divide by Zero");
 2938      break;
 2939       case '/': 
 2940      if( val2 ) this->value.data.lng = (val1 / val2); 
 2941      else       yyerror("Divide by Zero");
 2942      break;
 2943       case POWER:
 2944      this->value.data.lng = (long)pow((double)val1,(double)val2);
 2945      break;
 2946       case ACCUM:
 2947      this->value.data.lng = val1;
 2948      break;
 2949       case DIFF:
 2950      this->value.data.lng = 0;
 2951      break;
 2952       }
 2953       this->operation=CONST_OP;
 2954 
 2955    } else if ((this->operation == ACCUM) || (this->operation == DIFF)) {
 2956       long i, previous, curr;
 2957       long undef;
 2958       rows  = gParse.nRows;
 2959       nelem = this->value.nelem;
 2960       elem  = this->value.nelem * rows;
 2961       
 2962       Allocate_Ptrs( this );
 2963       
 2964       if( !gParse.status ) {
 2965     previous = that2->value.data.lng;
 2966     undef    = (long) that2->value.undef;
 2967     
 2968     if (this->operation == ACCUM) {
 2969       /* Cumulative sum of this chunk */
 2970       for (i=0; i<elem; i++) {
 2971         if (!that1->value.undef[i]) {
 2972           curr = that1->value.data.lngptr[i];
 2973           previous += curr;
 2974         }
 2975         this->value.data.lngptr[i] = previous;
 2976         this->value.undef[i] = 0;
 2977       }
 2978     } else {
 2979       /* Sequential difference for this chunk */
 2980       for (i=0; i<elem; i++) {
 2981         curr = that1->value.data.lngptr[i];
 2982         if (that1->value.undef[i] || undef) {
 2983           /* Either this, or previous, value was undefined */
 2984           this->value.data.lngptr[i] = 0;
 2985           this->value.undef[i] = 1;
 2986         } else {
 2987           /* Both defined, we are okay! */
 2988           this->value.data.lngptr[i] = curr - previous;
 2989           this->value.undef[i] = 0;
 2990         }
 2991 
 2992         previous = curr;
 2993         undef = that1->value.undef[i];
 2994       }
 2995     }     
 2996     
 2997     /* Store final cumulant for next pass */
 2998     that2->value.data.lng = previous;
 2999     that2->value.undef    = (char *) undef; /* XXX evil, but no harm here */
 3000       }
 3001       
 3002    } else {
 3003 
 3004       rows  = gParse.nRows;
 3005       nelem = this->value.nelem;
 3006       elem  = this->value.nelem * rows;
 3007 
 3008       Allocate_Ptrs( this );
 3009 
 3010       while( rows-- && !gParse.status ) {
 3011      while( nelem-- && !gParse.status ) {
 3012         elem--;
 3013 
 3014         if( vector1>1 ) {
 3015            val1  = that1->value.data.lngptr[elem];
 3016            null1 = that1->value.undef[elem];
 3017         } else if( vector1 ) {
 3018            val1  = that1->value.data.lngptr[rows];
 3019            null1 = that1->value.undef[rows];
 3020         }
 3021 
 3022         if( vector2>1 ) {
 3023            val2  = that2->value.data.lngptr[elem];
 3024            null2 = that2->value.undef[elem];
 3025         } else if( vector2 ) {
 3026            val2  = that2->value.data.lngptr[rows];
 3027            null2 = that2->value.undef[rows];
 3028         }
 3029 
 3030         this->value.undef[elem] = (null1 || null2);
 3031         switch( this->operation ) {
 3032         case '~':   /* Treat as == for LONGS */
 3033         case EQ:   this->value.data.logptr[elem] = (val1 == val2);   break;
 3034         case NE:   this->value.data.logptr[elem] = (val1 != val2);   break;
 3035         case GT:   this->value.data.logptr[elem] = (val1 >  val2);   break;
 3036         case LT:   this->value.data.logptr[elem] = (val1 <  val2);   break;
 3037         case LTE:  this->value.data.logptr[elem] = (val1 <= val2);   break;
 3038         case GTE:  this->value.data.logptr[elem] = (val1 >= val2);   break;
 3039            
 3040         case '+':  this->value.data.lngptr[elem] = (val1  + val2);   break;
 3041         case '-':  this->value.data.lngptr[elem] = (val1  - val2);   break;
 3042         case '*':  this->value.data.lngptr[elem] = (val1  * val2);   break;
 3043 
 3044         case '&':  this->value.data.lngptr[elem] = (val1  & val2);   break;
 3045         case '|':  this->value.data.lngptr[elem] = (val1  | val2);   break;
 3046         case '^':  this->value.data.lngptr[elem] = (val1  ^ val2);   break;
 3047 
 3048         case '%':   
 3049            if( val2 ) this->value.data.lngptr[elem] = (val1 % val2);
 3050            else {
 3051          this->value.data.lngptr[elem] = 0;
 3052          this->value.undef[elem] = 1;
 3053            }
 3054            break;
 3055         case '/': 
 3056            if( val2 ) this->value.data.lngptr[elem] = (val1 / val2); 
 3057            else {
 3058          this->value.data.lngptr[elem] = 0;
 3059          this->value.undef[elem] = 1;
 3060            }
 3061            break;
 3062         case POWER:
 3063            this->value.data.lngptr[elem] = (long)pow((double)val1,(double)val2);
 3064            break;
 3065         }
 3066      }
 3067      nelem = this->value.nelem;
 3068       }
 3069    }
 3070 
 3071    if( that1->operation>0 ) {
 3072       free( that1->value.data.ptr );
 3073    }
 3074    if( that2->operation>0 ) {
 3075       free( that2->value.data.ptr );
 3076    }
 3077 }
 3078 
 3079 static void Do_BinOp_dbl( Node *this )
 3080 {
 3081    Node   *that1, *that2;
 3082    int    vector1, vector2;
 3083    double val1=0.0, val2=0.0;
 3084    char   null1=0, null2=0;
 3085    long   rows, nelem, elem;
 3086 
 3087    that1 = gParse.Nodes + this->SubNodes[0];
 3088    that2 = gParse.Nodes + this->SubNodes[1];
 3089 
 3090    vector1 = ( that1->operation!=CONST_OP );
 3091    if( vector1 )
 3092       vector1 = that1->value.nelem;
 3093    else {
 3094       val1  = that1->value.data.dbl;
 3095    }
 3096 
 3097    vector2 = ( that2->operation!=CONST_OP );
 3098    if( vector2 )
 3099       vector2 = that2->value.nelem;
 3100    else {
 3101       val2  = that2->value.data.dbl;
 3102    } 
 3103 
 3104    if( !vector1 && !vector2 ) {  /*  Result is a constant  */
 3105 
 3106       switch( this->operation ) {
 3107       case '~':   this->value.data.log = ( fabs(val1-val2) < APPROX );   break;
 3108       case EQ:    this->value.data.log = (val1 == val2);   break;
 3109       case NE:    this->value.data.log = (val1 != val2);   break;
 3110       case GT:    this->value.data.log = (val1 >  val2);   break;
 3111       case LT:    this->value.data.log = (val1 <  val2);   break;
 3112       case LTE:   this->value.data.log = (val1 <= val2);   break;
 3113       case GTE:   this->value.data.log = (val1 >= val2);   break;
 3114 
 3115       case '+':   this->value.data.dbl = (val1  + val2);   break;
 3116       case '-':   this->value.data.dbl = (val1  - val2);   break;
 3117       case '*':   this->value.data.dbl = (val1  * val2);   break;
 3118 
 3119       case '%':
 3120      if( val2 ) this->value.data.dbl = val1 - val2*((int)(val1/val2));
 3121      else       yyerror("Divide by Zero");
 3122      break;
 3123       case '/': 
 3124      if( val2 ) this->value.data.dbl = (val1 / val2); 
 3125      else       yyerror("Divide by Zero");
 3126      break;
 3127       case POWER:
 3128      this->value.data.dbl = (double)pow(val1,val2);
 3129      break;
 3130       case ACCUM:
 3131      this->value.data.dbl = val1;
 3132      break;
 3133       case DIFF:
 3134     this->value.data.dbl = 0;
 3135      break;
 3136       }
 3137       this->operation=CONST_OP;
 3138 
 3139    } else if ((this->operation == ACCUM) || (this->operation == DIFF)) {
 3140       long i;
 3141       long undef;
 3142       double previous, curr;
 3143       rows  = gParse.nRows;
 3144       nelem = this->value.nelem;
 3145       elem  = this->value.nelem * rows;
 3146       
 3147       Allocate_Ptrs( this );
 3148       
 3149       if( !gParse.status ) {
 3150     previous = that2->value.data.dbl;
 3151     undef    = (long) that2->value.undef;
 3152     
 3153     if (this->operation == ACCUM) {
 3154       /* Cumulative sum of this chunk */
 3155       for (i=0; i<elem; i++) {
 3156         if (!that1->value.undef[i]) {
 3157           curr = that1->value.data.dblptr[i];
 3158           previous += curr;
 3159         }
 3160         this->value.data.dblptr[i] = previous;
 3161         this->value.undef[i] = 0;
 3162       }
 3163     } else {
 3164       /* Sequential difference for this chunk */
 3165       for (i=0; i<elem; i++) {
 3166         curr = that1->value.data.dblptr[i];
 3167         if (that1->value.undef[i] || undef) {
 3168           /* Either this, or previous, value was undefined */
 3169           this->value.data.dblptr[i] = 0;
 3170           this->value.undef[i] = 1;
 3171         } else {
 3172           /* Both defined, we are okay! */
 3173           this->value.data.dblptr[i] = curr - previous;
 3174           this->value.undef[i] = 0;
 3175         }
 3176 
 3177         previous = curr;
 3178         undef = that1->value.undef[i];
 3179       }
 3180     }     
 3181     
 3182     /* Store final cumulant for next pass */
 3183     that2->value.data.dbl = previous;
 3184     that2->value.undef    = (char *) undef; /* XXX evil, but no harm here */
 3185       }
 3186       
 3187    } else {
 3188 
 3189       rows  = gParse.nRows;
 3190       nelem = this->value.nelem;
 3191       elem  = this->value.nelem * rows;
 3192 
 3193       Allocate_Ptrs( this );
 3194 
 3195       while( rows-- && !gParse.status ) {
 3196      while( nelem-- && !gParse.status ) {
 3197         elem--;
 3198 
 3199         if( vector1>1 ) {
 3200            val1  = that1->value.data.dblptr[elem];
 3201            null1 = that1->value.undef[elem];
 3202         } else if( vector1 ) {
 3203            val1  = that1->value.data.dblptr[rows];
 3204            null1 = that1->value.undef[rows];
 3205         }
 3206 
 3207         if( vector2>1 ) {
 3208            val2  = that2->value.data.dblptr[elem];
 3209            null2 = that2->value.undef[elem];
 3210         } else if( vector2 ) {
 3211            val2  = that2->value.data.dblptr[rows];
 3212            null2 = that2->value.undef[rows];
 3213         }
 3214 
 3215         this->value.undef[elem] = (null1 || null2);
 3216         switch( this->operation ) {
 3217         case '~':   this->value.data.logptr[elem] =
 3218                                           ( fabs(val1-val2) < APPROX );   break;
 3219         case EQ:    this->value.data.logptr[elem] = (val1 == val2);   break;
 3220         case NE:    this->value.data.logptr[elem] = (val1 != val2);   break;
 3221         case GT:    this->value.data.logptr[elem] = (val1 >  val2);   break;
 3222         case LT:    this->value.data.logptr[elem] = (val1 <  val2);   break;
 3223         case LTE:   this->value.data.logptr[elem] = (val1 <= val2);   break;
 3224         case GTE:   this->value.data.logptr[elem] = (val1 >= val2);   break;
 3225            
 3226         case '+':   this->value.data.dblptr[elem] = (val1  + val2);   break;
 3227         case '-':   this->value.data.dblptr[elem] = (val1  - val2);   break;
 3228         case '*':   this->value.data.dblptr[elem] = (val1  * val2);   break;
 3229 
 3230         case '%':
 3231            if( val2 ) this->value.data.dblptr[elem] =
 3232                                 val1 - val2*((int)(val1/val2));
 3233            else {
 3234          this->value.data.dblptr[elem] = 0.0;
 3235          this->value.undef[elem] = 1;
 3236            }
 3237            break;
 3238         case '/': 
 3239            if( val2 ) this->value.data.dblptr[elem] = (val1 / val2); 
 3240            else {
 3241          this->value.data.dblptr[elem] = 0.0;
 3242          this->value.undef[elem] = 1;
 3243            }
 3244            break;
 3245         case POWER:
 3246            this->value.data.dblptr[elem] = (double)pow(val1,val2);
 3247            break;
 3248         }
 3249      }
 3250      nelem = this->value.nelem;
 3251       }
 3252    }
 3253 
 3254    if( that1->operation>0 ) {
 3255       free( that1->value.data.ptr );
 3256    }
 3257    if( that2->operation>0 ) {
 3258       free( that2->value.data.ptr );
 3259    }
 3260 }
 3261 
 3262 /*
 3263  *  This Quickselect routine is based on the algorithm described in
 3264  *  "Numerical recipes in C", Second Edition,
 3265  *  Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5
 3266  *  This code by Nicolas Devillard - 1998. Public domain.
 3267  * http://ndevilla.free.fr/median/median/src/quickselect.c
 3268  */
 3269 
 3270 #define ELEM_SWAP(a,b) { register long t=(a);(a)=(b);(b)=t; }
 3271 
 3272 /* 
 3273  * qselect_median_lng - select the median value of a long array
 3274  *
 3275  * This routine selects the median value of the long integer array
 3276  * arr[].  If there are an even number of elements, the "lower median"
 3277  * is selected.
 3278  *
 3279  * The array arr[] is scrambled, so users must operate on a scratch
 3280  * array if they wish the values to be preserved.
 3281  *
 3282  * long arr[] - array of values
 3283  * int n - number of elements in arr
 3284  *
 3285  * RETURNS: the lower median value of arr[]
 3286  *
 3287  */
 3288 long qselect_median_lng(long arr[], int n)
 3289 {
 3290     int low, high ;
 3291     int median;
 3292     int middle, ll, hh;
 3293 
 3294     low = 0 ; high = n-1 ; median = (low + high) / 2;
 3295     for (;;) {
 3296 
 3297         if (high <= low) { /* One element only */
 3298       return arr[median];     
 3299     }
 3300 
 3301         if (high == low + 1) {  /* Two elements only */
 3302             if (arr[low] > arr[high])
 3303                 ELEM_SWAP(arr[low], arr[high]) ;
 3304         return arr[median];
 3305         }
 3306 
 3307     /* Find median of low, middle and high items; swap into position low */
 3308     middle = (low + high) / 2;
 3309     if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]) ;
 3310     if (arr[low] > arr[high])       ELEM_SWAP(arr[low], arr[high]) ;
 3311     if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
 3312 
 3313     /* Swap low item (now in position middle) into position (low+1) */
 3314     ELEM_SWAP(arr[middle], arr[low+1]) ;
 3315 
 3316     /* Nibble from each end towards middle, swapping items when stuck */
 3317     ll = low + 1;
 3318     hh = high;
 3319     for (;;) {
 3320         do ll++; while (arr[low] > arr[ll]) ;
 3321         do hh--; while (arr[hh]  > arr[low]) ;
 3322 
 3323         if (hh < ll)
 3324         break;
 3325 
 3326         ELEM_SWAP(arr[ll], arr[hh]) ;
 3327     }
 3328 
 3329     /* Swap middle item (in position low) back into correct position */
 3330     ELEM_SWAP(arr[low], arr[hh]) ;
 3331 
 3332     /* Re-set active partition */
 3333     if (hh <= median)
 3334         low = ll;
 3335         if (hh >= median)
 3336         high = hh - 1;
 3337     }
 3338 }
 3339 
 3340 #undef ELEM_SWAP
 3341 
 3342 #define ELEM_SWAP(a,b) { register double t=(a);(a)=(b);(b)=t; }
 3343 
 3344 /* 
 3345  * qselect_median_dbl - select the median value of a double array
 3346  *
 3347  * This routine selects the median value of the double array
 3348  * arr[].  If there are an even number of elements, the "lower median"
 3349  * is selected.
 3350  *
 3351  * The array arr[] is scrambled, so users must operate on a scratch
 3352  * array if they wish the values to be preserved.
 3353  *
 3354  * double arr[] - array of values
 3355  * int n - number of elements in arr
 3356  *
 3357  * RETURNS: the lower median value of arr[]
 3358  *
 3359  */
 3360 double qselect_median_dbl(double arr[], int n)
 3361 {
 3362     int low, high ;
 3363     int median;
 3364     int middle, ll, hh;
 3365 
 3366     low = 0 ; high = n-1 ; median = (low + high) / 2;
 3367     for (;;) {
 3368         if (high <= low) { /* One element only */
 3369             return arr[median] ;
 3370     }
 3371 
 3372         if (high == low + 1) {  /* Two elements only */
 3373             if (arr[low] > arr[high])
 3374                 ELEM_SWAP(arr[low], arr[high]) ;
 3375             return arr[median] ;
 3376         }
 3377 
 3378     /* Find median of low, middle and high items; swap into position low */
 3379     middle = (low + high) / 2;
 3380     if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]) ;
 3381     if (arr[low] > arr[high])       ELEM_SWAP(arr[low], arr[high]) ;
 3382     if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
 3383 
 3384     /* Swap low item (now in position middle) into position (low+1) */
 3385     ELEM_SWAP(arr[middle], arr[low+1]) ;
 3386 
 3387     /* Nibble from each end towards middle, swapping items when stuck */
 3388     ll = low + 1;
 3389     hh = high;
 3390     for (;;) {
 3391         do ll++; while (arr[low] > arr[ll]) ;
 3392         do hh--; while (arr[hh]  > arr[low]) ;
 3393 
 3394         if (hh < ll)
 3395         break;
 3396 
 3397         ELEM_SWAP(arr[ll], arr[hh]) ;
 3398     }
 3399 
 3400     /* Swap middle item (in position low) back into correct position */
 3401     ELEM_SWAP(arr[low], arr[hh]) ;
 3402 
 3403     /* Re-set active partition */
 3404     if (hh <= median)
 3405         low = ll;
 3406         if (hh >= median)
 3407         high = hh - 1;
 3408     }
 3409 }
 3410 
 3411 #undef ELEM_SWAP
 3412 
 3413 /*
 3414  * angsep_calc - compute angular separation between celestial coordinates
 3415  *   
 3416  * This routine computes the angular separation between to coordinates
 3417  * on the celestial sphere (i.e. RA and Dec).  Note that all units are
 3418  * in DEGREES, unlike the other trig functions in the calculator.
 3419  *
 3420  * double ra1, dec1 - RA and Dec of the first position in degrees
 3421  * double ra2, dec2 - RA and Dec of the second position in degrees
 3422  * 
 3423  * RETURNS: (double) angular separation in degrees
 3424  *
 3425  */
 3426 double angsep_calc(double ra1, double dec1, double ra2, double dec2)
 3427 {
 3428 /*  double cd;  */
 3429   static double deg = 0;
 3430   double a, sdec, sra;
 3431   
 3432   if (deg == 0) deg = ((double)4)*atan((double)1)/((double)180);
 3433   /* deg = 1.0; **** UNCOMMENT IF YOU WANT RADIANS */
 3434 
 3435   /* The algorithm is the law of Haversines.  This algorithm is
 3436      stable even when the points are close together.  The normal
 3437      Law of Cosines fails for angles around 0.1 arcsec. */
 3438 
 3439   sra  = sin( (ra2 - ra1)*deg / 2 );
 3440   sdec = sin( (dec2 - dec1)*deg / 2);
 3441   a = sdec*sdec + cos(dec1*deg)*cos(dec2*deg)*sra*sra;
 3442 
 3443   /* Sanity checking to avoid a range error in the sqrt()'s below */
 3444   if (a < 0) { a = 0; }
 3445   if (a > 1) { a = 1; }
 3446 
 3447   return 2.0*atan2(sqrt(a), sqrt(1.0 - a)) / deg;
 3448 }
 3449 
 3450 static void Do_Func( Node *this )
 3451 {
 3452    Node *theParams[MAXSUBS];
 3453    int  vector[MAXSUBS], allConst;
 3454    lval pVals[MAXSUBS];
 3455    char pNull[MAXSUBS];
 3456    long   ival;
 3457    double dval;
 3458    int  i, valInit;
 3459    long row, elem, nelem;
 3460 
 3461    i = this->nSubNodes;
 3462    allConst = 1;
 3463    while( i-- ) {
 3464       theParams[i] = gParse.Nodes + this->SubNodes[i];
 3465       vector[i]   = ( theParams[i]->operation!=CONST_OP );
 3466       if( vector[i] ) {
 3467      allConst = 0;
 3468      vector[i] = theParams[i]->value.nelem;
 3469       } else {
 3470      if( theParams[i]->type==DOUBLE ) {
 3471         pVals[i].data.dbl = theParams[i]->value.data.dbl;
 3472      } else if( theParams[i]->type==LONG ) {
 3473         pVals[i].data.lng = theParams[i]->value.data.lng;
 3474      } else if( theParams[i]->type==BOOLEAN ) {
 3475         pVals[i].data.log = theParams[i]->value.data.log;
 3476      } else
 3477         strcpy(pVals[i].data.str, theParams[i]->value.data.str);
 3478      pNull[i] = 0;
 3479       }
 3480    }
 3481 
 3482    if( this->nSubNodes==0 ) allConst = 0; /* These do produce scalars */
 3483    /* Random numbers are *never* constant !! */
 3484    if( this->operation == poirnd_fct ) allConst = 0;
 3485    if( this->operation == gasrnd_fct ) allConst = 0;
 3486    if( this->operation == rnd_fct ) allConst = 0;
 3487 
 3488    if( allConst ) {
 3489 
 3490       switch( this->operation ) {
 3491 
 3492         /* Non-Trig single-argument functions */
 3493 
 3494      case sum_fct:
 3495         if( theParams[0]->type==BOOLEAN )
 3496            this->value.data.lng = ( pVals[0].data.log ? 1 : 0 );
 3497         else if( theParams[0]->type==LONG )
 3498            this->value.data.lng = pVals[0].data.lng;
 3499         else if( theParams[0]->type==DOUBLE )
 3500            this->value.data.dbl = pVals[0].data.dbl;
 3501         else if( theParams[0]->type==BITSTR )
 3502           strcpy(this->value.data.str, pVals[0].data.str);
 3503         break;
 3504          case average_fct:
 3505         if( theParams[0]->type==LONG )
 3506            this->value.data.dbl = pVals[0].data.lng;
 3507         else if( theParams[0]->type==DOUBLE )
 3508            this->value.data.dbl = pVals[0].data.dbl;
 3509         break;
 3510          case stddev_fct:
 3511         this->value.data.dbl = 0;  /* Standard deviation of a constant = 0 */
 3512         break;
 3513      case median_fct:
 3514         if( theParams[0]->type==BOOLEAN )
 3515            this->value.data.lng = ( pVals[0].data.log ? 1 : 0 );
 3516         else if( theParams[0]->type==LONG )
 3517            this->value.data.lng = pVals[0].data.lng;
 3518         else
 3519            this->value.data.dbl = pVals[0].data.dbl;
 3520         break;
 3521 
 3522      case poirnd_fct:
 3523         if( theParams[0]->type==DOUBLE )
 3524           this->value.data.lng = simplerng_getpoisson(pVals[0].data.dbl);
 3525         else
 3526           this->value.data.lng = simplerng_getpoisson(pVals[0].data.lng);
 3527         break;
 3528 
 3529      case abs_fct:
 3530         if( theParams[0]->type==DOUBLE ) {
 3531            dval = pVals[0].data.dbl;
 3532            this->value.data.dbl = (dval>0.0 ? dval : -dval);
 3533         } else {
 3534            ival = pVals[0].data.lng;
 3535            this->value.data.lng = (ival> 0  ? ival : -ival);
 3536         }
 3537         break;
 3538 
 3539             /* Special Null-Handling Functions */
 3540 
 3541          case nonnull_fct:
 3542         this->value.data.lng = 1; /* Constants are always 1-element and defined */
 3543         break;
 3544          case isnull_fct:  /* Constants are always defined */
 3545         this->value.data.log = 0;
 3546         break;
 3547          case defnull_fct:
 3548         if( this->type==BOOLEAN )
 3549            this->value.data.log = pVals[0].data.log;
 3550             else if( this->type==LONG )
 3551            this->value.data.lng = pVals[0].data.lng;
 3552             else if( this->type==DOUBLE )
 3553            this->value.data.dbl = pVals[0].data.dbl;
 3554             else if( this->type==STRING )
 3555            strcpy(this->value.data.str,pVals[0].data.str);
 3556         break;
 3557         case setnull_fct: /* Only defined for numeric expressions */
 3558             if( this->type==LONG )
 3559           this->value.data.lng = pVals[0].data.lng;
 3560             else if( this->type==DOUBLE )
 3561            this->value.data.dbl = pVals[0].data.dbl;
 3562         break;
 3563 
 3564         /* Math functions with 1 double argument */
 3565 
 3566      case sin_fct:
 3567         this->value.data.dbl = sin( pVals[0].data.dbl );
 3568         break;
 3569      case cos_fct:
 3570         this->value.data.dbl = cos( pVals[0].data.dbl );
 3571         break;
 3572      case tan_fct:
 3573         this->value.data.dbl = tan( pVals[0].data.dbl );
 3574         break;
 3575      case asin_fct:
 3576         dval = pVals[0].data.dbl;
 3577         if( dval<-1.0 || dval>1.0 )
 3578            yyerror("Out of range argument to arcsin");
 3579         else
 3580            this->value.data.dbl = asin( dval );
 3581         break;
 3582      case acos_fct:
 3583         dval = pVals[0].data.dbl;
 3584         if( dval<-1.0 || dval>1.0 )
 3585            yyerror("Out of range argument to arccos");
 3586         else
 3587            this->value.data.dbl = acos( dval );
 3588         break;
 3589      case atan_fct:
 3590         this->value.data.dbl = atan( pVals[0].data.dbl );
 3591         break;
 3592      case sinh_fct:
 3593         this->value.data.dbl = sinh( pVals[0].data.dbl );
 3594         break;
 3595      case cosh_fct:
 3596         this->value.data.dbl = cosh( pVals[0].data.dbl );
 3597         break;
 3598      case tanh_fct:
 3599         this->value.data.dbl = tanh( pVals[0].data.dbl );
 3600         break;
 3601      case exp_fct:
 3602         this->value.data.dbl = exp( pVals[0].data.dbl );
 3603         break;
 3604      case log_fct:
 3605         dval = pVals[0].data.dbl;
 3606         if( dval<=0.0 )
 3607            yyerror("Out of range argument to log");
 3608         else
 3609            this->value.data.dbl = log( dval );
 3610         break;
 3611      case log10_fct:
 3612         dval = pVals[0].data.dbl;
 3613         if( dval<=0.0 )
 3614            yyerror("Out of range argument to log10");
 3615         else
 3616            this->value.data.dbl = log10( dval );
 3617         break;
 3618      case sqrt_fct:
 3619         dval = pVals[0].data.dbl;
 3620         if( dval<0.0 )
 3621            yyerror("Out of range argument to sqrt");
 3622         else
 3623            this->value.data.dbl = sqrt( dval );
 3624         break;
 3625      case ceil_fct:
 3626         this->value.data.dbl = ceil( pVals[0].data.dbl );
 3627         break;
 3628      case floor_fct:
 3629         this->value.data.dbl = floor( pVals[0].data.dbl );
 3630         break;
 3631      case round_fct:
 3632         this->value.data.dbl = floor( pVals[0].data.dbl + 0.5 );
 3633         break;
 3634 
 3635         /* Two-argument Trig Functions */
 3636 
 3637      case atan2_fct:
 3638         this->value.data.dbl =
 3639            atan2( pVals[0].data.dbl, pVals[1].data.dbl );
 3640         break;
 3641 
 3642         /* Four-argument ANGSEP function */
 3643          case angsep_fct:
 3644         this->value.data.dbl = 
 3645           angsep_calc(pVals[0].data.dbl, pVals[1].data.dbl,
 3646               pVals[2].data.dbl, pVals[3].data.dbl);
 3647 
 3648         /*  Min/Max functions taking 1 or 2 arguments  */
 3649 
 3650          case min1_fct:
 3651         /* No constant vectors! */
 3652         if( this->type == DOUBLE )
 3653            this->value.data.dbl = pVals[0].data.dbl;
 3654         else if( this->type == LONG )
 3655            this->value.data.lng = pVals[0].data.lng;
 3656         else if( this->type == BITSTR )
 3657           strcpy(this->value.data.str, pVals[0].data.str);
 3658         break;
 3659          case min2_fct:
 3660         if( this->type == DOUBLE )
 3661            this->value.data.dbl =
 3662           minvalue( pVals[0].data.dbl, pVals[1].data.dbl );
 3663         else if( this->type == LONG )
 3664            this->value.data.lng =
 3665           minvalue( pVals[0].data.lng, pVals[1].data.lng );
 3666         break;
 3667          case max1_fct:
 3668         /* No constant vectors! */
 3669         if( this->type == DOUBLE )
 3670            this->value.data.dbl = pVals[0].data.dbl;
 3671         else if( this->type == LONG )
 3672            this->value.data.lng = pVals[0].data.lng;
 3673         else if( this->type == BITSTR )
 3674           strcpy(this->value.data.str, pVals[0].data.str);
 3675         break;
 3676          case max2_fct:
 3677         if( this->type == DOUBLE )
 3678            this->value.data.dbl =
 3679           maxvalue( pVals[0].data.dbl, pVals[1].data.dbl );
 3680         else if( this->type == LONG )
 3681            this->value.data.lng =
 3682           maxvalue( pVals[0].data.lng, pVals[1].data.lng );
 3683         break;
 3684 
 3685         /* Boolean SAO region Functions... scalar or vector dbls */
 3686 
 3687      case near_fct:
 3688         this->value.data.log = bnear( pVals[0].data.dbl, pVals[1].data.dbl,
 3689                       pVals[2].data.dbl );
 3690         break;
 3691      case circle_fct:
 3692         this->value.data.log = circle( pVals[0].data.dbl, pVals[1].data.dbl,
 3693                        pVals[2].data.dbl, pVals[3].data.dbl,
 3694                        pVals[4].data.dbl );
 3695         break;
 3696      case box_fct:
 3697         this->value.data.log = saobox( pVals[0].data.dbl, pVals[1].data.dbl,
 3698                        pVals[2].data.dbl, pVals[3].data.dbl,
 3699                        pVals[4].data.dbl, pVals[5].data.dbl,
 3700                        pVals[6].data.dbl );
 3701         break;
 3702      case elps_fct:
 3703         this->value.data.log =
 3704                                ellipse( pVals[0].data.dbl, pVals[1].data.dbl,
 3705                     pVals[2].data.dbl, pVals[3].data.dbl,
 3706                     pVals[4].data.dbl, pVals[5].data.dbl,
 3707                     pVals[6].data.dbl );
 3708         break;
 3709 
 3710             /* C Conditional expression:  bool ? expr : expr */
 3711 
 3712          case ifthenelse_fct:
 3713             switch( this->type ) {
 3714             case BOOLEAN:
 3715                this->value.data.log = ( pVals[2].data.log ?
 3716                                         pVals[0].data.log : pVals[1].data.log );
 3717                break;
 3718             case LONG:
 3719                this->value.data.lng = ( pVals[2].data.log ?
 3720                                         pVals[0].data.lng : pVals[1].data.lng );
 3721                break;
 3722             case DOUBLE:
 3723                this->value.data.dbl = ( pVals[2].data.log ?
 3724                                         pVals[0].data.dbl : pVals[1].data.dbl );
 3725                break;
 3726             case STRING:
 3727            strcpy(this->value.data.str, ( pVals[2].data.log ?
 3728                                               pVals[0].data.str :
 3729                                               pVals[1].data.str ) );
 3730                break;
 3731             }
 3732             break;
 3733 
 3734         /* String functions */
 3735          case strmid_fct:
 3736        cstrmid(this->value.data.str, this->value.nelem, 
 3737            pVals[0].data.str,    pVals[0].nelem,
 3738            pVals[1].data.lng);
 3739        break;
 3740          case strpos_fct:
 3741        {
 3742          char *res = strstr(pVals[0].data.str, pVals[1].data.str);
 3743          if (res == NULL) {
 3744            this->value.data.lng = 0; 
 3745          } else {
 3746            this->value.data.lng = (res - pVals[0].data.str) + 1;
 3747          }
 3748          break;
 3749        }
 3750 
 3751       }
 3752       this->operation = CONST_OP;
 3753 
 3754    } else {
 3755 
 3756       Allocate_Ptrs( this );
 3757 
 3758       row  = gParse.nRows;
 3759       elem = row * this->value.nelem;
 3760 
 3761       if( !gParse.status ) {
 3762      switch( this->operation ) {
 3763 
 3764         /* Special functions with no arguments */
 3765 
 3766      case row_fct:
 3767         while( row-- ) {
 3768            this->value.data.lngptr[row] = gParse.firstRow + row;
 3769            this->value.undef[row] = 0;
 3770         }
 3771         break;
 3772      case null_fct:
 3773             if( this->type==LONG ) {
 3774                while( row-- ) {
 3775                   this->value.data.lngptr[row] = 0;
 3776                   this->value.undef[row] = 1;
 3777                }
 3778             } else if( this->type==STRING ) {
 3779                while( row-- ) {
 3780                   this->value.data.strptr[row][0] = '\0';
 3781                   this->value.undef[row] = 1;
 3782                }
 3783             }
 3784         break;
 3785      case rnd_fct:
 3786        while( elem-- ) {
 3787          this->value.data.dblptr[elem] = simplerng_getuniform();
 3788          this->value.undef[elem] = 0;
 3789         }
 3790         break;
 3791 
 3792      case gasrnd_fct:
 3793         while( elem-- ) {
 3794            this->value.data.dblptr[elem] = simplerng_getnorm();
 3795            this->value.undef[elem] = 0;
 3796         }
 3797         break;
 3798 
 3799      case poirnd_fct:
 3800        if( theParams[0]->type==DOUBLE ) {
 3801           if (theParams[0]->operation == CONST_OP) {
 3802         while( elem-- ) {
 3803           this->value.undef[elem] = (pVals[0].data.dbl < 0);
 3804           if (! this->value.undef[elem]) {
 3805             this->value.data.lngptr[elem] = simplerng_getpoisson(pVals[0].data.dbl);
 3806           }
 3807         } 
 3808           } else {
 3809         while( elem-- ) {
 3810           this->value.undef[elem] = theParams[0]->value.undef[elem];
 3811           if (theParams[0]->value.data.dblptr[elem] < 0) 
 3812             this->value.undef[elem] = 1;
 3813           if (! this->value.undef[elem]) {
 3814             this->value.data.lngptr[elem] = 
 3815               simplerng_getpoisson(theParams[0]->value.data.dblptr[elem]);
 3816           }
 3817         } /* while */
 3818           } /* ! CONST_OP */
 3819        } else {
 3820          /* LONG */
 3821           if (theParams[0]->operation == CONST_OP) {
 3822         while( elem-- ) {
 3823           this->value.undef[elem] = (pVals[0].data.lng < 0);
 3824           if (! this->value.undef[elem]) {
 3825             this->value.data.lngptr[elem] = simplerng_getpoisson(pVals[0].data.lng);
 3826           }
 3827         } 
 3828           } else {
 3829         while( elem-- ) {
 3830           this->value.undef[elem] = theParams[0]->value.undef[elem];
 3831           if (theParams[0]->value.data.lngptr[elem] < 0) 
 3832             this->value.undef[elem] = 1;
 3833           if (! this->value.undef[elem]) {
 3834             this->value.data.lngptr[elem] = 
 3835               simplerng_getpoisson(theParams[0]->value.data.lngptr[elem]);
 3836           }
 3837         } /* while */
 3838           } /* ! CONST_OP */
 3839        } /* END LONG */
 3840        break;
 3841 
 3842 
 3843         /* Non-Trig single-argument functions */
 3844         
 3845      case sum_fct:
 3846         elem = row * theParams[0]->value.nelem;
 3847         if( theParams[0]->type==BOOLEAN ) {
 3848            while( row-- ) {
 3849           this->value.data.lngptr[row] = 0;
 3850           /* Default is UNDEF until a defined value is found */
 3851           this->value.undef[row] = 1;
 3852           nelem = theParams[0]->value.nelem;
 3853           while( nelem-- ) {
 3854              elem--;
 3855              if ( ! theParams[0]->value.undef[elem] ) {
 3856                this->value.data.lngptr[row] +=
 3857              ( theParams[0]->value.data.logptr[elem] ? 1 : 0 );
 3858                this->value.undef[row] = 0;
 3859              }
 3860           }
 3861            }
 3862         } else if( theParams[0]->type==LONG ) {
 3863            while( row-- ) {
 3864           this->value.data.lngptr[row] = 0;
 3865           /* Default is UNDEF until a defined value is found */
 3866           this->value.undef[row] = 1;
 3867           nelem = theParams[0]->value.nelem;
 3868           while( nelem-- ) {
 3869              elem--;
 3870              if ( ! theParams[0]->value.undef[elem] ) {
 3871                this->value.data.lngptr[row] +=
 3872              theParams[0]->value.data.lngptr[elem];
 3873                this->value.undef[row] = 0;
 3874              }
 3875           }
 3876            }          
 3877         } else if( theParams[0]->type==DOUBLE ){
 3878            while( row-- ) {
 3879           this->value.data.dblptr[row] = 0.0;
 3880           /* Default is UNDEF until a defined value is found */
 3881           this->value.undef[row] = 1;
 3882           nelem = theParams[0]->value.nelem;
 3883           while( nelem-- ) {
 3884              elem--;
 3885              if ( ! theParams[0]->value.undef[elem] ) {
 3886                this->value.data.dblptr[row] +=
 3887              theParams[0]->value.data.dblptr[elem];
 3888                this->value.undef[row] = 0;
 3889              }
 3890           }
 3891            }          
 3892         } else { /* BITSTR */
 3893            nelem = theParams[0]->value.nelem;
 3894            while( row-- ) {
 3895           char *sptr1 = theParams[0]->value.data.strptr[row];
 3896           this->value.data.lngptr[row] = 0;
 3897           this->value.undef[row] = 0;
 3898           while (*sptr1) {
 3899             if (*sptr1 == '1') this->value.data.lngptr[row] ++;
 3900             sptr1++;
 3901           }
 3902            }          
 3903         }
 3904         break;
 3905 
 3906      case average_fct:
 3907         elem = row * theParams[0]->value.nelem;
 3908         if( theParams[0]->type==LONG ) {
 3909            while( row-- ) {
 3910           int count = 0;
 3911           this->value.data.dblptr[row] = 0;
 3912           nelem = theParams[0]->value.nelem;
 3913           while( nelem-- ) {
 3914              elem--;
 3915              if (theParams[0]->value.undef[elem] == 0) {
 3916                this->value.data.dblptr[row] +=
 3917              theParams[0]->value.data.lngptr[elem];
 3918                count ++;
 3919              }
 3920           }
 3921           if (count == 0) {
 3922             this->value.undef[row] = 1;
 3923           } else {
 3924             this->value.undef[row] = 0;
 3925             this->value.data.dblptr[row] /= count;
 3926           }
 3927            }          
 3928         } else if( theParams[0]->type==DOUBLE ){
 3929            while( row-- ) {
 3930           int count = 0;
 3931           this->value.data.dblptr[row] = 0;
 3932           nelem = theParams[0]->value.nelem;
 3933           while( nelem-- ) {
 3934              elem--;
 3935              if (theParams[0]->value.undef[elem] == 0) {
 3936                this->value.data.dblptr[row] +=
 3937              theParams[0]->value.data.dblptr[elem];
 3938                count ++;
 3939              }
 3940           }
 3941           if (count == 0) {
 3942             this->value.undef[row] = 1;
 3943           } else {
 3944             this->value.undef[row] = 0;
 3945             this->value.data.dblptr[row] /= count;
 3946           }
 3947            }          
 3948         }
 3949         break;
 3950      case stddev_fct:
 3951         elem = row * theParams[0]->value.nelem;
 3952         if( theParams[0]->type==LONG ) {
 3953 
 3954            /* Compute the mean value */
 3955            while( row-- ) {
 3956           int count = 0;
 3957           double sum = 0, sum2 = 0;
 3958 
 3959           nelem = theParams[0]->value.nelem;
 3960           while( nelem-- ) {
 3961              elem--;
 3962              if (theParams[0]->value.undef[elem] == 0) {
 3963                sum += theParams[0]->value.data.lngptr[elem];
 3964                count ++;
 3965              }
 3966           }
 3967           if (count > 1) {
 3968             sum /= count;
 3969 
 3970             /* Compute the sum of squared deviations */
 3971             nelem = theParams[0]->value.nelem;
 3972             elem += nelem;  /* Reset elem for second pass */
 3973             while( nelem-- ) {
 3974               elem--;
 3975               if (theParams[0]->value.undef[elem] == 0) {
 3976             double dx = (theParams[0]->value.data.lngptr[elem] - sum);
 3977             sum2 += (dx*dx);
 3978               }
 3979             }
 3980 
 3981             sum2 /= (double)count-1;
 3982 
 3983             this->value.undef[row] = 0;
 3984             this->value.data.dblptr[row] = sqrt(sum2);
 3985           } else {
 3986             this->value.undef[row] = 0;       /* STDDEV => 0 */
 3987             this->value.data.dblptr[row] = 0;
 3988           }
 3989            }
 3990         } else if( theParams[0]->type==DOUBLE ){
 3991 
 3992            /* Compute the mean value */
 3993            while( row-- ) {
 3994           int count = 0;
 3995           double sum = 0, sum2 = 0;
 3996 
 3997           nelem = theParams[0]->value.nelem;
 3998           while( nelem-- ) {
 3999              elem--;
 4000              if (theParams[0]->value.undef[elem] == 0) {
 4001                sum += theParams[0]->value.data.dblptr[elem];
 4002                count ++;
 4003              }
 4004           }
 4005           if (count > 1) {
 4006             sum /= count;
 4007 
 4008             /* Compute the sum of squared deviations */
 4009             nelem = theParams[0]->value.nelem;
 4010             elem += nelem;  /* Reset elem for second pass */
 4011             while( nelem-- ) {
 4012               elem--;
 4013               if (theParams[0]->value.undef[elem] == 0) {
 4014             double dx = (theParams[0]->value.data.dblptr[elem] - sum);
 4015             sum2 += (dx*dx);
 4016               }
 4017             }
 4018 
 4019             sum2 /= (double)count-1;
 4020 
 4021             this->value.undef[row] = 0;
 4022             this->value.data.dblptr[row] = sqrt(sum2);
 4023           } else {
 4024             this->value.undef[row] = 0;       /* STDDEV => 0 */
 4025             this->value.data.dblptr[row] = 0;
 4026           }
 4027            }
 4028         }
 4029         break;
 4030 
 4031      case median_fct:
 4032        elem = row * theParams[0]->value.nelem;
 4033        nelem = theParams[0]->value.nelem;
 4034        if( theParams[0]->type==LONG ) {
 4035            long *dptr = theParams[0]->value.data.lngptr;
 4036            char *uptr = theParams[0]->value.undef;
 4037            long *mptr = (long *) malloc(sizeof(long)*nelem);
 4038            int irow;
 4039 
 4040            /* Allocate temporary storage for this row, since the
 4041                   quickselect function will scramble the contents */
 4042            if (mptr == 0) {
 4043          yyerror("Could not allocate temporary memory in median function");
 4044          free( this->value.data.ptr );
 4045          break;
 4046            }
 4047 
 4048            for (irow=0; irow<row; irow++) {
 4049           long *p = mptr;
 4050           int nelem1 = nelem;
 4051 
 4052 
 4053           while ( nelem1-- ) { 
 4054             if (*uptr == 0) {
 4055               *p++ = *dptr;   /* Only advance the dest pointer if we copied */
 4056             }
 4057             dptr ++;  /* Advance the source pointer ... */
 4058             uptr ++;  /* ... and source "undef" pointer */
 4059           }
 4060           
 4061           nelem1 = (p - mptr);  /* Number of accepted data points */
 4062           if (nelem1 > 0) {
 4063             this->value.undef[irow] = 0;
 4064             this->value.data.lngptr[irow] = qselect_median_lng(mptr, nelem1);
 4065           } else {
 4066             this->value.undef[irow] = 1;
 4067             this->value.data.lngptr[irow] = 0;
 4068           }
 4069             
 4070            }          
 4071 
 4072            free(mptr);
 4073         } else {
 4074            double *dptr = theParams[0]->value.data.dblptr;
 4075            char   *uptr = theParams[0]->value.undef;
 4076            double *mptr = (double *) malloc(sizeof(double)*nelem);
 4077            int irow;
 4078 
 4079            /* Allocate temporary storage for this row, since the
 4080                   quickselect function will scramble the contents */
 4081            if (mptr == 0) {
 4082          yyerror("Could not allocate temporary memory in median function");
 4083          free( this->value.data.ptr );
 4084          break;
 4085            }
 4086 
 4087            for (irow=0; irow<row; irow++) {
 4088           double *p = mptr;
 4089           int nelem1 = nelem;
 4090 
 4091           while ( nelem1-- ) { 
 4092             if (*uptr == 0) {
 4093               *p++ = *dptr;   /* Only advance the dest pointer if we copied */
 4094             }
 4095             dptr ++;  /* Advance the source pointer ... */
 4096             uptr ++;  /* ... and source "undef" pointer */
 4097           }
 4098 
 4099           nelem1 = (p - mptr);  /* Number of accepted data points */
 4100           if (nelem1 > 0) {
 4101             this->value.undef[irow] = 0;
 4102             this->value.data.dblptr[irow] = qselect_median_dbl(mptr, nelem1);
 4103           } else {
 4104             this->value.undef[irow] = 1;
 4105             this->value.data.dblptr[irow] = 0;
 4106           }
 4107 
 4108            }
 4109            free(mptr);
 4110         }
 4111         break;
 4112      case abs_fct:
 4113         if( theParams[0]->type==DOUBLE )
 4114            while( elem-- ) {
 4115           dval = theParams[0]->value.data.dblptr[elem];
 4116           this->value.data.dblptr[elem] = (dval>0.0 ? dval : -dval);
 4117           this->value.undef[elem] = theParams[0]->value.undef[elem];
 4118            }
 4119         else
 4120            while( elem-- ) {
 4121           ival = theParams[0]->value.data.lngptr[elem];
 4122           this->value.data.lngptr[elem] = (ival> 0  ? ival : -ival);
 4123           this->value.undef[elem] = theParams[0]->value.undef[elem];
 4124            }
 4125         break;
 4126 
 4127             /* Special Null-Handling Functions */
 4128 
 4129      case nonnull_fct:
 4130        nelem = theParams[0]->value.nelem;
 4131        if ( theParams[0]->type==STRING ) nelem = 1;
 4132        elem = row * nelem;
 4133        while( row-- ) {
 4134          int nelem1 = nelem;
 4135 
 4136          this->value.undef[row] = 0;        /* Initialize to 0 (defined) */
 4137          this->value.data.lngptr[row] = 0;
 4138          while( nelem1-- ) {    
 4139            elem --;
 4140            if ( theParams[0]->value.undef[elem] == 0 ) this->value.data.lngptr[row] ++;
 4141          }
 4142        }
 4143        break;
 4144      case isnull_fct:
 4145         if( theParams[0]->type==STRING ) elem = row;
 4146         while( elem-- ) {
 4147            this->value.data.logptr[elem] = theParams[0]->value.undef[elem];
 4148            this->value.undef[elem] = 0;
 4149         }
 4150         break;
 4151          case defnull_fct:
 4152         switch( this->type ) {
 4153         case BOOLEAN:
 4154            while( row-- ) {
 4155           nelem = this->value.nelem;
 4156           while( nelem-- ) {
 4157              elem--;
 4158              i=2; while( i-- )
 4159             if( vector[i]>1 ) {
 4160                pNull[i] = theParams[i]->value.undef[elem];
 4161                pVals[i].data.log =
 4162                   theParams[i]->value.data.logptr[elem];
 4163             } else if( vector[i] ) {
 4164                pNull[i] = theParams[i]->value.undef[row];
 4165                pVals[i].data.log =
 4166                   theParams[i]->value.data.logptr[row];
 4167             }
 4168              if( pNull[0] ) {
 4169             this->value.undef[elem] = pNull[1];
 4170             this->value.data.logptr[elem] = pVals[1].data.log;
 4171              } else {
 4172             this->value.undef[elem] = 0;
 4173             this->value.data.logptr[elem] = pVals[0].data.log;
 4174              }
 4175           }
 4176            }
 4177            break;
 4178         case LONG:
 4179            while( row-- ) {
 4180           nelem = this->value.nelem;
 4181           while( nelem-- ) {
 4182              elem--;
 4183              i=2; while( i-- )
 4184             if( vector[i]>1 ) {
 4185                pNull[i] = theParams[i]->value.undef[elem];
 4186                pVals[i].data.lng =
 4187                   theParams[i]->value.data.lngptr[elem];
 4188             } else if( vector[i] ) {
 4189                pNull[i] = theParams[i]->value.undef[row];
 4190                pVals[i].data.lng =
 4191                   theParams[i]->value.data.lngptr[row];
 4192             }
 4193              if( pNull[0] ) {
 4194             this->value.undef[elem] = pNull[1];
 4195             this->value.data.lngptr[elem] = pVals[1].data.lng;
 4196              } else {
 4197             this->value.undef[elem] = 0;
 4198             this->value.data.lngptr[elem] = pVals[0].data.lng;
 4199              }
 4200           }
 4201            }
 4202            break;
 4203         case DOUBLE:
 4204            while( row-- ) {
 4205           nelem = this->value.nelem;
 4206           while( nelem-- ) {
 4207              elem--;
 4208              i=2; while( i-- )
 4209             if( vector[i]>1 ) {
 4210                pNull[i] = theParams[i]->value.undef[elem];
 4211                pVals[i].data.dbl =
 4212                   theParams[i]->value.data.dblptr[elem];
 4213             } else if( vector[i] ) {
 4214                pNull[i] = theParams[i]->value.undef[row];
 4215                pVals[i].data.dbl =
 4216                   theParams[i]->value.data.dblptr[row];
 4217             }
 4218              if( pNull[0] ) {
 4219             this->value.undef[elem] = pNull[1];
 4220             this->value.data.dblptr[elem] = pVals[1].data.dbl;
 4221              } else {
 4222             this->value.undef[elem] = 0;
 4223             this->value.data.dblptr[elem] = pVals[0].data.dbl;
 4224              }
 4225           }
 4226            }
 4227            break;
 4228         case STRING:
 4229            while( row-- ) {
 4230           i=2; while( i-- )
 4231              if( vector[i] ) {
 4232             pNull[i] = theParams[i]->value.undef[row];
 4233             strcpy(pVals[i].data.str,
 4234                    theParams[i]->value.data.strptr[row]);
 4235              }
 4236           if( pNull[0] ) {
 4237              this->value.undef[row] = pNull[1];
 4238              strcpy(this->value.data.strptr[row],pVals[1].data.str);
 4239           } else {
 4240              this->value.undef[elem] = 0;
 4241              strcpy(this->value.data.strptr[row],pVals[0].data.str);
 4242           }
 4243            }
 4244         }
 4245         break;
 4246          case setnull_fct:
 4247         switch( this->type ) {
 4248         case LONG:
 4249           while( elem-- ) {
 4250         if ( theParams[1]->value.data.lng == 
 4251              theParams[0]->value.data.lngptr[elem] ) {
 4252           this->value.data.lngptr[elem] = 0;
 4253           this->value.undef[elem] = 1;
 4254         } else {
 4255           this->value.data.lngptr[elem] = theParams[0]->value.data.lngptr[elem];
 4256           this->value.undef[elem] = theParams[0]->value.undef[elem];
 4257         }
 4258           }
 4259           break;
 4260         case DOUBLE:
 4261           while( elem-- ) {
 4262         if ( theParams[1]->value.data.dbl == 
 4263              theParams[0]->value.data.dblptr[elem] ) {
 4264           this->value.data.dblptr[elem] = 0;
 4265           this->value.undef[elem] = 1;
 4266         } else {
 4267           this->value.data.dblptr[elem] = theParams[0]->value.data.dblptr[elem];
 4268           this->value.undef[elem] = theParams[0]->value.undef[elem];
 4269         }
 4270           }
 4271           break;
 4272         }
 4273         break;
 4274 
 4275         /* Math functions with 1 double argument */
 4276 
 4277      case sin_fct:
 4278         while( elem-- )
 4279            if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
 4280           this->value.data.dblptr[elem] = 
 4281              sin( theParams[0]->value.data.dblptr[elem] );
 4282            }
 4283         break;
 4284      case cos_fct:
 4285         while( elem-- )
 4286            if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
 4287           this->value.data.dblptr[elem] = 
 4288              cos( theParams[0]->value.data.dblptr[elem] );
 4289            }
 4290         break;
 4291      case tan_fct:
 4292         while( elem-- )
 4293            if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
 4294           this->value.data.dblptr[elem] = 
 4295              tan( theParams[0]->value.data.dblptr[elem] );
 4296            }
 4297         break;
 4298      case asin_fct:
 4299         while( elem-- )
 4300            if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
 4301           dval = theParams[0]->value.data.dblptr[elem];
 4302           if( dval<-1.0 || dval>1.0 ) {
 4303              this->value.data.dblptr[elem] = 0.0;
 4304              this->value.undef[elem] = 1;
 4305           } else
 4306              this->value.data.dblptr[elem] = asin( dval );
 4307            }
 4308         break;
 4309      case acos_fct:
 4310         while( elem-- )
 4311            if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
 4312           dval = theParams[0]->value.data.dblptr[elem];
 4313           if( dval<-1.0 || dval>1.0 ) {
 4314              this->value.data.dblptr[elem] = 0.0;
 4315              this->value.undef[elem] = 1;
 4316           } else
 4317              this->value.data.dblptr[elem] = acos( dval );
 4318            }
 4319         break;
 4320      case atan_fct:
 4321         while( elem-- )
 4322            if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
 4323           dval = theParams[0]->value.data.dblptr[elem];
 4324           this->value.data.dblptr[elem] = atan( dval );
 4325            }
 4326         break;
 4327      case sinh_fct:
 4328         while( elem-- )
 4329            if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
 4330           this->value.data.dblptr[elem] = 
 4331              sinh( theParams[0]->value.data.dblptr[elem] );
 4332            }
 4333         break;
 4334      case cosh_fct:
 4335         while( elem-- )
 4336            if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
 4337           this->value.data.dblptr[elem] = 
 4338              cosh( theParams[0]->value.data.dblptr[elem] );
 4339            }
 4340         break;
 4341      case tanh_fct:
 4342         while( elem-- )
 4343            if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
 4344           this->value.data.dblptr[elem] = 
 4345              tanh( theParams[0]->value.data.dblptr[elem] );
 4346            }
 4347         break;
 4348      case exp_fct:
 4349         while( elem-- )
 4350            if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
 4351           dval = theParams[0]->value.data.dblptr[elem];
 4352           this->value.data.dblptr[elem] = exp( dval );
 4353            }
 4354         break;
 4355      case log_fct:
 4356         while( elem-- )
 4357            if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
 4358           dval = theParams[0]->value.data.dblptr[elem];
 4359           if( dval<=0.0 ) {
 4360              this->value.data.dblptr[elem] = 0.0;
 4361              this->value.undef[elem] = 1;
 4362           } else
 4363              this->value.data.dblptr[elem] = log( dval );
 4364            }
 4365         break;
 4366      case log10_fct:
 4367         while( elem-- )
 4368            if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
 4369           dval = theParams[0]->value.data.dblptr[elem];
 4370           if( dval<=0.0 ) {
 4371              this->value.data.dblptr[elem] = 0.0;
 4372              this->value.undef[elem] = 1;
 4373           } else
 4374              this->value.data.dblptr[elem] = log10( dval );
 4375            }
 4376         break;
 4377      case sqrt_fct:
 4378         while( elem-- )
 4379            if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
 4380           dval = theParams[0]->value.data.dblptr[elem];
 4381           if( dval<0.0 ) {
 4382              this->value.data.dblptr[elem] = 0.0;
 4383              this->value.undef[elem] = 1;
 4384           } else
 4385              this->value.data.dblptr[elem] = sqrt( dval );
 4386            }
 4387         break;
 4388      case ceil_fct:
 4389         while( elem-- )
 4390            if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
 4391           this->value.data.dblptr[elem] = 
 4392              ceil( theParams[0]->value.data.dblptr[elem] );
 4393            }
 4394         break;
 4395      case floor_fct:
 4396         while( elem-- )
 4397            if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
 4398           this->value.data.dblptr[elem] = 
 4399              floor( theParams[0]->value.data.dblptr[elem] );
 4400            }
 4401         break;
 4402      case round_fct:
 4403         while( elem-- )
 4404            if( !(this->value.undef[elem] = theParams[0]->value.undef[elem]) ) {
 4405           this->value.data.dblptr[elem] = 
 4406              floor( theParams[0]->value.data.dblptr[elem] + 0.5);
 4407            }
 4408         break;
 4409 
 4410         /* Two-argument Trig Functions */
 4411         
 4412      case atan2_fct:
 4413         while( row-- ) {
 4414            nelem = this->value.nelem;
 4415            while( nelem-- ) {
 4416           elem--;
 4417           i=2; while( i-- )
 4418              if( vector[i]>1 ) {
 4419             pVals[i].data.dbl =
 4420                theParams[i]->value.data.dblptr[elem];
 4421             pNull[i] = theParams[i]->value.undef[elem];
 4422              } else if( vector[i] ) {
 4423             pVals[i].data.dbl =
 4424                theParams[i]->value.data.dblptr[row];
 4425             pNull[i] = theParams[i]->value.undef[row];
 4426              }
 4427           if( !(this->value.undef[elem] = (pNull[0] || pNull[1]) ) )
 4428              this->value.data.dblptr[elem] =
 4429             atan2( pVals[0].data.dbl, pVals[1].data.dbl );
 4430            }
 4431         }
 4432         break;
 4433 
 4434         /* Four-argument ANGSEP Function */
 4435         
 4436      case angsep_fct:
 4437         while( row-- ) {
 4438            nelem = this->value.nelem;
 4439            while( nelem-- ) {
 4440           elem--;
 4441           i=4; while( i-- )
 4442              if( vector[i]>1 ) {
 4443             pVals[i].data.dbl =
 4444                theParams[i]->value.data.dblptr[elem];
 4445             pNull[i] = theParams[i]->value.undef[elem];
 4446              } else if( vector[i] ) {
 4447             pVals[i].data.dbl =
 4448                theParams[i]->value.data.dblptr[row];
 4449             pNull[i] = theParams[i]->value.undef[row];
 4450              }
 4451           if( !(this->value.undef[elem] = (pNull[0] || pNull[1] ||
 4452                            pNull[2] || pNull[3]) ) )
 4453              this->value.data.dblptr[elem] =
 4454                angsep_calc(pVals[0].data.dbl, pVals[1].data.dbl,
 4455                    pVals[2].data.dbl, pVals[3].data.dbl);
 4456            }
 4457         }
 4458         break;
 4459 
 4460 
 4461 
 4462         /*  Min/Max functions taking 1 or 2 arguments  */
 4463 
 4464          case min1_fct:
 4465         elem = row * theParams[0]->value.nelem;
 4466         if( this->type==LONG ) {
 4467            long minVal=0;
 4468            while( row-- ) {
 4469           valInit = 1;
 4470           this->value.undef[row] = 1;
 4471           nelem = theParams[0]->value.nelem;
 4472           while( nelem-- ) {
 4473              elem--;
 4474              if ( !theParams[0]->value.undef[elem] ) {
 4475                if ( valInit ) {
 4476              valInit = 0;
 4477              minVal  = theParams[0]->value.data.lngptr[elem];
 4478                } else {
 4479              minVal  = minvalue( minVal,
 4480                          theParams[0]->value.data.lngptr[elem] );
 4481                }
 4482                this->value.undef[row] = 0;
 4483              }
 4484           }  
 4485           this->value.data.lngptr[row] = minVal;
 4486            }          
 4487         } else if( this->type==DOUBLE ) {
 4488            double minVal=0.0;
 4489            while( row-- ) {
 4490           valInit = 1;
 4491           this->value.undef[row] = 1;
 4492           nelem = theParams[0]->value.nelem;
 4493           while( nelem-- ) {
 4494              elem--;
 4495              if ( !theParams[0]->value.undef[elem] ) {
 4496                if ( valInit ) {
 4497              valInit = 0;
 4498              minVal  = theParams[0]->value.data.dblptr[elem];
 4499                } else {
 4500              minVal  = minvalue( minVal,
 4501                          theParams[0]->value.data.dblptr[elem] );
 4502                }
 4503                this->value.undef[row] = 0;
 4504              }
 4505           }  
 4506           this->value.data.dblptr[row] = minVal;
 4507            }          
 4508         } else if( this->type==BITSTR ) {
 4509            char minVal;
 4510            while( row-- ) {
 4511           char *sptr1 = theParams[0]->value.data.strptr[row];
 4512           minVal = '1';
 4513           while (*sptr1) {
 4514             if (*sptr1 == '0') minVal = '0';
 4515             sptr1++;
 4516           }
 4517           this->value.data.strptr[row][0] = minVal;
 4518           this->value.data.strptr[row][1] = 0;     /* Null terminate */
 4519            }          
 4520         }
 4521         break;
 4522          case min2_fct:
 4523         if( this->type==LONG ) {
 4524            while( row-- ) {
 4525           nelem = this->value.nelem;
 4526           while( nelem-- ) {
 4527              elem--;
 4528              i=2; while( i-- )
 4529             if( vector[i]>1 ) {
 4530                pVals[i].data.lng =
 4531                   theParams[i]->value.data.lngptr[elem];
 4532                pNull[i] = theParams[i]->value.undef[elem];
 4533             } else if( vector[i] ) {
 4534                pVals[i].data.lng =
 4535                   theParams[i]->value.data.lngptr[row];
 4536                pNull[i] = theParams[i]->value.undef[row];
 4537             }
 4538              if( pNull[0] && pNull[1] ) {
 4539                this->value.undef[elem] = 1;
 4540                this->value.data.lngptr[elem] = 0;
 4541              } else if (pNull[0]) {
 4542                this->value.undef[elem] = 0;
 4543                this->value.data.lngptr[elem] = pVals[1].data.lng;
 4544              } else if (pNull[1]) {
 4545                this->value.undef[elem] = 0;
 4546                this->value.data.lngptr[elem] = pVals[0].data.lng;
 4547              } else {
 4548                this->value.undef[elem] = 0;
 4549                this->value.data.lngptr[elem] =
 4550              minvalue( pVals[0].data.lng, pVals[1].data.lng );
 4551              }
 4552           }
 4553            }
 4554         } else if( this->type==DOUBLE ) {
 4555            while( row-- ) {
 4556           nelem = this->value.nelem;
 4557           while( nelem-- ) {
 4558              elem--;
 4559              i=2; while( i-- )
 4560             if( vector[i]>1 ) {
 4561                pVals[i].data.dbl =
 4562                   theParams[i]->value.data.dblptr[elem];
 4563                pNull[i] = theParams[i]->value.undef[elem];
 4564             } else if( vector[i] ) {
 4565                pVals[i].data.dbl =
 4566                   theParams[i]->value.data.dblptr[row];
 4567                pNull[i] = theParams[i]->value.undef[row];
 4568             }
 4569              if( pNull[0] && pNull[1] ) {
 4570                this->value.undef[elem] = 1;
 4571                this->value.data.dblptr[elem] = 0;
 4572              } else if (pNull[0]) {
 4573                this->value.undef[elem] = 0;
 4574                this->value.data.dblptr[elem] = pVals[1].data.dbl;
 4575              } else if (pNull[1]) {
 4576                this->value.undef[elem] = 0;
 4577                this->value.data.dblptr[elem] = pVals[0].data.dbl;
 4578              } else {
 4579                this->value.undef[elem] = 0;
 4580                this->value.data.dblptr[elem] =
 4581              minvalue( pVals[0].data.dbl, pVals[1].data.dbl );
 4582              }
 4583           }
 4584            }
 4585         }
 4586         break;
 4587 
 4588          case max1_fct:
 4589         elem = row * theParams[0]->value.nelem;
 4590         if( this->type==LONG ) {
 4591            long maxVal=0;
 4592            while( row-- ) {
 4593           valInit = 1;
 4594           this->value.undef[row] = 1;
 4595           nelem = theParams[0]->value.nelem;
 4596           while( nelem-- ) {
 4597              elem--;
 4598              if ( !theParams[0]->value.undef[elem] ) {
 4599                if ( valInit ) {
 4600              valInit = 0;
 4601              maxVal  = theParams[0]->value.data.lngptr[elem];
 4602                } else {
 4603              maxVal  = maxvalue( maxVal,
 4604                          theParams[0]->value.data.lngptr[elem] );
 4605                }
 4606                this->value.undef[row] = 0;
 4607              }
 4608           }
 4609           this->value.data.lngptr[row] = maxVal;
 4610            }          
 4611         } else if( this->type==DOUBLE ) {
 4612            double maxVal=0.0;
 4613            while( row-- ) {
 4614           valInit = 1;
 4615           this->value.undef[row] = 1;
 4616           nelem = theParams[0]->value.nelem;
 4617           while( nelem-- ) {
 4618              elem--;
 4619              if ( !theParams[0]->value.undef[elem] ) {
 4620                if ( valInit ) {
 4621              valInit = 0;
 4622              maxVal  = theParams[0]->value.data.dblptr[elem];
 4623                } else {
 4624              maxVal  = maxvalue( maxVal,
 4625                          theParams[0]->value.data.dblptr[elem] );
 4626                }
 4627                this->value.undef[row] = 0;
 4628              }
 4629           }
 4630           this->value.data.dblptr[row] = maxVal;
 4631            }          
 4632         } else if( this->type==BITSTR ) {
 4633            char maxVal;
 4634            while( row-- ) {
 4635           char *sptr1 = theParams[0]->value.data.strptr[row];
 4636           maxVal = '0';
 4637           while (*sptr1) {
 4638             if (*sptr1 == '1') maxVal = '1';
 4639             sptr1++;
 4640           }
 4641           this->value.data.strptr[row][0] = maxVal;
 4642           this->value.data.strptr[row][1] = 0;     /* Null terminate */
 4643            }          
 4644         }
 4645         break;
 4646          case max2_fct:
 4647         if( this->type==LONG ) {
 4648            while( row-- ) {
 4649           nelem = this->value.nelem;
 4650           while( nelem-- ) {
 4651              elem--;
 4652              i=2; while( i-- )
 4653             if( vector[i]>1 ) {
 4654                pVals[i].data.lng =
 4655                   theParams[i]->value.data.lngptr[elem];
 4656                pNull[i] = theParams[i]->value.undef[elem];
 4657             } else if( vector[i] ) {
 4658                pVals[i].data.lng =
 4659                   theParams[i]->value.data.lngptr[row];
 4660                pNull[i] = theParams[i]->value.undef[row];
 4661             }
 4662              if( pNull[0] && pNull[1] ) {
 4663                this->value.undef[elem] = 1;
 4664                this->value.data.lngptr[elem] = 0;
 4665              } else if (pNull[0]) {
 4666                this->value.undef[elem] = 0;
 4667                this->value.data.lngptr[elem] = pVals[1].data.lng;
 4668              } else if (pNull[1]) {
 4669                this->value.undef[elem] = 0;
 4670                this->value.data.lngptr[elem] = pVals[0].data.lng;
 4671              } else {
 4672                this->value.undef[elem] = 0;
 4673                this->value.data.lngptr[elem] =
 4674              maxvalue( pVals[0].data.lng, pVals[1].data.lng );
 4675              }
 4676           }
 4677            }
 4678         } else if( this->type==DOUBLE ) {
 4679            while( row-- ) {
 4680           nelem = this->value.nelem;
 4681           while( nelem-- ) {
 4682              elem--;
 4683              i=2; while( i-- )
 4684             if( vector[i]>1 ) {
 4685                pVals[i].data.dbl =
 4686                   theParams[i]->value.data.dblptr[elem];
 4687                pNull[i] = theParams[i]->value.undef[elem];
 4688             } else if( vector[i] ) {
 4689                pVals[i].data.dbl =
 4690                   theParams[i]->value.data.dblptr[row];
 4691                pNull[i] = theParams[i]->value.undef[row];
 4692             }
 4693              if( pNull[0] && pNull[1] ) {
 4694                this->value.undef[elem] = 1;
 4695                this->value.data.dblptr[elem] = 0;
 4696              } else if (pNull[0]) {
 4697                this->value.undef[elem] = 0;
 4698                this->value.data.dblptr[elem] = pVals[1].data.dbl;
 4699              } else if (pNull[1]) {
 4700                this->value.undef[elem] = 0;
 4701                this->value.data.dblptr[elem] = pVals[0].data.dbl;
 4702              } else {
 4703                this->value.undef[elem] = 0;
 4704                this->value.data.dblptr[elem] =
 4705              maxvalue( pVals[0].data.dbl, pVals[1].data.dbl );
 4706              }
 4707           }
 4708            }
 4709         }
 4710         break;
 4711 
 4712         /* Boolean SAO region Functions... scalar or vector dbls */
 4713 
 4714      case near_fct:
 4715         while( row-- ) {
 4716            nelem = this->value.nelem;
 4717            while( nelem-- ) {
 4718           elem--;
 4719           i=3; while( i-- )
 4720              if( vector[i]>1 ) {
 4721             pVals[i].data.dbl =
 4722                theParams[i]->value.data.dblptr[elem];
 4723             pNull[i] = theParams[i]->value.undef[elem];
 4724              } else if( vector[i] ) {
 4725             pVals[i].data.dbl =
 4726                theParams[i]->value.data.dblptr[row];
 4727             pNull[i] = theParams[i]->value.undef[row];
 4728              }
 4729           if( !(this->value.undef[elem] = (pNull[0] || pNull[1] ||
 4730                            pNull[2]) ) )
 4731             this->value.data.logptr[elem] =
 4732               bnear( pVals[0].data.dbl, pVals[1].data.dbl,
 4733                  pVals[2].data.dbl );
 4734            }
 4735         }
 4736         break;
 4737 
 4738      case circle_fct:
 4739         while( row-- ) {
 4740            nelem = this->value.nelem;
 4741            while( nelem-- ) {
 4742           elem--;
 4743           i=5; while( i-- )
 4744              if( vector[i]>1 ) {
 4745             pVals[i].data.dbl =
 4746                theParams[i]->value.data.dblptr[elem];
 4747             pNull[i] = theParams[i]->value.undef[elem];
 4748              } else if( vector[i] ) {
 4749             pVals[i].data.dbl =
 4750                theParams[i]->value.data.dblptr[row];
 4751             pNull[i] = theParams[i]->value.undef[row];
 4752              }
 4753           if( !(this->value.undef[elem] = (pNull[0] || pNull[1] ||
 4754                            pNull[2] || pNull[3] ||
 4755                            pNull[4]) ) )
 4756             this->value.data.logptr[elem] =
 4757              circle( pVals[0].data.dbl, pVals[1].data.dbl,
 4758                  pVals[2].data.dbl, pVals[3].data.dbl,
 4759                  pVals[4].data.dbl );
 4760            }
 4761         }
 4762         break;
 4763 
 4764      case box_fct:
 4765         while( row-- ) {
 4766            nelem = this->value.nelem;
 4767            while( nelem-- ) {
 4768           elem--;
 4769           i=7; while( i-- )
 4770              if( vector[i]>1 ) {
 4771             pVals[i].data.dbl =
 4772                theParams[i]->value.data.dblptr[elem];
 4773             pNull[i] = theParams[i]->value.undef[elem];
 4774              } else if( vector[i] ) {
 4775             pVals[i].data.dbl =
 4776                theParams[i]->value.data.dblptr[row];
 4777             pNull[i] = theParams[i]->value.undef[row];
 4778              }
 4779           if( !(this->value.undef[elem] = (pNull[0] || pNull[1] ||
 4780                            pNull[2] || pNull[3] ||
 4781                            pNull[4] || pNull[5] ||
 4782                            pNull[6] ) ) )
 4783             this->value.data.logptr[elem] =
 4784              saobox( pVals[0].data.dbl, pVals[1].data.dbl,
 4785                  pVals[2].data.dbl, pVals[3].data.dbl,
 4786                  pVals[4].data.dbl, pVals[5].data.dbl,
 4787                  pVals[6].data.dbl );   
 4788            }
 4789         }
 4790         break;
 4791 
 4792      case elps_fct:
 4793         while( row-- ) {
 4794            nelem = this->value.nelem;
 4795            while( nelem-- ) {
 4796           elem--;
 4797           i=7; while( i-- )
 4798              if( vector[i]>1 ) {
 4799             pVals[i].data.dbl =
 4800                theParams[i]->value.data.dblptr[elem];
 4801             pNull[i] = theParams[i]->value.undef[elem];
 4802              } else if( vector[i] ) {
 4803             pVals[i].data.dbl =
 4804                theParams[i]->value.data.dblptr[row];
 4805             pNull[i] = theParams[i]->value.undef[row];
 4806              }
 4807           if( !(this->value.undef[elem] = (pNull[0] || pNull[1] ||
 4808                            pNull[2] || pNull[3] ||
 4809                            pNull[4] || pNull[5] ||
 4810                            pNull[6] ) ) )
 4811             this->value.data.logptr[elem] =
 4812              ellipse( pVals[0].data.dbl, pVals[1].data.dbl,
 4813                   pVals[2].data.dbl, pVals[3].data.dbl,
 4814                   pVals[4].data.dbl, pVals[5].data.dbl,
 4815                   pVals[6].data.dbl );
 4816            }
 4817         }
 4818         break;
 4819 
 4820             /* C Conditional expression:  bool ? expr : expr */
 4821 
 4822          case ifthenelse_fct:
 4823             switch( this->type ) {
 4824             case BOOLEAN:
 4825            while( row-- ) {
 4826           nelem = this->value.nelem;
 4827           while( nelem-- ) {
 4828              elem--;
 4829                      if( vector[2]>1 ) {
 4830                         pVals[2].data.log =
 4831                            theParams[2]->value.data.logptr[elem];
 4832                         pNull[2] = theParams[2]->value.undef[elem];
 4833                      } else if( vector[2] ) {
 4834                         pVals[2].data.log =
 4835                            theParams[2]->value.data.logptr[row];
 4836                         pNull[2] = theParams[2]->value.undef[row];
 4837                      }
 4838              i=2; while( i-- )
 4839             if( vector[i]>1 ) {
 4840                pVals[i].data.log =
 4841                   theParams[i]->value.data.logptr[elem];
 4842                pNull[i] = theParams[i]->value.undef[elem];
 4843             } else if( vector[i] ) {
 4844                pVals[i].data.log =
 4845                   theParams[i]->value.data.logptr[row];
 4846                pNull[i] = theParams[i]->value.undef[row];
 4847             }
 4848              if( !(this->value.undef[elem] = pNull[2]) ) {
 4849                         if( pVals[2].data.log ) {
 4850                            this->value.data.logptr[elem] = pVals[0].data.log;
 4851                            this->value.undef[elem]       = pNull[0];
 4852                         } else {
 4853                            this->value.data.logptr[elem] = pVals[1].data.log;
 4854                            this->value.undef[elem]       = pNull[1];
 4855                         }
 4856                      }
 4857           }
 4858            }
 4859                break;
 4860             case LONG:
 4861            while( row-- ) {
 4862           nelem = this->value.nelem;
 4863           while( nelem-- ) {
 4864              elem--;
 4865                      if( vector[2]>1 ) {
 4866                         pVals[2].data.log =
 4867                            theParams[2]->value.data.logptr[elem];
 4868                         pNull[2] = theParams[2]->value.undef[elem];
 4869                      } else if( vector[2] ) {
 4870                         pVals[2].data.log =
 4871                            theParams[2]->value.data.logptr[row];
 4872                         pNull[2] = theParams[2]->value.undef[row];
 4873                      }
 4874              i=2; while( i-- )
 4875             if( vector[i]>1 ) {
 4876                pVals[i].data.lng =
 4877                   theParams[i]->value.data.lngptr[elem];
 4878                pNull[i] = theParams[i]->value.undef[elem];
 4879             } else if( vector[i] ) {
 4880                pVals[i].data.lng =
 4881                   theParams[i]->value.data.lngptr[row];
 4882                pNull[i] = theParams[i]->value.undef[row];
 4883             }
 4884              if( !(this->value.undef[elem] = pNull[2]) ) {
 4885                         if( pVals[2].data.log ) {
 4886                            this->value.data.lngptr[elem] = pVals[0].data.lng;
 4887                            this->value.undef[elem]       = pNull[0];
 4888                         } else {
 4889                            this->value.data.lngptr[elem] = pVals[1].data.lng;
 4890                            this->value.undef[elem]       = pNull[1];
 4891                         }
 4892                      }
 4893           }
 4894            }
 4895                break;
 4896             case DOUBLE:
 4897            while( row-- ) {
 4898           nelem = this->value.nelem;
 4899           while( nelem-- ) {
 4900              elem--;
 4901                      if( vector[2]>1 ) {
 4902                         pVals[2].data.log =
 4903                            theParams[2]->value.data.logptr[elem];
 4904                         pNull[2] = theParams[2]->value.undef[elem];
 4905                      } else if( vector[2] ) {
 4906                         pVals[2].data.log =
 4907                            theParams[2]->value.data.logptr[row];
 4908                         pNull[2] = theParams[2]->value.undef[row];
 4909                      }
 4910              i=2; while( i-- )
 4911             if( vector[i]>1 ) {
 4912                pVals[i].data.dbl =
 4913                   theParams[i]->value.data.dblptr[elem];
 4914                pNull[i] = theParams[i]->value.undef[elem];
 4915             } else if( vector[i] ) {
 4916                pVals[i].data.dbl =
 4917                   theParams[i]->value.data.dblptr[row];
 4918                pNull[i] = theParams[i]->value.undef[row];
 4919             }
 4920              if( !(this->value.undef[elem] = pNull[2]) ) {
 4921                         if( pVals[2].data.log ) {
 4922                            this->value.data.dblptr[elem] = pVals[0].data.dbl;
 4923                            this->value.undef[elem]       = pNull[0];
 4924                         } else {
 4925                            this->value.data.dblptr[elem] = pVals[1].data.dbl;
 4926                            this->value.undef[elem]       = pNull[1];
 4927                         }
 4928                      }
 4929           }
 4930            }
 4931                break;
 4932             case STRING:
 4933            while( row-- ) {
 4934                   if( vector[2] ) {
 4935                      pVals[2].data.log = theParams[2]->value.data.logptr[row];
 4936                      pNull[2] = theParams[2]->value.undef[row];
 4937                   }
 4938                   i=2; while( i-- )
 4939                      if( vector[i] ) {
 4940                         strcpy( pVals[i].data.str,
 4941                                 theParams[i]->value.data.strptr[row] );
 4942                         pNull[i] = theParams[i]->value.undef[row];
 4943                      }
 4944                   if( !(this->value.undef[row] = pNull[2]) ) {
 4945                      if( pVals[2].data.log ) {
 4946                         strcpy( this->value.data.strptr[row],
 4947                                 pVals[0].data.str );
 4948                         this->value.undef[row]       = pNull[0];
 4949                      } else {
 4950                         strcpy( this->value.data.strptr[row],
 4951                                 pVals[1].data.str );
 4952                         this->value.undef[row]       = pNull[1];
 4953                      }
 4954                   } else {
 4955                      this->value.data.strptr[row][0] = '\0';
 4956                   }
 4957            }
 4958                break;
 4959 
 4960             }
 4961             break;
 4962 
 4963         /* String functions */
 4964             case strmid_fct:
 4965           {
 4966         int strconst = theParams[0]->operation == CONST_OP;
 4967         int posconst = theParams[1]->operation == CONST_OP;
 4968         int lenconst = theParams[2]->operation == CONST_OP;
 4969         int dest_len = this->value.nelem;
 4970         int src_len  = theParams[0]->value.nelem;
 4971 
 4972         while (row--) {
 4973           int pos;
 4974           int len;
 4975           char *str;
 4976           int undef = 0;
 4977 
 4978           if (posconst) {
 4979             pos = theParams[1]->value.data.lng;
 4980           } else {
 4981             pos = theParams[1]->value.data.lngptr[row];
 4982             if (theParams[1]->value.undef[row]) undef = 1;
 4983           }
 4984           if (strconst) {
 4985             str = theParams[0]->value.data.str;
 4986             if (src_len == 0) src_len = strlen(str);
 4987           } else {
 4988             str = theParams[0]->value.data.strptr[row];
 4989             if (theParams[0]->value.undef[row]) undef = 1;
 4990           }
 4991           if (lenconst) {
 4992             len = dest_len;
 4993           } else {
 4994             len = theParams[2]->value.data.lngptr[row];
 4995             if (theParams[2]->value.undef[row]) undef = 1;
 4996           }
 4997           this->value.data.strptr[row][0] = '\0';
 4998           if (pos == 0) undef = 1;
 4999           if (! undef ) {
 5000             if (cstrmid(this->value.data.strptr[row], len,
 5001                 str, src_len, pos) < 0) break;
 5002           }
 5003           this->value.undef[row] = undef;
 5004         }
 5005           }           
 5006           break;
 5007 
 5008         /* String functions */
 5009             case strpos_fct:
 5010           {
 5011         int const1 = theParams[0]->operation == CONST_OP;
 5012         int const2 = theParams[1]->operation == CONST_OP;
 5013 
 5014         while (row--) {
 5015           char *str1, *str2;
 5016           int undef = 0;
 5017