home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / PLOT / PLOT3D_A.ZIP / CALC.C next >
Encoding:
C/C++ Source or Header  |  1991-02-24  |  18.5 KB  |  572 lines
  1. /*
  2.  
  3. Plot3d -- calc.c, the function evaluation code
  4.  
  5. By Karl Crary -- kc2m+@andrew.cmu.edu
  6.  
  7.  
  8. Copyright (c) 1991 by Karl Crary
  9.  
  10. You may use and distribute this program as much as you like so long as you do 
  11. not charge for this service.  I am not liable for failure of this program to 
  12. perform in any way.  
  13.  
  14. */
  15.  
  16.  
  17. #include "plot3d.h"
  18.  
  19.    int convert(char string[],char **errormsg,char stackno);  /* convert infix to RPN in fstack */
  20.    void process(char operator,unsigned char *findex,char cstack[],char *cstacktop,term fstack[]);
  21.    char precedence(char operator, char push);
  22.    void simplify(char stackno);
  23.    double gamma(double arg);
  24.  
  25.    extern char stop;
  26.    char ferr;
  27.  
  28. double ff(double arg1, double arg2, char stackno)
  29. {
  30.    static unsigned char begin = 0;
  31.    long double stack[25], intermediate;
  32.    unsigned char findex;
  33.    char height, ii, recurse;
  34.    if ((arg1 < -MAXDOUBLE) || (arg1 > MAXDOUBLE)) {
  35.       ferr = 1;
  36.       return(0);
  37.    }
  38.    if (range(-MINDOUBLE,arg1,MINDOUBLE)) arg1 = 0;
  39.    ferr =  0;
  40.    height = recurse = 0;
  41.    if (begin) recurse = 1;
  42.  
  43.    for ((findex = begin),(begin = 0); fstack[stackno][findex].operator != '='; findex++) {
  44.       if (!stop && kbhit()) {
  45.          ii = getch();
  46.          if (ii == ' ') stop = 2;
  47.          if (ii == 24) terminate();
  48.          if (ii == 27) {
  49.             stop = 1;
  50.             return(0);
  51.          }
  52.       }
  53.       switch (fstack[stackno][findex].operator) {
  54.          case '#' :
  55.             if (height >= 25) {
  56.                ferr = 1;
  57.                break;
  58.             }
  59.             for (ii = height; ii > 0; ii--) stack[ii] = stack[ii-1];
  60.             height++;
  61.             stack[0] = fstack[stackno][findex].operand;
  62.             break;
  63.          case 'y' :
  64.                 if (height >= 25) {
  65.                ferr = 1;
  66.                break;
  67.             }
  68.             for (ii = height; ii > 0; ii--) stack[ii] = stack[ii-1];
  69.             height++;
  70.             stack[0] = arg2;
  71.             break;
  72.  
  73.          case 'x' :
  74.             if (height >= 25) {
  75.                ferr = 1;
  76.                break;
  77.             }
  78.             for (ii = height; ii > 0; ii--) stack[ii] = stack[ii-1];
  79.             height++;
  80.             stack[0] = arg1;
  81.             break;
  82.          case '^' :
  83.             if (((stack[1] == 0) && (stack[0] <= 0)) || ((stack[1] < 0) && (stack[0] != floor(stack[0])))) ferr = 1;
  84.             else intermediate = pow(stack[1],stack[0]);
  85.             if (intermediate == HUGE_VAL) ferr = 1;
  86.             else stack[0] = intermediate;
  87.             goto popstack;
  88.          case '\\' :  /* root */
  89.             if ((stack[0] == 0) && (stack[1] <= 0)) ferr = 1;  /* division by zero */
  90.             if (stack[1] == 0) ferr = 1;  /* zeroth root */
  91.             if ((stack[0] < 0) && (floor((stack[1]+1)/2) != (stack[1]+1)/2)) ferr = 1;
  92.             if (!ferr) stack[0] = (stack[0] < 0 ? -1 : 1)*pow(fabs(stack[0]),1/stack[1]);
  93.             goto popstack;
  94.          case '*' :
  95.             stack[0] = stack[0] * stack[1];
  96.             goto popstack;
  97.          case '/' :
  98.         if (stack[0] == 0) ferr = 1;
  99.             else stack[0] = stack[1] / stack[0];
  100.             goto popstack;
  101.          case '+' :
  102.             stack[0] = stack[0] + stack[1];
  103.             goto popstack;
  104.          case '-' :
  105.             stack[0] = stack[1] - stack[0];
  106.             goto popstack;
  107.          case 'm' :
  108.             stack[0] = -stack[0];
  109.             break;
  110.          case 's' :
  111.             stack[0] = sin(stack[0]);
  112.             break;
  113.          case 'c' :
  114.             stack[0] = cos(stack[0]);
  115.             break;
  116.          case 't' :
  117.             stack[0] = tan(stack[0]);
  118.             break;
  119.          case 'S' :
  120.             if (!range(-1,stack[0],1)) ferr = 1;
  121.             else stack[0] = asin(stack[0]);
  122.             break;
  123.          case 'C' :
  124.             if (!range(-1,stack[0],1)) ferr = 1;
  125.             else stack[0] = acos(stack[0]);
  126.             break;
  127.          case 'T' :
  128.             stack[0] = atan(stack[0]);
  129.             break;
  130.          case 1 :  /* sinh */
  131.             stack[0] = sinh(stack[0]);
  132.             break;
  133.          case 2 :  /* cosh */
  134.             stack[0] = cosh(stack[0]);
  135.             break;
  136.          case 3 :  /* tanh */
  137.             stack[0] = tanh(stack[0]);
  138.             break;
  139.          case 4 :  /* invsinh */
  140.             stack[0] = log(stack[0] + sqrt(pow(stack[0],2) + 1));
  141.             break;
  142.          case 5 :  /* invcosh */
  143.             intermediate = pow(stack[0],2);
  144.             if (intermediate < 1) ferr = 1;
  145.             else {
  146.                intermediate = stack[0] + sqrt(intermediate - 1);
  147.                if (intermediate <= 0) ferr = 1;
  148.                else stack[0] = log(stack[0] + sqrt(pow(stack[0],2) - 1));
  149.             }
  150.             break;
  151.          case 6 :  /* invtanh */
  152.             if (stack[0] == 1) ferr = 1;
  153.             else {
  154.                intermediate = (1 + stack[0])/(1 - stack[0]);
  155.                if (intermediate <= 0) ferr = 1;
  156.                else stack[0] = log(intermediate)/2;
  157.             }
  158.             break;
  159.          case 'l' :
  160.             if (stack[0] <= 0) ferr = 1;
  161.             else stack[0] = log10(stack[0]);
  162.             break;
  163.          case 'n' :
  164.             if (stack[0] <= 0) ferr = 1;
  165.             else stack[0] = log(stack[0]);
  166.             break;
  167.          case 'a' :  /* absolute value */
  168.             stack[0] = fabs(stack[0]);
  169.             break;
  170.          case 'g' :  /* gamma */
  171.             stack[0] = gamma(stack[0]);
  172.             break;
  173.          case 'f' :  /* user function */
  174.         stack[0] = ff(stack[0],0,1);
  175.             break;
  176.          popstack:
  177.             height--;
  178.             for (ii = 1; ii < height; ii++) stack[ii] = stack[ii+1];
  179.             break;
  180.       }
  181.       if ((stack[0] < -2 * MAXDOUBLE) || (stack[0] > 2 * MAXDOUBLE)) ferr = 1;  /* roundoff makes MAXDOUBLE > MAXDOUBLE */
  182.       if (ferr) return(0);
  183.    }
  184.    if (recurse) begin = findex;  /* return where ended */
  185.    return(stack[0]);
  186. }
  187.  
  188.  
  189.  
  190. int matherr(struct exception *e)
  191. {
  192.    ferr = 1;
  193.    return(1);
  194. }
  195.  
  196.  
  197.  
  198. double gamma(double arg)
  199. {
  200.    double xx, total;
  201.    char ii;
  202.    /* magic gamma generating coefficients */
  203.    static double coef[7] = { 0.988205891,-0.897056937, 0.918206857,
  204.                             -0.756704078, 0.482199394,-0.193527818,
  205.                              0.035868343};
  206.  
  207.    if (!stop && kbhit()) {
  208.       ii = getch();
  209.       if (ii == ' ') stop = 2;
  210.       if (ii == 24) terminate();
  211.       if (ii == 27) {
  212.          stop = 1;
  213.          return(0);
  214.       }
  215.    }
  216.    if (arg <= 0) {
  217.       if (arg == floor(arg)) {
  218.          ferr = 1;
  219.          return(0);
  220.       }
  221.       else {
  222.          arg = -arg + 1;  /* recurse if x < 0 */
  223.          return(M_PI/(gamma(arg)*sin(M_PI*arg)));
  224.       }
  225.    }
  226.    else if (arg < 1) return(gamma(arg+1)/arg);  /* recurse if x < 1 */
  227.    else if (arg <= 2) {
  228.       arg--;
  229.       total = 1 - 0.577191652 * arg;
  230.       xx = arg;
  231.       for (ii = 0; ii <= 6; ii++) {
  232.          xx *= arg;
  233.          total += coef[ii] * xx;
  234.       }
  235.       return(total);
  236.    }
  237.    else {  /* arg > 2 */
  238.       xx = arg * arg;  /* another approximation for larger x */
  239.       total = (arg - 0.5) * log(arg);
  240.       total -= arg;
  241.       total += log(2 * M_PI) / 2;
  242.       total += 1 / (12 * arg);
  243.       total -= 1 / (360 * arg * xx);
  244.       total += 1 / (1260 * arg * xx * xx);
  245.       total -= arg / (1680 * xx * xx * xx * xx);
  246.       return(exp(total));
  247.    }
  248. }
  249.  
  250.  
  251.  
  252.  
  253.  
  254. /* errormessages */
  255.    char errorlist[6][45] = {"Values must precede each binary operator.",\
  256.                             "Unmatched right parenthesis.",\
  257.                             "Expression must end with value.",\
  258.                             "Stack overflow.",\
  259.                             "Function ends in operator.",\
  260.                             "Unbalanced parentheses."};
  261.  
  262.  
  263.  
  264. int convert(char string[], char **errormsg, char stackno)
  265. {
  266.    unsigned char findex;
  267.    char *cursor, value, inverse;
  268.    char cstack[100], cstacktop;
  269.    #define PROCESS(operator) process((operator),&findex,cstack,&cstacktop,fstack[stackno])
  270.  
  271.    cursor = string;
  272.    findex = value = inverse = 0;
  273.    cstack[0] = '(';  /* initialize cstack */
  274.    cstacktop = 0;
  275.    while (*cursor != 0) {
  276.       switch (*cursor) {
  277.          case '=' :  /* explicit definition */
  278.             cursor++;
  279.             findex = value = inverse = 0;  /* discard everything and start over */
  280.             cstack[0] = '(';
  281.             cstacktop = 0;
  282.             break;
  283.          case 'h' :  /* hyperbolics already taken care of */
  284.             cursor++;
  285.             break;
  286.          case '0' :  /* numbers */
  287.          case '1' :
  288.          case '2' :
  289.          case '3' :
  290.          case '4' :
  291.          case '5' :
  292.          case '6' :
  293.          case '7' :
  294.          case '8' :
  295.          case '9' :
  296.          case '.' :
  297.             if (value) PROCESS('$');  /* process implicit multiplication */
  298.             fstack[stackno][findex].operator = '#';  /* term is a number */
  299.             fstack[stackno][findex].operand = strtod(cursor,&cursor);  /* store number */
  300.             findex++;
  301.             value = 1;
  302.             break;
  303.          case 'x' :  /* variables */
  304.          case 'y' :
  305.          case 'o' :  /* theta */
  306.          case 'r' :
  307.          case 'u' :
  308.          case 'p' :  /* pi */
  309.          case 'e' :  /* e */
  310.          case '!' :  /* infinity */
  311.             if (value) PROCESS('$');
  312.             switch (*cursor) {
  313.                case 'p' :
  314.                   fstack[stackno][findex].operator = '#';
  315.                   fstack[stackno][findex].operand = M_PI;
  316.                   break;
  317.                case 'e' :
  318.                   fstack[stackno][findex].operator = '#';
  319.                   fstack[stackno][findex].operand = M_E;
  320.                   break;
  321.                case '!' :
  322.                   fstack[stackno][findex].operator = '#';
  323.                   fstack[stackno][findex].operand = MAXDOUBLE;
  324.                   break;
  325.                case 'y' :
  326.                   fstack[stackno][findex].operator = 'y';
  327.                   break;
  328.                default :  /* variables */
  329.                   fstack[stackno][findex].operator = 'x';
  330.                   break;
  331.             }
  332.             findex++;
  333.             cursor++;
  334.             value = 1;
  335.             break;
  336.          case '^' :  /* operators */
  337.          case '\\' :
  338.          case '*' :
  339.          case '/' :
  340.          case '+' :
  341.          case '-' :
  342.             if (!value) {  /* process unary operators */
  343.                if (*cursor == '+') {  /* unary plus, do nothing */
  344.                   cursor++;
  345.                   break;
  346.                }
  347.                else if (*cursor == '-') {  /* unary minus */
  348.                   PROCESS('m');
  349.                   cursor++;
  350.                   break;
  351.                }
  352.                if (errormsg != NULL) (*errormsg) = errorlist[0];
  353.                return(cursor-string);  /* return location of error */
  354.             }
  355.             PROCESS(*cursor);  /* process operator */
  356.             cursor++;
  357.             value = 0;
  358.             break;
  359.          case '(' :
  360.             if (value) PROCESS('$');  /* process implicit multiplication */
  361.             PROCESS('(');  /* process paren */
  362.             cursor++;
  363.             value = 0;
  364.             break;
  365.          case ')' :
  366.             if (!value) {  /* expression ends in operator */
  367.                if (errormsg != NULL) *errormsg = errorlist[2];
  368.                return(cursor-string);
  369.             }
  370.             PROCESS(')');
  371.             cstacktop--;  /* zap open paren */
  372.             if (cstacktop < 0) {
  373.                if (errormsg != NULL) (*errormsg) = errorlist[1];
  374.                return(cursor-string);  /* unbalanced parens */
  375.             }
  376.             cursor++;
  377.             value = 1;
  378.             break;
  379.          case 'd' :
  380.             if (value) PROCESS('$');  /* process implicit multiplication */
  381.             fstack[stackno][findex].operator = 'd';  /* mark start of derivative */
  382.             findex++;
  383.             PROCESS('=');  /* process end of derivative */
  384.             cursor++; /* will be followed by parenthesis */
  385.             value = 0;
  386.             break;
  387.          case 'A' :
  388.             cursor++;
  389.             /* inverse of non-trig prevented in IO */
  390.             inverse = 32;
  391.          case 's' :
  392.          case 'c' :
  393.          case 't' :
  394.             if (value) PROCESS('$');  /* process implicit multiplication */
  395.             if (cursor[1] != 'h') PROCESS(*cursor - inverse);  /* inverse is lowercase */
  396.             else switch (*cursor) {
  397.                case 's' :
  398.                   PROCESS(inverse ? 4 : 1);
  399.                   break;
  400.                case 'c' :
  401.                   PROCESS(inverse ? 5 : 2);
  402.                   break;
  403.                case 't' :
  404.                   PROCESS(inverse ? 6 : 3);
  405.                   break;
  406.             }  /* advance cursor on next pass */
  407.             cursor++;
  408.             value = inverse = 0;
  409.             break;
  410.          case ' ' :
  411.             cursor++;
  412.             break;
  413.          default :  /* functions */
  414.             if (value) PROCESS('$');  /* process implicit multiplication */
  415.             PROCESS(*cursor);
  416.             cursor++;
  417.             value = 0;
  418.             break;
  419.       }
  420.       if ((cstacktop >= 100) || (findex >= 99)) {
  421.          if (errormsg != NULL) *errormsg = errorlist[3];
  422.          return(cursor-string);  /* stack overflow */
  423.       }
  424.    }
  425.    if (!value) {
  426.       if (errormsg != NULL) *errormsg = errorlist[4];
  427.       return(cursor-1-string);  /* function ends with operator, return location error */
  428.    }
  429.    PROCESS(')');  /* flush conversion stack */
  430.    if (cstacktop != 0) {  /* unbalanced parens */
  431.       if (errormsg != NULL) *errormsg = errorlist[5];
  432.       return(cursor-string);
  433.    }
  434.    if (findex >= 99) {
  435.       if (errormsg != NULL) *errormsg = errorlist[3];
  436.       return(cursor-string);  /* stack overflow or unbalanced parens */
  437.    }
  438.    fstack[stackno][findex].operator = '=';
  439.    return(-1);
  440. }
  441.  
  442. void process(char operator,unsigned char *findex,char cstack[],char *cstacktop,term fstack[])
  443. {
  444.    char pushprec;
  445.    pushprec = precedence(operator,1);
  446.    while (precedence(cstack[*cstacktop],0) >= pushprec) {
  447.       if (cstack[*cstacktop] == '$') cstack[*cstacktop] = '*';
  448.       fstack[*findex].operator = cstack[*cstacktop];
  449.       if (++(*findex) >= 99) break;  /* stack overflow */
  450.       (*cstacktop)--;
  451.    }
  452.    if ((operator != ')') && (*cstacktop < 99)) {
  453.       (*cstacktop)++;
  454.       cstack[*cstacktop] = operator;
  455.    }
  456. }
  457.  
  458. char precedence(char operator, char push)
  459. {
  460.    switch (operator) {
  461.       case '(' :
  462.          return(push ? 10 : 0);
  463.       case ')' :
  464.          return(1);
  465.       case '=' :  /* end of derivative */
  466.          return(9);
  467.       case '$' :  /* implicit multiplication */
  468.          return(6);
  469.       case '\\' :
  470.       case '^' :
  471.          return(push ? 7 : 5);
  472.       case '*' :
  473.       case '/' :
  474.          return(3);
  475.       case '+' :
  476.       case '-' :
  477.          return(2);
  478.       default :  /* functions */
  479.          return(push ? 8 : 4);
  480.    }
  481. }
  482.  
  483. void simplify(char stackno)
  484. {
  485.    int ii, jj, restart;
  486.    char last, prevlast;
  487.    restart = -1;
  488.    last = prevlast = 0;
  489.    for (ii = 0; fstack[stackno][ii].operator != '='; ii++) switch (fstack[stackno][ii].operator) {
  490.       case '+' :  /* binary operators */
  491.       case '-' :
  492.       case '/' :
  493.       case '*' :
  494.       case '\\' :
  495.       case '^' :
  496.          if (restart == ii - 1) {
  497.             restart = ii;
  498.             break;
  499.          }
  500.          if (last && prevlast) {
  501.             /* load auxiliary function to calculate simplified value */
  502.             fstack[maxffs][0].operator = fstack[stackno][ii-2].operator;  /* might not be # */
  503.             fstack[maxffs][0].operand = fstack[stackno][ii-2].operand;
  504.             fstack[maxffs][1].operator = fstack[stackno][ii-1].operator;  /* might not be # */
  505.             fstack[maxffs][1].operand = fstack[stackno][ii-1].operand;
  506.             fstack[maxffs][2].operator = fstack[stackno][ii].operator;
  507.             fstack[maxffs][3].operator = '=';
  508.             /* simplify stack */
  509.             ii -= 2;  /* update stack position */
  510.         fstack[stackno][ii].operand = ff(0,0,maxffs);  /* calculate new value */
  511.             fstack[stackno][ii].operator = '#';
  512.             for (jj = ii+1; fstack[stackno][jj+1].operator != '='; jj++) {  /* shift stack left by 2 */
  513.                fstack[stackno][jj].operator = fstack[stackno][jj+2].operator;
  514.                fstack[stackno][jj].operand = fstack[stackno][jj+2].operand;
  515.             }
  516.             last = prevlast = 0;
  517.             ii = restart;  /* restart */
  518.             continue;
  519.          }
  520.          else last = 0;  /* prevlast irrelevant */
  521.          break;
  522.       case '#' :
  523.       case 'p' :
  524.       case 'e' :
  525.       case '!' :
  526.          prevlast = last;
  527.          last = 1;
  528.          break;
  529.       case 'x' :
  530.       case 'y' :
  531.          restart = ii;  /* no need to restart before unsimplifyable character */
  532.          last = 0;  /* prevlast irrelevant */
  533.          break;
  534.       case 'm' :
  535.          if (fstack[stackno][ii-1].operator == 'm') {  /* nuke double negate */
  536.             for (jj = ii-1; fstack[stackno][jj+1].operator != '='; jj++) {  /* shift stack left by 2 */
  537.                fstack[stackno][jj].operator = fstack[stackno][jj+2].operator;
  538.                fstack[stackno][jj].operand = fstack[stackno][jj+2].operand;
  539.             }
  540.             last = prevlast = 0;
  541.             ii = restart;
  542.             continue;
  543.          }
  544.          /* no break */
  545.       default :  /* unary operators */
  546.          if (restart == ii - 1) {
  547.             restart = ii;
  548.             break;
  549.          }
  550.          if (last) {
  551.             /* load auxiliary function to calculate simplified value */
  552.             fstack[maxffs][0].operator = fstack[stackno][ii-1].operator;  /* might not be # */
  553.             fstack[maxffs][0].operand = fstack[stackno][ii-1].operand;
  554.             fstack[maxffs][1].operator = fstack[stackno][ii].operator;
  555.             fstack[maxffs][2].operator = '=';
  556.             /* simplify stack */
  557.             ii -= 1;  /* update stack position */
  558.         fstack[stackno][ii].operand = ff(0,0,maxffs);  /* calculate new value */
  559.             fstack[stackno][ii].operator = '#';
  560.             for (jj = ii+1; fstack[stackno][jj].operator != '='; jj++) {  /* shift stack left */
  561.                fstack[stackno][jj].operator = fstack[stackno][jj+1].operator;
  562.                fstack[stackno][jj].operand = fstack[stackno][jj+1].operand;
  563.             }
  564.             last = prevlast = 0;
  565.             ii = restart;  /* restart */
  566.             continue;
  567.          }
  568.          else last = 0;  /* prevlast irrelevant */
  569.          break;
  570.    }
  571. }
  572.