home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / formu220.zip / formulc.c < prev    next >
C/C++ Source or Header  |  1995-05-17  |  35KB  |  1,341 lines
  1. /*FORMULC.C 2.2           as of 12/13/94        */
  2. /* definitive  version */
  3. /*A fast interpreter of mathematical functions */
  4.  
  5. /*Copyright (c) 1995 by Harald Helfgott        */
  6. /* This program must be distributed with its corresponding README.DOC */
  7. /* The full copyright and availability notice is in README.DOC          */
  8. /*     This program is provided "as is", without any explicit or */
  9. /* implicit warranty. */
  10.  
  11.  
  12. /* Programmer's Address:
  13.         Harald Helfgott
  14.         MB 1807, Brandeis University
  15.         P.O. Box 9110
  16.         Waltham, MA 02254-9110
  17.         U.S.A.
  18.         hhelf@cs.brandeis.edu
  19.            OR
  20.          (during the summer)
  21.         2606 Willett Apt. 427
  22.         Laramie, Wyoming 82070
  23.         seiere@uwyo.edu */
  24. /* Differences between versions 2.1 and 2.2:
  25.        1) FORMULC now runs in machines which don't like unaligned doubles.
  26.        It was necessary to encapsulate the coded-function type;
  27.        now, you must use formu instead of UCHAR.
  28.        2) FORMULC 2.2 has a main function which can be activated
  29.        by defining STAND_ALONE.
  30.        3) FORMULC 2.2 gives error messages.
  31.            
  32.        Modification 2) is due to Ralf Grosse Kunstleve;
  33.        Modification 1) is due to Mr. Grosse and the author.
  34.        Mr. Grosse's address: ralf@kristall.erdw.ethz.ch
  35.        Thanks, Ralf!                                               
  36.  
  37.    Differences between the September 1994 and the December 1994 versions:
  38.      read CHANGES.DOC                                     */
  39.  
  40. #include <math.h>
  41. #include <string.h>
  42. #include <stdlib.h>
  43. #include <stdio.h>
  44. #include <stdarg.h>
  45. #include <errno.h>
  46. #include <ctype.h>
  47. #include <time.h>
  48. #include "formulc.h"
  49.  
  50. #undef FORMU_DEBUG_ON
  51. /* Substitute #define for #undef to trace evaluation */
  52. #ifdef FORMU_DEBUG_ON
  53.  #define DBG(anything)  anything
  54. #else
  55.  #define DBG(anything)
  56.   /* nothing */
  57. #endif
  58.  
  59. static double pi(void);
  60.  
  61. static double value(formu function);
  62.  
  63. static const char *i_error; /*pointer to the character in source[]
  64.             that causes an error */
  65. #define Max_ctable 255
  66.    /*maximum number of items in a table of constants */
  67.    /* Max_ctable must be less than 256 */
  68. static int i_pctable; /* number of items in a table of constants -
  69.                          used only by the translating functions */
  70. static double *i_ctable; /*current table of constants -
  71.                            used only by the translating functions */
  72. static UCHAR *i_trans(UCHAR *function, char *begin, char *end);
  73. static char  *my_strtok(char *s);
  74. static UCHAR *comp_time(UCHAR *function, UCHAR *fend, int npars);
  75.  
  76. static char *errmes = NULL;
  77. static void fset_error(char *);
  78.  
  79. static double param['z'-'a'+1];
  80. typedef struct {
  81.   char *name;
  82.   Func f;    /* pointer to function*/
  83.   int n_pars; /* number of parameters (0, 1, 2 or 3) */
  84.   int varying; /* Does the result of the function vary
  85.           even when the parameters stay the same?
  86.           varying=1 for e.g. random-number generators. */ 
  87. } formu_item;
  88. #define TABLESIZE 256
  89. #define STD_LIB_NUM 13
  90. static formu_item ftable[TABLESIZE]=
  91. {
  92.   {"exp", exp,1,0},
  93.   {"ln",  log,1,0},
  94.   {"sin", sin,1,0},
  95.   {"cos", cos,1,0},
  96.   {"tan", tan,1,0},
  97.   {"asin", asin,1,0},
  98.   {"acos", acos,1,0},
  99.   {"atan", atan,1,0},
  100.   {"atan2",(Func) atan2,2,0},
  101.   {"abs",  fabs,1,0},
  102.   {"sqrt",  sqrt,1,0},
  103.   {"pi", (Func) pi,0,0},
  104.   {"rnd", (Func) rnd, 0, 1}, /*returns a random number from 0 to 1 */
  105.   {NULL,NULL,0}};
  106. /* please, don't use this array directly;
  107.    methods are provided for its manipulation */
  108.  
  109. /*********************************************************/
  110. /* The following routines manipulate the table of functions */
  111.  
  112. int read_table(int i, char *name, int *n_pars, int *varying)
  113. /* returns 1 if succesful */
  114. /* returns 0 otherwise */
  115. {
  116.  if(!ftable[i].f) {
  117.   fset_error("index out of bounds"); 
  118.   return 0;
  119.  }
  120.  else {
  121.   strcpy(name,ftable[i].name);
  122.   *n_pars = ftable[i].n_pars;
  123.   *varying = ftable[i].varying;
  124.   fset_error(NULL);
  125.   return 1;
  126.  }
  127. }
  128.  
  129. int where_table(char *name)
  130. /* If the function exists, where_table() returns the index of its name
  131.     in the table. Otherwise, it returns -1. */
  132. {
  133.  formu_item *table_p;
  134.  
  135.  for(table_p=ftable; table_p->f != NULL &&
  136.     strcmp(name,table_p->name); table_p++)
  137.    ;
  138.  if(table_p->f == NULL) /*The end of the table has been reached,
  139.          but name[] is not there. */
  140.   {
  141.     fset_error("function not found");
  142.     return -1;
  143.   }
  144.  else {
  145.    fset_error(NULL);
  146.    return table_p - ftable;
  147.  }
  148. }
  149.  
  150. int fdel(char *name)
  151. /* If the function exists, it is deleted and a non-negative value
  152.     is returned. */
  153. /* Otherwise, -1 is returned. */
  154. /* Original library functions may not be deleted. */
  155. {
  156.  int place;
  157.  formu_item *scan;
  158.  
  159.  if((place=where_table(name)) == -1) 
  160.   return -1; /* there is an error message already */
  161.  if(place<STD_LIB_NUM) {
  162.   fset_error("original functions may not be deleted");
  163.   return -1;
  164.  }
  165.  free(ftable[place].name);
  166.  for(scan = &ftable[place]; scan->f!=NULL; scan++) {
  167.   DBG(printf("%s \t",scan->name));
  168.   scan->name  =  (scan+1)->name;
  169.   scan->f     =  (scan+1) -> f;
  170.   scan->n_pars = (scan+1) -> n_pars;
  171.  }
  172.  fset_error(NULL);
  173.  return scan-ftable;
  174. } /*end of fdel */
  175.  
  176. int fnew(char *name, Func f, int n_pars, int varying)
  177. /* 0 is rendered if there is an error */
  178. /* 1 is rendered otherwise */
  179. {
  180.  formu_item *where;
  181.  
  182.  if(n_pars<0 || n_pars>3) {
  183.   fset_error("invalid number of parameters"); 
  184.   return 0;
  185.  }
  186.  for(where=ftable; where->f != NULL && strcmp(name,where->name); where++);
  187.  if(where->f != NULL) {
  188.   where->f=f;
  189.   where->varying = varying;
  190.   where->n_pars = n_pars;   /*old function is superseded */
  191.   fset_error(NULL);
  192.   return 1;
  193.  } else if((where-ftable) >= TABLESIZE-1) {
  194.   fset_error("function table full");
  195.   return 0; 
  196.  }
  197.  else {
  198.   where->name = (char *) calloc(strlen(name)+1, sizeof(char));
  199.   if(where->name==NULL) {
  200.     fset_error("no memory");
  201.     return 0;
  202.   }
  203.   strcpy(where->name,name);
  204.   where->f=f;
  205.   where->varying = varying;
  206.   where->n_pars = n_pars;
  207.   fset_error(NULL);
  208.   return 1;
  209.  }
  210. }  /* end of fnew */
  211.  
  212. /***********************************************************/
  213. /* Error functions                                         */
  214.  
  215. static void fset_error(char *s)
  216. /* fset_error(NULL) and fset_error("") erase 
  217.    any previous error message */
  218. {
  219.  if (s == NULL || *s == '\0') errmes = NULL; /* an empty error message means
  220.                                    that there is no error */
  221.  else errmes = s;
  222. }
  223.  
  224. const char *fget_error(void)
  225. {
  226.  return errmes;
  227. }
  228.  
  229. /**********************************************************/
  230. /* Evaluating functions                                   */
  231.  
  232. double fval_at(formu function)
  233. {
  234.   fset_error(NULL);
  235.   return value(function);
  236. }
  237.  
  238. void make_var(char var, double value)
  239. /*for use with fval_at */
  240. /* make_var('x',3); makes x=3 */
  241. {
  242.   param[var-'a']=value;
  243. }
  244.  
  245. double f_x_val(formu function, double x)
  246. {
  247.  fset_error(NULL);
  248.  param['x'-'a']=x;
  249.  return value(function);
  250. }
  251.  
  252. double fval(formu function, char *args, ...)
  253. {
  254.  va_list ap;
  255.  double result;
  256.  
  257.  DBG(puts("Enter fval"));
  258.  fset_error(NULL);
  259.  va_start(ap, args);
  260.  while(*args)
  261.   param[(*args++)-'a'] = va_arg(ap, double);
  262.  va_end(ap);
  263.  result=value(function);
  264.  return result;
  265. }
  266.  
  267.  
  268. #define BUFSIZE 500
  269. /* bufsize is the size of the stack for double value(formu func) */
  270. static double value(formu func)
  271. {
  272.  double buffer[BUFSIZE];
  273.  register double *bufp = buffer;
  274.       /* points to the first free space in the buffer */
  275.  double x,y,z;
  276.  register double result;
  277.  register UCHAR *function=func.code;
  278.  register double *ctable=func.ctable;
  279.  
  280.  DBG(puts("Entering value"));
  281.  if(!function) {
  282.    fset_error("empty coded function");
  283.    return 0; /* non-existent function; result of
  284.                 an unsuccesful call to translate */
  285.  }
  286.  for(;;) {
  287.    switch(*function++) {
  288.     case '\0':goto finish; /* there is a reason for this "goto":
  289.                   this function must be as fast as possible */
  290.     case 'D': *bufp++ = ctable[*function++];
  291.           DBG(printf("%g ",ctable[*(function-1)]));
  292.           break;
  293.     case 'V': *bufp++ = param[(*function++)-'a'];
  294.           DBG( printf("(%c = %g)   ",*(function-1),*(bufp-1)) );
  295.           break;
  296.     case 'M':DBG(printf("(Unary -) "));
  297.          result = -(*--bufp);
  298.          *bufp++ = result;
  299.          break;
  300.     case '+':DBG(printf("+ "));
  301.          y = *(--bufp);
  302.          result = y + *(--bufp);
  303.          *bufp++ = result;
  304.       break;
  305.     case '-':DBG(printf("- "));
  306.          y = *--bufp;
  307.          result= *(--bufp) - y;
  308.          *bufp++ = result;
  309.          break;
  310.     case '*':DBG(printf("* "));
  311.          y = *(--bufp);
  312.          result = *(--bufp) * y;
  313.          *bufp++ = result;
  314.          break;
  315.     case '/':DBG(printf("/ "));
  316.          y = *--bufp;
  317.          result = *(--bufp) / y;
  318.          *bufp++ = result;
  319.          break;
  320.     case '^':DBG(printf("^ "));
  321.          y = *--bufp;
  322.          result = pow(*(--bufp),y);
  323.          *bufp++ = result;
  324.          break;
  325.     case 'F':DBG(printf("%s ",ftable[*function].name));
  326.          switch(ftable[*function].n_pars) {
  327.            case 0:*bufp++ = ((Func0)ftable[*function++].f)();
  328.               break;
  329.            case 1:x = *--bufp;
  330.               *bufp++ = ftable[*function++].f(x);
  331.               break;
  332.         case 2:y = *--bufp;
  333.            x = *--bufp;
  334.            *bufp++ = ((Func2)ftable[*function++].f)(x,y);
  335.               break;
  336.            case 3:z = *--bufp;
  337.               y = *--bufp;
  338.               x = *--bufp;
  339.               *bufp++ = ((Func3)ftable[*function++].f)(x,y,z);
  340.               break;
  341.            default:fset_error("I2: too many parameters\n");
  342.                return 0;
  343.           }
  344.          break;
  345.     default:fset_error("I1: unrecognizable operator");
  346.      return 0;
  347.    }
  348.  }
  349.  finish: if((bufp-buffer)!=1)
  350.       {
  351.         DBG(putchar('\n'));
  352.         fset_error("I3: corrupted buffer");
  353.         DBG(printf("Buffer: "));
  354.         DBG(while(bufp-- > buffer))
  355.           DBG(printf("%g ",*bufp));
  356.        DBG(putchar('\n'));
  357.      }
  358.  DBG(else putchar('\n'));
  359.  return buffer[0];
  360. } /* end of value */
  361.  
  362. /**********************************************************/
  363. /* Manipulation of data of type formu                     */
  364.  
  365. void destrf(formu old)
  366. {
  367.  fset_error(NULL);
  368.  free(old.code);
  369.  free(old.ctable);
  370. }
  371.  
  372. void make_empty(formu f)
  373. {
  374.  fset_error(NULL);
  375.  f.code=NULL;
  376.  f.ctable=NULL;
  377. }
  378.  
  379. int fnot_empty(formu f)
  380. {
  381.  fset_error(NULL);
  382.  return(f.code!=NULL);
  383. }
  384.  
  385. /*********************************************************/
  386. /* Interpreting functions                                */
  387.  
  388. static int isoper(char c)
  389. {
  390.  return ((c == '+') || (c == '-') || (c == '*') || (c == '/')
  391.             || (c == '^'));
  392. }
  393.  
  394. static int is_code_oper(UCHAR c)
  395. {
  396.  return ((c == '+') || (c == '-') || (c == '*') || (c == '/')
  397.             || (c == '^') || (c == 'M'));
  398. }
  399. static int isin_real(char c)
  400. /* + and - are not included */
  401. {
  402.  return (isdigit(c) || c=='.' || c=='E');
  403. }
  404.  
  405. size_t max_size(const char *source)
  406. /* gives an upper estimate of the size required for
  407.    the coded form of source (including the final '\0') */
  408. /* Take care when modifying: the upper estimate
  409.    returned by max_size must not also accomodate
  410.    *proper* output, but also *improper* output
  411.    which takes place before the translator detects an error. */
  412. {
  413.  int numbers=0;
  414.  int functions=0;
  415.  int operators=0;
  416.  int variables=0;
  417.  
  418. /* const size_t func_size=2*sizeof(UCHAR); */ /* not needed */
  419.  const size_t var_size=2*sizeof(UCHAR);
  420.  const size_t num_size=sizeof(UCHAR)+sizeof(double);
  421.  const size_t op_size=sizeof(UCHAR);
  422.  const size_t end_size=sizeof('\0');
  423.  
  424.  const char *scan;
  425.  
  426.  for(scan=source; *scan; scan++)
  427.   if(isalpha(*scan) && (*scan != 'E'))
  428.   {
  429.     if(isalpha(*(scan+1))) ; /* it is a function name,
  430.                 it will be counted later on */
  431.     else 
  432.      if(*(scan+1) == '(')  functions++;
  433.      else variables++;
  434.   }
  435.  
  436.  if(isoper(*source)) operators++;
  437.  if(*source != '\0')
  438.   for(scan = source+1; *scan; scan++)
  439.    if(isoper(*scan) && *(scan-1) != 'E') operators++;
  440.  
  441.  /* counting numbers.. */
  442.  scan=source;
  443.  while(*scan)
  444.   if(isin_real(*scan) || ((*scan == '+' || *scan == '-') &&
  445.                scan>source && *(scan-1)=='E'))
  446.    {numbers++;
  447.     scan++;
  448.     while(isin_real(*scan) || ((*scan == '+' || *scan == '-') &&
  449.                 scan>source && *(scan-1)=='E'))
  450.      scan++;
  451.    }
  452.   else scan++;
  453.  
  454.  return(numbers*num_size + operators*op_size + functions*num_size
  455.              + variables*var_size + end_size);
  456.  /*Do you wonder why "function" is multiplied with "num_size"
  457.    and not with func_size? This function calculates an upper-bound
  458.    (i.e. pessimistic) estimate. It supposes that all functions are
  459.    converted into doubles by comp_time. For example, pi() actually
  460.    becomes a double. */
  461. }
  462.  
  463. /***********************************************************/
  464. /* Interface for interpreting functions                     */
  465.  
  466. formu translate(const char *sourc, const char *args, int *leng,          
  467.                                                         int *error)
  468. {
  469.  UCHAR *result;
  470.  char *source;
  471.  const char *scan, *scarg;
  472.  UCHAR *function;
  473.  UCHAR *nfunc; /* used to free unused heap space */
  474.  size_t size_estim; /* upper bound for the size of the
  475.                     coded function */
  476.  double *ctable;
  477.  formu returned; /*the value to be returned by the function
  478.            is stored here */
  479.  
  480.  
  481.  i_error=NULL;
  482.  
  483.  source = (char *) malloc(strlen(sourc) + 1);
  484.  if(source==NULL) {
  485.    fset_error("no memory");
  486.    *leng = 0;
  487.    *error = 0; /* first character */
  488.    returned.code = NULL;
  489.    returned.ctable = NULL;
  490.    return(returned);
  491.  }
  492.  strcpy(source,sourc);
  493.  /* FORMULC's routines must have their own copy of sourc
  494.     because the copy could be temporarily modified.
  495.     Modifying a string constant can make some computers crash. */
  496.  
  497.  /* search for undeclared parameters */
  498.  for(scan=source; *scan != '\0'; scan++) {
  499.   if(islower(*scan) && !isalpha(*(scan+1)) &&
  500.       (scan==source || !isalpha(*(scan-1))) ) {
  501.    for(scarg=args; *scarg != '\0' && *scarg != *scan; scarg++)
  502.      ;
  503.    if(*scarg == '\0') /*parameter not found */
  504.     {
  505.      i_error = scan;
  506.  
  507.      fset_error("undeclared parameter");
  508.      *leng = 0;
  509.      *error = i_error - source;
  510.      returned.code=NULL;
  511.      returned.ctable=NULL;
  512.      free(source);
  513.      return(returned);
  514.     }
  515.   }
  516.  }  /* end of search for undeclared... */
  517.  
  518.  size_estim=max_size(source); /* upper estimate of the size
  519.                  of the coded function,
  520.                  which doesn't exist yet */
  521.  
  522.  if(!(function = (UCHAR *) malloc(size_estim))) {
  523.   /* out of memory */
  524.   fset_error("no memory");
  525.   *leng = 0;
  526.   *error = -1;
  527.   returned.code=NULL;
  528.   returned.ctable=NULL;
  529.   free(source);
  530.   return (returned);
  531.  }
  532.  
  533.  /*table of memory is cleaned: */
  534.  i_pctable=0;   
  535.  if(!(i_ctable = (double *) malloc(Max_ctable * sizeof(double)) )) {
  536.   /* out of memory */
  537.   fset_error("no memory"); 
  538.   free(function);
  539.   *leng = 0;
  540.   *error = -1;
  541.   returned.code=NULL;
  542.   returned.ctable=NULL;
  543.   free(source);
  544.   return (returned);
  545.  }
  546.  ctable = i_ctable;
  547.  
  548.  fset_error(NULL);
  549.  /* THIS IS THE CORE STATEMENT */
  550.  result=i_trans(function,(char *) source,(char *) source+strlen(source));
  551.  
  552.  if(!result || fget_error()) { 
  553.   free(function);
  554.   free(i_ctable);
  555.   *leng = 0;
  556.   if(i_error)
  557.    *error = i_error-source;
  558.   else *error = -1; /* internal error or out of memory */
  559.   returned.code=NULL;
  560.   returned.ctable=NULL;
  561.   free(source);
  562.   return (returned);
  563.  }
  564.  else { /* OK */
  565.   *result = '\0';
  566.   *error = -1;
  567.   *leng = result-function;
  568.  
  569.   /* free unused heap space.. */
  570.   if(((*leng)+1) * sizeof(UCHAR) > size_estim)
  571.    /* one must use (*leng)+1 instead of *leng because '\0'
  572.       has not been counted */
  573.    { 
  574.     fset_error("I4: size estimate too small");
  575.     returned.code=NULL;
  576.     returned.ctable=NULL;
  577.     free(source);
  578.     return (returned);
  579.    }
  580.   else if(((*leng)+1) * sizeof(UCHAR) < size_estim) {
  581.     nfunc = (UCHAR *) malloc(((*leng)+1) * sizeof(UCHAR));
  582.       if(nfunc) {
  583.     memcpy( nfunc, function, ((*leng)+1) * sizeof(UCHAR) );
  584.     free(function);
  585.     function=nfunc;
  586.       }
  587.   } /* end of if-else stairs */
  588.  
  589.   /* free heap space hoarded by i_ctable.. */
  590.   if(i_pctable<Max_ctable) {
  591.     ctable = (double *) malloc(i_pctable * sizeof(double));
  592.     if(ctable) {
  593.       memcpy(ctable, i_ctable, i_pctable * sizeof(double));
  594.       free(i_ctable);
  595.     } else ctable = i_ctable;
  596.   } else ctable = i_ctable;
  597.  
  598.   returned.code=function;
  599.   returned.ctable=ctable;
  600.   fset_error(NULL);
  601.   free(source);
  602.   return(returned);
  603.  } /* end of OK */
  604. }  /* end of translate */
  605.  
  606. static UCHAR *comp_time(UCHAR *function, UCHAR *fend, int npars)
  607.   /* calculates at "compile time" */
  608.   /* Postconditions: If the coded expression in *function..*(fend-1)
  609.       can be calculated, its value is stored in *function..*(fend-1) */
  610.   /* comp_time returns a pointer to the first character after the
  611.      end of the coded function; if this function cannot be evaluated
  612.      at compile time, comp_time returns fend, of course.  */
  613.   /* Only memory positions from *function to *comp_time are touched. */
  614. {
  615.   UCHAR *scan;
  616.   UCHAR temp;
  617.   double tempd;
  618.   int i;
  619.   formu trans;
  620.  
  621.   DBG(puts("Entering comp_time"));
  622.   scan=function;
  623.   for(i=0; i<npars; i++) {
  624.    if(*scan++ != 'D') return fend;
  625.    scan++;
  626.   }
  627.  
  628.   if(!( ( scan == fend - (sizeof((UCHAR) 'F')+sizeof(UCHAR))
  629.        && *(fend-2) == 'F' && ftable[*(fend-1)].varying == 0) ||
  630.      ( scan == fend - sizeof(UCHAR)
  631.        && is_code_oper(*(fend-1)) ) )
  632.     )
  633.     /* compile-time evaluation is done only if 
  634.        1) everything but the ending function consists of doubles
  635.        AND
  636.        2) the function does not vary when its parameters remain the same
  637.           (i.e. random-number generators are not evaluated at compile time)
  638.       */
  639.    return fend;
  640.  
  641.   temp = *fend;
  642.   *fend = '\0';
  643.    
  644.   trans.code=function;
  645.   trans.ctable=i_ctable;
  646.   tempd = value(trans);
  647.   *fend = temp;
  648.   *function++ = 'D';
  649.   i_pctable -= npars;
  650.   *function++ = (UCHAR) i_pctable;
  651.   i_ctable[i_pctable++] = tempd;
  652.   
  653.   DBG(puts("Exiting comp_time succesfully"));
  654.   return function;
  655. } /* end of comp_time */
  656.  
  657. static char *my_strtok(char *s)
  658. /* a version of strtok that respects parentheses */
  659. /* token delimiter = comma */
  660. {
  661.  int pars;
  662.  static char *token=NULL;
  663.  char *next_token;
  664.  
  665.  if(s!=NULL) token=s;
  666.  else if(token!=NULL) s=token;
  667.  else return NULL;
  668.  
  669.  for(pars=0; *s != '\0' && (*s != ',' || pars!=0); s++) {
  670.    if(*s == '(') ++pars;
  671.    if(*s == ')') --pars;
  672.  }
  673.  if(*s=='\0') {
  674.   next_token=NULL;
  675.   s=token;
  676.  
  677.   token=next_token;
  678.   DBG(printf("The token is: %s\n",s));
  679.   return s;
  680.  } else {
  681.   *s = '\0';
  682.   next_token=s+1;
  683.   s=token;
  684.  
  685.   token=next_token;
  686.   DBG(printf("The token is: %s\n",s));
  687.   return s;
  688.  }
  689. } /* end of my_strtok */
  690.  
  691.  
  692. /************************************************************/
  693. /* Here begins the core of interpretation                   */
  694.  
  695. #define TWO_OP {                                 \
  696.     if((tempu=i_trans(function,begin,scan)) &&      \
  697.        (temp3=i_trans(tempu,scan+1,end)) ) {       \
  698.     *temp3++ = *scan; /* copies operator */                 \
  699.     temp3 = comp_time(function,temp3,2); /*tries to simplify expression*/ \
  700.    if(fget_error()) return NULL; /* internal error in comp_time */  \
  701.    else return temp3; /* expression has been translated */ \
  702.   } else return NULL; /* something is wrong with the operands */ \
  703.  }
  704.  
  705. #define ERROR_MEM {    \
  706.    fset_error("no memory"); \
  707.    i_error=NULL;  \
  708.    return NULL;   \
  709.   }
  710. static UCHAR *i_trans(UCHAR *function, char *begin, char *end)
  711.  /* the source is *begin .. *(end-1) */
  712.  /* returns NULL if a normal error or an internal error occured; 
  713.     otherwise, returns a pointer 
  714.     to the first character after the end of function[] */
  715.  /* i_trans() does not write a '\0' at the end of function[], */
  716.  /* but it MAY touch its end (i.e. *i_trans) without changing it.*/
  717. {
  718.  int pars;     /* parentheses */
  719.  char *scan;
  720.  UCHAR *tempu, *temp3;
  721.  char *temps;
  722.  char tempch;
  723.  double tempd;
  724.  char *endf;     /* points to the opening
  725.             parenthesis of a function (e.g. of sin(x) ) */
  726.  int n_function;
  727.  int space;
  728.  int i;
  729.  
  730.  char *paramstr[MAXPAR];
  731.  char *par_buf;
  732.  
  733.  if(begin>=end) {
  734.   fset_error("missing operand"); 
  735.   i_error = begin;
  736.   return NULL;
  737.  }
  738.  
  739.  DBG(tempch = *end);
  740.  DBG(*end = '\0');
  741.  DBG(puts(begin));
  742.  DBG(*end = tempch);
  743.  
  744.  /* test paired parentheses */
  745.  for(pars=0, scan=begin; scan<end && pars>=0; scan++) {
  746.   if(*scan == '(') pars++;
  747.   else if(*scan == ')') pars--;
  748.  }
  749.  if(pars<0 || pars>0) {
  750.   fset_error("unmatched parentheses"); 
  751.   i_error = scan-1;
  752.   return NULL;
  753.  }
  754.  
  755.  /* plus and binary minus */
  756.  for(pars=0, scan=end-1; scan>=begin; scan--) {
  757.   if(*scan == '(') pars++;
  758.   else if(*scan == ')') pars--;
  759.   else if(!pars && (*scan == '+' || ((*scan == '-') && scan!=begin))
  760.                       /* recognizes unary
  761.                          minuses */
  762.          && (scan==begin || *(scan-1)!='E') )
  763.       /* be wary of misunderstanding exponential notation */
  764.    break;
  765.  }
  766.  
  767.  if(scan >= begin) TWO_OP
  768.  
  769.  /* multiply and divide */
  770.  for(pars=0, scan=end-1; scan>=begin; scan--) {
  771.   if(*scan == '(') pars++;
  772.   else if(*scan == ')') pars--;
  773.   else if(!pars && (*scan == '*' || *scan == '/' ))
  774.    break;
  775.  }
  776.  
  777.  if(scan >= begin) TWO_OP
  778.  
  779.  /* unary minus */
  780.  if(*begin == '-') {
  781.    tempu=i_trans(function,begin+1,end);
  782.    if(tempu) {
  783.      *tempu++ = 'M';
  784.      tempu=comp_time(function,tempu,1); /*tries to simplify
  785.                       expression*/
  786.      if(fget_error()) return NULL; /* internal error in comp_time */
  787.      else return tempu;
  788.    } else return NULL;
  789.  }
  790.  
  791.  /* power */
  792.  for(pars=0, scan=end-1; scan>=begin; scan--) {
  793.   if(*scan == '(') pars++;
  794.   else if(*scan == ')') pars--;
  795.   else if(!pars && (*scan == '^'))
  796.    break;
  797.  }
  798.  
  799.  if(scan >= begin) TWO_OP
  800.  
  801.  /* erase white space */
  802.  while(isspace(*begin))
  803.   begin++;
  804.  while(isspace(*(end-1)))
  805.   end--;
  806.  
  807.  /* parentheses around the expression */
  808.  if(*begin == '(' && *(end-1) == ')')
  809.   return i_trans(function,begin+1,end-1);
  810.  
  811.  /* variable */
  812.  if(end == begin+1 && islower(*begin)) {
  813.   *function++ = 'V';
  814.   *function++ = *begin;
  815.   return function;
  816.  }
  817.  
  818.  /* number */
  819.  tempch = *end;
  820.  *end = '\0';
  821.  tempd=strtod(begin,(char**) &tempu);
  822.  *end = tempch;
  823.  if((char*) tempu == end) {
  824.   *function++ = 'D';
  825.   if (i_pctable < Max_ctable)
  826.     {
  827.       i_ctable[i_pctable] = tempd;
  828.       *function++ = (UCHAR) i_pctable++;
  829.     }
  830.   else
  831.     {
  832.       fset_error("too many constants");
  833.       i_error=begin;
  834.       return NULL;
  835.     }
  836.   return function;
  837.  }
  838.  
  839.  /*function*/
  840.  if(!isalpha(*begin) && *begin != '_')
  841.             /* underscores are allowed */
  842.  {
  843.   fset_error("syntax error"); 
  844.   i_error=begin;
  845.   return NULL;
  846.  }
  847.  for(endf = begin+1; endf<end && (isalnum(*endf) || *endf=='_');
  848.                                endf++);
  849.  tempch = *endf;
  850.  *endf = '\0';
  851.  if((n_function=where_table(begin)) == -1) {
  852.   *endf = tempch;
  853.   i_error=begin;
  854.   /* error message has already been created */
  855.   return NULL;
  856.  }
  857.  *endf = tempch;
  858.  if(*endf != '(' || *(end-1) != ')') {
  859.   fset_error("improper function syntax"); 
  860.   i_error=endf;
  861.   return NULL;
  862.  }
  863.  if(ftable[n_function].n_pars==0) {
  864.   /*function without parameters (e.g. pi() ) */
  865.    space=1;
  866.    for(scan=endf+1; scan<(end-1); scan++)
  867.     if(!isspace(*scan)) space=0;
  868.    if(space) {
  869.     *function++ = 'F';
  870.     *function++ = n_function;
  871.     function = comp_time(function-2,function,0);
  872.     if(fget_error()) return NULL; /* internal error in comp_time */
  873.     else return function;
  874.    } else {
  875.     i_error=endf+1;
  876.     fset_error("too many parameters");
  877.     return NULL;
  878.    }
  879.  } else {    /*function with parameters*/
  880.     tempch = *(end-1);
  881.     *(end-1) = '\0';
  882.     par_buf = (char *) malloc(strlen(endf+1)+1);
  883.     if(!par_buf)
  884.      ERROR_MEM;
  885.     strcpy(par_buf, endf+1);
  886.     *(end-1) = tempch;
  887.     /* look at the first parameter */
  888.     for(i=0; i<ftable[n_function].n_pars; i++) {
  889.      if( ( temps=my_strtok((i==0) ? par_buf : NULL) ) == NULL )
  890.       break; /* too few parameters */
  891.      paramstr[i]=temps;
  892.     }
  893.     if(temps==NULL) {
  894.      /* too few parameters */
  895.      free(par_buf);
  896.      i_error=end-2;
  897.      fset_error("too few parameters");
  898.      return NULL;
  899.     }
  900.     if((temps=my_strtok(NULL))!=NULL) {
  901.      /* too many parameters */
  902.      free(par_buf);
  903.      i_error=(temps-par_buf)+(endf+1); /* points to the first character
  904.                       of the first superfluous
  905.                       parameter */
  906.      fset_error("too many parameters");
  907.      return NULL;
  908.     }
  909.  
  910.     tempu=function;
  911.     for(i=0; i<ftable[n_function].n_pars; i++)
  912.      if(!(tempu=i_trans( tempu, paramstr[i],
  913.                  paramstr[i]+strlen(paramstr[i]) ) ) )
  914.      {
  915.       i_error=(i_error-par_buf)+(endf+1); /* moves i_error to
  916.                        the permanent copy of the
  917.                        parameter */
  918.       free(par_buf);
  919.       return NULL; /* error in one of the parameters */
  920.      }
  921.     /* OK */
  922.     free(par_buf);
  923.     *tempu++ = 'F';
  924.     *tempu++ = n_function;
  925.     tempu = comp_time(function,tempu,ftable[n_function].n_pars);
  926.     if(fget_error()) return NULL; /* internal error in comp_time */
  927.     else return tempu;
  928.  }
  929. }
  930.  
  931. /**************************************************************/
  932. /* Here begins the stand-alone part */ 
  933.  
  934.  #ifdef STAND_ALONE
  935.   /* this part of FORMULC enables it to work as a UNIX filter. */
  936.   
  937.  static char *progn = "formulc";
  938.  
  939.  
  940.  #define CtrlZ 0x001A
  941.  
  942.  static int fgetline(FILE *fpin, char s[], int size_s)
  943.  {
  944.    int         last_s, c, i;
  945.  
  946.    last_s = size_s - 1;
  947.  
  948.    i = 0;
  949.  #ifdef __MSDOS__
  950.    while ((c = getc(fpin)) != EOF && c != CtrlZ && c != '\n')
  951.  #else
  952.    while ((c = getc(fpin)) != EOF && c != '\n')
  953.  #endif
  954.      if (i < last_s) s[i++] = c;
  955.  
  956.    s[i] = '\0';
  957.  
  958.  #ifdef __MSDOS__
  959.    if (i == 0 && (c == EOF || c == CtrlZ))
  960.  #else
  961.    if (i == 0 && c == EOF)
  962.  #endif
  963.      return 0;
  964.    else
  965.      return 1;
  966.  }
  967.  
  968.  static int strarg(char *arg, char *s, int n)
  969.                                        /* - string argument -            */
  970.  {                                     /* return n-th argument of string */
  971.    char *a = arg;                      /* s in arg.                      */
  972.                                        /* n = 0 for the first argument.  */
  973.    while (*s == ' ' || *s == '\t') s++;   /* Blank and Tab are seperators.*/
  974.  
  975.    while (n--)
  976.    { while (*s != '\0' && *s != ' ' && *s != '\t') s++;
  977.      while (*s == ' ' || *s == '\t') s++;
  978.    }
  979.  
  980.    while (*s != '\0' && *s != ' ' && *s != '\t') *a++ = *s++;
  981.    *a = '\0';
  982.  
  983.    return (a != arg);
  984.  }
  985.  
  986.  static int is_legal_c_format(char *cfmt)
  987.  {
  988.    int   count, fmt_on, fmt_width;
  989.    char  *start_digits;
  990.  
  991.  
  992.    count = 0;
  993.    start_digits = NULL;
  994.  
  995.    for (fmt_on = fmt_width = 0; *cfmt; cfmt++)
  996.    {
  997.      if (isdigit(*cfmt))
  998.      {
  999.        if (start_digits == NULL)
  1000.            start_digits = cfmt;
  1001.      }
  1002.  
  1003.      if (fmt_on == 0)
  1004.      {
  1005.        if (*cfmt == '%')
  1006.        {
  1007.          fmt_on = 1;
  1008.          fmt_width = 0;
  1009.        }
  1010.      }
  1011.      else
  1012.      {
  1013.        if      (*cfmt == '%')
  1014.        {
  1015.          if (fmt_width) return 0;
  1016.          fmt_on = 0;
  1017.        }
  1018.        else if (*cfmt == '$')
  1019.        {
  1020.          if (start_digits == NULL) return 0;
  1021.          if (cfmt - start_digits != 1) return 0;
  1022.          if (*start_digits != '1') return 0;
  1023.        }
  1024.        else if (*cfmt == '*')
  1025.          return 0;
  1026.        else if (isalpha(*cfmt))
  1027.        {
  1028.          if (strchr("feEgG", *cfmt) == NULL) return 0;
  1029.          fmt_on = 0;
  1030.          count++;
  1031.        }
  1032.  
  1033.        fmt_width++;
  1034.      }
  1035.  
  1036.      if (isdigit(*cfmt) == 0)
  1037.        start_digits = NULL;
  1038.    }
  1039.  
  1040.    if (fmt_on || count != 1) return 0;
  1041.  
  1042.    return 1;
  1043.  }
  1044.  
  1045.  
  1046.  static void usage(void)
  1047.  {
  1048.    fprintf(stderr,
  1049.      "usage: %s \"function\" \"arguments\" [\"c-format\"]\n",
  1050.      progn);
  1051.    exit(1);
  1052.  }
  1053.  
  1054.  
  1055.  int main(int argc, char *argv[])
  1056.  {
  1057.    int            error, ivar, n, i;
  1058.    formu          function;
  1059.    char           *args, *e_args, buf[256], svar[sizeof buf], xtrac;
  1060.     /* args = table of variables of function;
  1061.        this table is an argument of translate and of fval */
  1062.    double         var, result;
  1063.    long           linec;
  1064.  
  1065.  
  1066.    if (argc != 3 && argc != 4) usage();
  1067.    if(argc == 4 && is_legal_c_format(argv[3]) == 0)
  1068.    {
  1069.      fprintf(stderr, "%s: Invalid format string\n", progn);
  1070.      exit(1);
  1071.    }
  1072.  
  1073.    n = strlen(argv[1]);
  1074.    if (n == 0)
  1075.      usage();   
  1076.     
  1077.    n = 0;
  1078.    e_args = argv[2];
  1079.    while (*e_args)
  1080.    {
  1081.      if (*e_args != ' ' && *e_args != '\t')
  1082.      {
  1083.        i = tolower(*e_args) - 'a';
  1084.        if (i < 0 || i > 'z' - 'a')
  1085.        {
  1086.          fprintf(stderr, "%s: Illegal argument character\n", progn);
  1087.          exit(1);
  1088.        }
  1089.  
  1090.        n++;
  1091.      }
  1092.  
  1093.      e_args++;
  1094.    }
  1095.  
  1096.    args = (char *) malloc((n + 1) * sizeof *args);
  1097.    if (args == NULL)
  1098.    {
  1099.      fprintf(stderr, "%s: Not enough core\n", progn);
  1100.      exit(1);
  1101.    }
  1102.  
  1103.    n = 0;
  1104.    e_args = argv[2];
  1105.    while (*e_args)
  1106.    {
  1107.      if (*e_args != ' ' && *e_args != '\t')
  1108.        args[n++] = tolower(*e_args);
  1109.  
  1110.      e_args++;
  1111.    }
  1112.  
  1113.    args[n] = '\0';
  1114.  
  1115.  
  1116.    function = translate(argv[1], args, &n, &error);
  1117.    if (n == 0)
  1118.    {
  1119.      fprintf(stderr, "   %s\n", argv[1]);
  1120.      for (i = 0; i < error; i++) putc('-', stderr);
  1121.      fprintf(stderr, "---^\n");
  1122.      fprintf(stderr, "%s: %s\n", progn, fget_error());
  1123.      exit(1);
  1124.    }
  1125.  
  1126.    rnd_init();
  1127.    linec = 0;
  1128.  
  1129.    while (fgetline(stdin, buf, sizeof buf))
  1130.    {
  1131.      linec++;
  1132.  
  1133.      ivar = 0;
  1134.             e_args = args;
  1135.      while(*e_args)
  1136.      {
  1137.        if (strarg(svar, buf, ivar++) == 0)
  1138.        {
  1139.          fprintf(stderr, "%s: Missing variable %d on line %ld\n",
  1140.            progn, ivar, linec);
  1141.          exit(1);
  1142.        }
  1143.  
  1144.            n = sscanf(svar, "%lf%c", &var, &xtrac);
  1145.        if (n != 1)
  1146.        {
  1147.          fprintf(stderr, "%s: Illegal variable %d on line %ld\n",
  1148.            progn, ivar, linec);
  1149.      exit(1);
  1150.        }
  1151.  
  1152.        make_var(*e_args++ , var);
  1153.      }
  1154.  
  1155.      result = fval_at(function);
  1156.  
  1157.      if (argc == 4)
  1158.       fprintf(stdout, argv[3], result);
  1159.      else
  1160.       fprintf(stdout, "%.10g", result);
  1161.  
  1162.      putc('\n', stdout);
  1163.    }
  1164.  
  1165.    return 0;
  1166.  }
  1167.  
  1168. #endif
  1169.  
  1170. /* Here is the definition of some functions in the FORMULC standard
  1171.    library */
  1172. static double pi(void)
  1173. {
  1174.  return 3.14159265358979323846264;
  1175. }
  1176.  
  1177. #ifndef MY_RND
  1178.  
  1179. void r250_init(int seed);
  1180. double dr250(void);
  1181.  
  1182. double rnd(void)
  1183. {
  1184.   return dr250();
  1185. }
  1186.  
  1187. void rnd_init(void)
  1188. {
  1189.   r250_init(time(NULL));
  1190. } 
  1191. /*The following functions were not written by Harald Helfgott. */
  1192. /******************************************************************************
  1193. *  Module:  r250.c   
  1194. *  Description: implements R250 random number generator, from S. 
  1195. *  Kirkpatrick and E.  Stoll, Journal of Computational Physics, 40, p. 
  1196. *  517 (1981).
  1197. *  Written by:    W. L. Maier
  1198. *
  1199. *    METHOD...
  1200. *        16 parallel copies of a linear shift register with
  1201. *        period 2^250 - 1.  FAR longer period than the usual
  1202. *        linear congruent generator, and commonly faster as
  1203. *        well.  (For details see the above paper, and the
  1204. *        article in DDJ referenced below.)
  1205. *
  1206. *    HISTORY...
  1207. *        Sep 92: Number returned by dr250() is in range [0,1) instead 
  1208. *            of [0,1], so for example a random angle in the
  1209. *            interval [0, 2*PI) can be calculated
  1210. *            conveniently.  (J. R. Van Zandt <jrv@mbunix.mitre.org>)
  1211. *        Aug 92: Initialization is optional.  Default condition is 
  1212. *            equivalent to initializing with the seed 12345,
  1213. *            so that the sequence of random numbers begins:
  1214. *            1173, 53403, 52352, 35341...  (9996 numbers
  1215. *            skipped) ...57769, 14511, 46930, 11942, 7978,
  1216. *            56163, 46506, 45768, 21162, 43113...  Using ^=
  1217. *            operator rather than two separate statements. 
  1218. *            Initializing with own linear congruent
  1219. *            pseudorandom number generator for portability. 
  1220. *            Function prototypes moved to a header file. 
  1221. *            Implemented r250n() to generate numbers
  1222. *            uniformly distributed in a specific range
  1223. *            [0,n), because the common expedient rand()%n is
  1224. *            incorrect.  (J. R. Van Zandt <jrv@mbunix.mitre.org>)
  1225. *        May 91: Published by W. L. Maier, "A Fast Pseudo Random Number 
  1226. *            Generator", Dr. Dobb's Journal #176.
  1227. ******************************************************************************/
  1228.  
  1229. #include <stdlib.h>
  1230.  
  1231. /**** Static variables ****/
  1232. static int r250_index = 0;
  1233. static unsigned int r250_buffer[250] = {
  1234.     15301,57764,10921,56345,19316,43154,54727,49252,32360,49582,
  1235.     26124,25833,34404,11030,26232,13965,16051,63635,55860,5184,
  1236.     15931,39782,16845,11371,38624,10328,9139,1684,48668,59388,
  1237.     13297,1364,56028,15687,63279,27771,5277,44628,31973,46977,
  1238.     16327,23408,36065,52272,33610,61549,58364,3472,21367,56357,
  1239.     56345,54035,7712,55884,39774,10241,50164,47995,1718,46887,
  1240.     47892,6010,29575,54972,30458,21966,54449,10387,4492,644,
  1241.     57031,41607,61820,54588,40849,54052,59875,43128,50370,44691,
  1242.     286,12071,3574,61384,15592,45677,9711,23022,35256,45493,
  1243.     48913,146,9053,5881,36635,43280,53464,8529,34344,64955,
  1244.     38266,12730,101,16208,12607,58921,22036,8221,31337,11984,
  1245.     20290,26734,19552,48,31940,43448,34762,53344,60664,12809,
  1246.     57318,17436,44730,19375,30,17425,14117,5416,23853,55783,
  1247.     57995,32074,26526,2192,11447,11,53446,35152,64610,64883,
  1248.     26899,25357,7667,3577,39414,51161,4,58427,57342,58557,
  1249.     53233,1066,29237,36808,19370,17493,37568,3,61468,38876,
  1250.     17586,64937,21716,56472,58160,44955,55221,63880,1,32200,
  1251.     62066,22911,24090,10438,40783,36364,14999,2489,43284,9898,
  1252.     39612,9245,593,34857,41054,30162,65497,53340,27209,45417,
  1253.     37497,4612,58397,52910,56313,62716,22377,40310,15190,34471,
  1254.     64005,18090,11326,50839,62901,59284,5580,15231,9467,13161,
  1255.     58500,7259,317,50968,2962,23006,32280,6994,18751,5148,
  1256.     52739,49370,51892,18552,52264,54031,2804,17360,1919,19639,
  1257.     2323,9448,43821,11022,45500,31509,49180,35598,38883,19754,
  1258.     987,11521,55494,38056,20664,2629,50986,31009,54043,59743
  1259.     };
  1260.  
  1261. static unsigned myrand(void);
  1262. static void mysrand(unsigned newseed);
  1263.  
  1264. /**** Function: r250_init  
  1265.     Description: initializes r250 random number generator. ****/
  1266. void r250_init(int seed)
  1267. {
  1268. /*---------------------------------------------------------------------------*/
  1269.     int        j, k;
  1270.     unsigned int mask;
  1271.     unsigned int msb;
  1272. /*---------------------------------------------------------------------------*/
  1273.     mysrand(seed);
  1274.     r250_index = 0;
  1275.     for (j = 0; j < 250; j++)     /* Fill the r250 buffer with 15-bit values */
  1276.         r250_buffer[j] = myrand();
  1277.     for (j = 0; j < 250; j++)     /* Set some of the MS bits to 1 */
  1278.         if (myrand() > 16384)
  1279.             r250_buffer[j] |= 0x8000;
  1280.     msb = 0x8000;       /* To turn on the diagonal bit   */
  1281.     mask = 0xffff;      /* To turn off the leftmost bits */
  1282.     for (j = 0; j < 16; j++)
  1283.         {
  1284.         k = 11 * j + 3;             /* Select a word to operate on        */
  1285.         r250_buffer[k] &= mask;     /* Turn off bits left of the diagonal */
  1286.         r250_buffer[k] |= msb;      /* Turn on the diagonal bit           */
  1287.         mask >>= 1;
  1288.         msb >>= 1;
  1289.         }
  1290. }
  1291.  
  1292. /**** Function: dr250 
  1293.         Description: returns a random double z in range 0 <= z < 1.  ****/
  1294. double dr250(void)
  1295. {
  1296. /*---------------------------------------------------------------------------*/
  1297.     register int    j;
  1298.     register unsigned int new_rand;
  1299. /*---------------------------------------------------------------------------*/
  1300.     if (r250_index >= 147)
  1301.         j = r250_index - 147;     /* Wrap pointer around */
  1302.     else
  1303.         j = r250_index + 103;
  1304.  
  1305.     new_rand = r250_buffer[r250_index] ^= r250_buffer[j];
  1306.  
  1307.     if (r250_index >= 249)      /* Increment pointer for next time */
  1308.         r250_index = 0;
  1309.     else
  1310.         r250_index++;
  1311.     return new_rand / 65536.;   /* Return a number in [0.0 to 1.0) */
  1312. }
  1313.  
  1314.  
  1315. /***  linear congruent pseudorandom number generator for initialization ***/
  1316.  
  1317. static unsigned long seed=1;
  1318.  
  1319.     /*    Return a pseudorandom number in the interval 0 <= n < 32768.
  1320.         This produces the following sequence of pseudorandom 
  1321.         numbers:
  1322.         346, 130, 10982, 1090...  (9996 numbers skipped) ...23369,
  1323.         2020, 5703, 12762, 10828, 16252, 28648, 27041, 23444, 6604... 
  1324.     */ 
  1325. static unsigned myrand(void)
  1326. {
  1327.     seed = seed*0x015a4e35L + 1;
  1328.     return (seed>>16)&0x7fff;
  1329. }
  1330.  
  1331.     /*    Initialize above linear congruent pseudorandom number generator */
  1332. static void mysrand(unsigned newseed)
  1333. {    seed = newseed;
  1334. }
  1335.  
  1336. #endif
  1337. /* end of random-number generator */
  1338.  
  1339.  
  1340.  
  1341.