doublereal fname( doublereal * x )    \
   {  float  xx , val ;                  \
      xx     = (float) (*x) ;            \
      val = sfunc( xx , fcode , NULL ) ; \
      return (doublereal) val ;          \
   }

doublereal fname( doublereal * x , doublereal * a ) \
   {  float aux[1] , xx , val ;                        \
      xx     = (float) (*x) ;                          \
      aux[0] = (float) (*a) ;                          \
      val = sfunc( xx , fcode , aux ) ;                \
      return (doublereal) val ;                        \
   }

doublereal fname( doublereal * x, doublereal * a, doublereal * b ) \
   {  float aux[2] , xx , val ;                                       \
      xx     = (float) (*x) ;                                         \
      aux[0] = (float) (*a) ;                                         \
      aux[1] = (float) (*b) ;                                         \
      val = sfunc( xx , fcode , aux ) ;                               \
      return (doublereal) val ;                                       \
   }

doublereal fname( doublereal * x, doublereal * a, doublereal * b,  \
                                                     doublereal * c ) \
   {  float aux[3] , xx , val ;                                       \
      xx     = (float) (*x) ;                                         \
      aux[0] = (float) (*a) ;                                         \
      aux[1] = (float) (*b) ;                                         \
      aux[2] = (float) (*c) ;                                         \
      val = sfunc( xx , fcode , aux ) ;                               \
      return (doublereal) val ;                                       \
   }

00229 { return (doublereal) j0( (double) *x ) ; }

00232 { return (doublereal) j1( (double) *x ) ; }

00235 { return (doublereal) (*x>0) ? y0( (double) *x ) : 0.0 ; }

00238 { return (doublereal) (*x>0) ? y1( (double) *x ) : 0.0 ; }

00241 { return (doublereal) erf( (double) *x ) ; }

00244 { return (doublereal) erfc( (double) *x ) ; }

01180 {
01181     /* System generated locals */
01182     integer ret_val, i__1;
01183 
01184     /* Builtin functions */
01185     integer s_cmp(char *, char *, ftnlen, ftnlen);
01186 
01187     /* Local variables */
01188     static integer ncode;
01189     static char sss[1];
01190 
01191 
01192 
01193 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
01194  */
01195 
01196     /* Parameter adjustments */
01197     c_code__ -= 8;
01198 
01199     /* Function Body */
01200     ret_val = 0;
01201     if (*num_code__ <= 0) {
01202         return ret_val;
01203     }
01204     *(unsigned char *)sss = *(unsigned char *)sym;
01205 
01206     i__1 = *num_code__;
01207     for (ncode = 1; ncode <= i__1; ++ncode) {
01208         if (s_cmp(c_code__ + (ncode << 3), "PUSHSYM", 8L, 7L) == 0) {
01209             if (*(unsigned char *)&c_code__[(ncode + 1) * 8] == *(unsigned 
01210                     char *)sss) {
01211                 ret_val = 1;
01212                 return ret_val;
01213             }
01214         }
01215 /* L1000: */
01216     }
01217 
01218     return ret_val;
01219 } /* hassym_ */

00261 {
00262    double x = *xx ;
00263    int    m = (int)(*mm) ;
00264 
00265    if( m < 0 ) return 1.0 ;    /* bad input */
00266 
00267    switch( m ){
00268     case 0: return 1.0 ;
00269     case 1: return x ;
00270     case 2: return (3.0*x*x-1.0)/2.0 ;
00271     case 3: return (5.0*x*x-3.0)*x/2.0 ;
00272     case 4: return ((35.0*x*x-30.0)*x*x+3.0)/8.0 ;
00273     case 5: return ((63.0*x*x-70.0)*x*x+15.0)*x/8.0 ;
00274     case 6: return (((231.0*x*x-315.0)*x*x+105.0)*x*x-5.0)/16.0 ;
00275     case 7: return (((429.0*x*x-693.0)*x*x+315.0)*x*x-35.0)*x/16.0 ;
00276     case 8: return ((((6435.0*x*x-12012.0)*x*x+6930.0)*x*x-1260.0)*x*x+35.0)/128.0;
00277 
00278            /** 07 Feb 2005: this part generated by Maple, then hand massaged **/
00279 
00280     case 9:
00281       return (0.24609375e1 + (-0.3609375e2 + (0.140765625e3 + (-0.20109375e3
00282               + 0.949609375e2 * x * x) * x * x) * x * x) * x * x) * x;
00283 
00284     case 10:
00285       return -0.24609375e0 + (0.1353515625e2 + (-0.1173046875e3 +
00286               (0.3519140625e3 + (-0.42732421875e3 + 0.18042578125e3 * x * x)
00287              * x * x) * x * x) * x * x) * x * x;
00288 
00289     case 11:
00290       return (-0.270703125e1 + (0.5865234375e2 + (-0.3519140625e3 +
00291              (0.8546484375e3 + (-0.90212890625e3 + 0.34444921875e3 * x * x)
00292              * x * x) * x * x) * x * x) * x * x) * x;
00293 
00294     case 12:
00295       return 0.2255859375e0 + (-0.17595703125e2 + (0.2199462890625e3 +
00296              (-0.99708984375e3 + (0.20297900390625e4 + (-0.1894470703125e4
00297              + 0.6601943359375e3 * x * x) * x * x) * x * x) * x * x) * x * x)
00298              * x * x;
00299 
00300     case 13:
00301       return (0.29326171875e1 + (-0.87978515625e2 + (0.7478173828125e3 +
00302              (-0.270638671875e4 + (0.47361767578125e4 + (-0.3961166015625e4
00303              + 0.12696044921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
00304             * x * x) * x;
00305 
00306     case 14:
00307       return -0.20947265625e0 + (0.2199462890625e2 + (-0.37390869140625e3 +
00308              (0.236808837890625e4 + (-0.710426513671875e4 +
00309              (0.1089320654296875e5 + (-0.825242919921875e4 +
00310             0.244852294921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
00311            * x * x) * x * x;
00312 
00313     case 15:
00314       return (-0.314208984375e1 + (0.12463623046875e3 + (-0.142085302734375e4
00315             + (0.710426513671875e4 + (-0.1815534423828125e5 +
00316               (0.2475728759765625e5 + (-0.1713966064453125e5 +
00317                0.473381103515625e4 * x * x) * x * x) * x * x) * x * x)
00318              * x * x) * x * x) * x * x) * x;
00319 
00320     case 16:
00321       return 0.196380615234375e0 + (-0.26707763671875e2 + (0.5920220947265625e3
00322             + (-0.4972985595703125e4 + (0.2042476226806641e5 +
00323               (-0.4538836059570312e5 + (0.5570389709472656e5 +
00324                (-0.3550358276367188e5 + 0.9171758880615234e4 * x * x) * x * x)
00325             * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
00326 
00327     case 17:
00328       return (0.3338470458984375e1 + (-0.169149169921875e3 +
00329              (0.2486492797851562e4 + (-0.1633980981445312e5 +
00330              (0.5673545074462891e5 + (-0.1114077941894531e6 +
00331              (0.1242625396728516e6 + (-0.7337407104492188e5 +
00332               0.1780400253295898e5 * x * x) * x * x) * x * x) * x * x)
00333            * x * x) * x * x) * x * x) * x * x) * x;
00334 
00335     case 18:
00336       return -0.1854705810546875e0 + (0.3171546936035156e2 +
00337             (-0.8880331420898438e3 + (0.9531555725097656e4 +
00338             (-0.5106190567016602e5 + (0.153185717010498e6 +
00339             (-0.2692355026245117e6 + (0.275152766418457e6 +
00340             (-0.1513340215301514e6 + 0.3461889381408691e5 * x * x) * x * x)
00341            * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
00342 
00343     case 19:
00344       return (-0.3523941040039062e1 + (0.2220082855224609e3 +
00345              (-0.4084952453613281e4 + (0.3404127044677734e5 +
00346              (-0.153185717010498e6 + (0.4038532539367676e6 +
00347              (-0.6420231216430664e6 + (0.6053360861206055e6 +
00348              (-0.3115700443267822e6 + 0.6741574058532715e5 * x * x) * x * x)
00349           * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x;
00350 
00351     case 20:
00352       return 0.1761970520019531e0 + (-0.3700138092041016e2 +
00353             (0.127654764175415e4 + (-0.1702063522338867e5 +
00354             (0.1148892877578735e6 + (-0.4442385793304443e6 +
00355             (0.1043287572669983e7 + (-0.1513340215301514e7 +
00356             (0.1324172688388824e7 + (-0.6404495355606079e6 +
00357              0.1314606941413879e6 * x * x) * x * x) * x * x) * x * x) * x * x)
00358             * x * x) * x * x) * x * x) * x * x) * x * x;
00359    }
00360 
00361         /** if here, m > 20 ==> use recurrence relation **/
00362 
00363    { double k , pk, pkm1, pkm2 ;
00364 
00365      k = 19.0; pkm2 = legendre_( &k , xx ) ;
00366      k = 20.0; pkm1 = legendre_( &k , xx ) ;
00367 
00368      while( k < m ){
00369        k += 1.0 ;
00370        pk = ((2.0*k-1.0)*x*pkm1 - (k-1.0)*pkm2)/k ;
00371        pkm2 = pkm1 ; pkm1 = pk ;
00372      }
00373      return pk ;
00374    }
00375 }

00171 {
00172    PARSER_code * pcode = NULL ;
00173    char sym[4] ;
00174    double atoz[26] ;
00175    int ii , kvar ;
00176 
00177    if( expr == NULL || nt <= 0 || vec == NULL ) return 0 ;  /* bad */
00178 
00179    pcode = PARSER_generate_code( expr ) ;        /* compile */
00180    if( pcode == NULL ) return 0 ;                /* bad news */
00181 
00182    kvar = -1 ;                                   /* find symbol */
00183    for( ii=0 ; ii < 26 ; ii++ ){
00184       sym[0] = 'A' + ii ; sym[1] = '\0' ;
00185       if( PARSER_has_symbol(sym,pcode) ){ kvar = ii ; break ; }
00186    }
00187 
00188    for( ii=0 ; ii < 26 ; ii++ ) atoz[ii] = 0.0 ; /* initialize */
00189 
00190    if( kvar >= 0 ){                              /* the normal case */
00191       for( ii=0 ; ii < nt ; ii++ ){
00192          atoz[kvar] = tz + ii*dt ;
00193          vec[ii]    = PARSER_evaluate_one( pcode , atoz ) ;
00194       }
00195    } else {                                      /* no variable found! */
00196       vec[0] = PARSER_evaluate_one( pcode , atoz ) ;
00197       for( ii=1 ; ii < nt ; ii++ ) vec[ii] = vec[0] ;
00198    }
00199 
00200    free(pcode) ; return 1 ;
00201 }

00066 {
00067    integer num_code ;
00068    double  value ;
00069 
00070    if( pc == NULL || pc->num_code <= 0 ) return 0.0 ;
00071 
00072    num_code = (integer) pc->num_code ;
00073 
00074    value = (double) pareval_( &num_code, pc->c_code,
00075                               (doublereal *) atoz , (ftnlen) 8 ) ;
00076    return value ;
00077 }

00131 {
00132    integer num_code , nvar , ivar , lvec , ldvec ;
00133 
00134    if( pc == NULL || pc->num_code <= 0 ) return ;
00135 
00136    num_code = (integer) pc->num_code ;
00137    lvec     = (integer) nv ;
00138 
00139    parevec_( &num_code , pc->c_code ,
00140              (doublereal *) atoz[0]  , (doublereal *) atoz[1] ,
00141              (doublereal *) atoz[2]  , (doublereal *) atoz[3] ,
00142              (doublereal *) atoz[4]  , (doublereal *) atoz[5] ,
00143              (doublereal *) atoz[6]  , (doublereal *) atoz[7] ,
00144              (doublereal *) atoz[8]  , (doublereal *) atoz[9] ,
00145              (doublereal *) atoz[10] , (doublereal *) atoz[11] ,
00146              (doublereal *) atoz[12] , (doublereal *) atoz[13] ,
00147              (doublereal *) atoz[14] , (doublereal *) atoz[15] ,
00148              (doublereal *) atoz[16] , (doublereal *) atoz[17] ,
00149              (doublereal *) atoz[18] , (doublereal *) atoz[19] ,
00150              (doublereal *) atoz[20] , (doublereal *) atoz[21] ,
00151              (doublereal *) atoz[22] , (doublereal *) atoz[23] ,
00152              (doublereal *) atoz[24] , (doublereal *) atoz[25] ,
00153          &lvec , (doublereal *) vout , (ftnlen) 8 ) ;
00154 
00155    return ;
00156 }

00025 {
00026    logical pr ;
00027    integer num_code ;
00028    int nexp ;
00029    PARSER_code * pc ;
00030    char *exp,cc ; int ii,jj ;  /* 22 Jul 2003 */
00031 
00032    if( expression == NULL ) return NULL ;
00033    nexp = strlen( expression ) ;
00034    if( nexp == 0 ) return NULL ;
00035 
00036    /* 22 Jul 2003: copy into local string, tossing bad stuff */
00037 
00038    exp = AFMALL(char, nexp+4) ;
00039    for( ii=jj=0 ; ii < nexp ; ii++ ){
00040      cc = expression[ii] ;
00041      if( !isspace(cc) && !iscntrl(cc) ) exp[jj++] = cc ;
00042    }
00043    exp[jj] = '\0' ;
00044    nexp = strlen(exp) ; if( nexp == 0 ) return NULL ;
00045 
00046    pc = (PARSER_code *) malloc( sizeof(PARSER_code) ) ;
00047 
00048    pr = (printout) ? TRUE_ : FALSE_ ;
00049 
00050    parser_( exp, &pr, &num_code, pc->c_code, (ftnlen) nexp, (ftnlen) 8 ) ;
00051 
00052    free(exp) ;  /* 22 Jul 2003 */
00053 
00054    if( num_code <= 0 ){ free(pc) ; return NULL ; }
00055 
00056    pc->num_code = (int) num_code ;
00057    return pc ;
00058 }

00088 {
00089    int hh ;
00090    char sss[8] ;
00091    integer num_code ;
00092 
00093    if( !isalpha(sym[0]) ) return 0 ;          /* not alphabetic */
00094 
00095    sss[0] = toupper(sym[0]) ; sss[1] = '\0' ; /* uppercase it */
00096 
00097    num_code = (integer) pc->num_code ;
00098 
00099    hh = (int) hassym_( sss , &num_code , pc->c_code ,
00100                        (ftnlen) 8 , (ftnlen) 8       ) ;
00101 
00102    return hh ;
00103 }

00106 {
00107    int ii ;
00108    static char abet[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZ" ;
00109 
00110    if( pc == NULL || sl == NULL ) return ;
00111 
00112    for( ii=0 ; ii < 26 ; ii++ )
00113       sl[ii] = PARSER_has_symbol( abet+ii , pc ) ;
00114 
00115    return ;
00116 }

00015 { printout = p ; }

00208 {
00209    PARSER_code * pcode = NULL ;
00210    double atoz[26] , val ;
00211    int ii ;
00212 
00213    if( expr == NULL ) return 0 ;                 /* bad */
00214 
00215    pcode = PARSER_generate_code( expr ) ;        /* compile */
00216    if( pcode == NULL ) return 0 ;                /* bad news */
00217 
00218    for( ii=0 ; ii < 26 ; ii++ ) atoz[ii] = 0.0 ; /* initialize */
00219 
00220    val = PARSER_evaluate_one( pcode , atoz ) ;
00221 
00222    free(pcode) ; return val ;
00223 }

00069 {
00070    float stat = -1.0 ;   /* error flag */
00071 
00072    if( pval >= 0.999999 ) return 0.0 ;  /* WTF */
00073 
00074    if( stataux == NULL && statcode != FUNC_ZT_TYPE ) return pval ;
00075 
00076    switch( statcode ){  /* if statcode is illegal, will return -1 */
00077 
00078       /** the routines below are in mri_stats.c **/
00079 
00080       case FUNC_COR_TYPE:
00081          stat = correl_p2t( pval , stataux[0] , stataux[1] , stataux[2] ) ;
00082       break ;
00083 
00084       case FUNC_TT_TYPE:
00085          stat = student_p2t( pval , stataux[0] ) ;
00086       break ;
00087 
00088       case FUNC_FT_TYPE:
00089          stat = fstat_p2t( pval , stataux[0] , stataux[1] ) ;
00090       break ;
00091 
00092       case FUNC_ZT_TYPE:                /* only type that doesn't */
00093          stat = normal_p2t( pval ) ;    /* use stataux parameters */
00094       break ;
00095 
00096       case FUNC_CT_TYPE:
00097          stat = chisq_p2t( pval , stataux[0] ) ;
00098       break ;
00099 
00100       case FUNC_BT_TYPE:
00101          stat = beta_p2t( pval , stataux[0] , stataux[1] ) ;
00102       break ;
00103 
00104       case FUNC_BN_TYPE:
00105          stat = binomial_p2t( pval , stataux[0] , stataux[1] ) ;
00106       break ;
00107 
00108       case FUNC_GT_TYPE:
00109          stat = gamma_p2t( pval , stataux[0] , stataux[1] ) ;
00110       break ;
00111 
00112       case FUNC_PT_TYPE:
00113          stat = poisson_p2t( pval , stataux[0] ) ;
00114       break ;
00115    }
00116 
00117    return stat ;
00118 }

00017 {
00018    float pval = -1.0 ;   /* error flag */
00019 
00020    if( stataux == NULL && statcode != FUNC_ZT_TYPE ) return pval ;
00021 
00022    if( thr == 0.0 ) return 1.0 ;
00023 
00024    switch( statcode ){  /* if statcode is illegal, will return -1 */
00025 
00026       case FUNC_COR_TYPE:
00027          pval = correl_t2p( thr , stataux[0] , stataux[1] , stataux[2] ) ;
00028       break ;
00029 
00030       case FUNC_TT_TYPE:
00031          pval = student_t2p( thr , stataux[0] ) ;
00032       break ;
00033 
00034       case FUNC_FT_TYPE:
00035          pval = fstat_t2p( thr , stataux[0] , stataux[1] ) ;
00036       break ;
00037 
00038       case FUNC_ZT_TYPE:               /* only type that doesn't */
00039          pval = normal_t2p( thr ) ;    /* use stataux parameters */
00040       break ;
00041 
00042       case FUNC_CT_TYPE:
00043          pval = chisq_t2p( thr , stataux[0] ) ;
00044       break ;
00045 
00046       case FUNC_BT_TYPE:
00047          pval = beta_t2p( thr , stataux[0] , stataux[1] ) ;
00048       break ;
00049 
00050       case FUNC_BN_TYPE:
00051          pval = binomial_t2p( thr , stataux[0] , stataux[1] ) ;
00052       break ;
00053 
00054       case FUNC_GT_TYPE:
00055          pval = gamma_t2p( thr , stataux[0] , stataux[1] ) ;
00056       break ;
00057 
00058       case FUNC_PT_TYPE:
00059          pval = poisson_t2p( thr , stataux[0] ) ;
00060       break ;
00061    }
00062 
00063    return pval ;
00064 }

00123 {
00124    float zscore = thr ;
00125 
00126    if( stataux == NULL && statcode != FUNC_ZT_TYPE ) return zscore ;
00127 
00128    switch( statcode ){  /* if statcode is illegal, will return -1 */
00129 
00130       /** the routines below are in mri_stats.c **/
00131 
00132       case FUNC_COR_TYPE:
00133          zscore = correl_t2z( thr , stataux[0] , stataux[1] , stataux[2] ) ;
00134       break ;
00135 
00136       case FUNC_TT_TYPE:
00137          zscore = student_t2z( thr , stataux[0] ) ;
00138       break ;
00139 
00140       case FUNC_FT_TYPE:
00141          zscore = fstat_t2z( thr , stataux[0] , stataux[1] ) ;
00142       break ;
00143 
00144       case FUNC_ZT_TYPE:                 /* only type that doesn't */
00145          zscore = normal_t2z( thr ) ;    /* use stataux parameters */
00146       break ;
00147 
00148       case FUNC_CT_TYPE:
00149          zscore = chisq_t2z( thr , stataux[0] ) ;
00150       break ;
00151 
00152       case FUNC_BT_TYPE:
00153          zscore = beta_t2z( thr , stataux[0] , stataux[1] ) ;
00154       break ;
00155 
00156       case FUNC_BN_TYPE:
00157          zscore = binomial_t2z( thr , stataux[0] , stataux[1] ) ;
00158       break ;
00159 
00160       case FUNC_GT_TYPE:
00161          zscore = gamma_t2z( thr , stataux[0] , stataux[1] ) ;
00162       break ;
00163 
00164       case FUNC_PT_TYPE:
00165          zscore = poisson_t2z( thr , stataux[0] ) ;
00166       break ;
00167    }
00168 
00169    return zscore ;
00170 }

00249 {
00250    static first=1 ;
00251    doublereal val ;
00252    if( first ){ srand48((long)time(NULL)); first=0; }
00253    val = (doublereal) drand48() ;
00254    return val ;
00255 }