home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 18 REXX / 18-REXX.zip / rximc175.zip / rxmathfn.c < prev    next >
Text File  |  2002-08-06  |  9KB  |  333 lines
  1. /* A sample function package for REXX/imc to provide the math functions */
  2. /*                                                 (C) Ian Collier 1992 */
  3.  
  4. #include<stdio.h>
  5. #include<math.h>
  6. #include<stdlib.h>
  7.  
  8. #if defined(Solaris) && !defined(__STDC__)
  9. #define const        /* Fix broken Solaris headers */
  10. #endif
  11.  
  12. #include"const.h"
  13. #include"functions.h"
  14.  
  15. /* Below is the dictionary which informs REXX what functions are available
  16.    and where to find them */
  17.  
  18. int rxsin(),rxcos(),rxtan(),rxexp(),rxatn(),rxln(),rxpower(),rxasn(),
  19.     rxacs(),rxsqr();
  20.  
  21. #ifdef sgi
  22. /* The SGI doesn't seem to like the dictionary defined below because
  23.    although the function references get relocated correctly, the
  24.    names do not.  If anyone can fix this I would like to know...
  25.    It probably involves specifying a particular flag on ld.
  26.    Instead, here is an alternative dictionary which does work. */
  27.  
  28. char names[10][8]={"ACOS","ASIN","ATAN","COS","EXP","LN","SIN",
  29.                    "SQRT","TAN","TOPOWER"};
  30. dictionary rxdictionary[]=
  31. {  names[0],rxacs,
  32.    names[1],rxasn,
  33.    names[2],rxatn,
  34.    names[3],rxcos,
  35.    names[4],rxexp,
  36.    names[5],rxln,
  37.    names[6],rxsin,
  38.    names[7],rxsqr,
  39.    names[8],rxtan,
  40.    names[9],rxpower,
  41.    0,       0
  42. };
  43.  
  44. #else
  45.  
  46. dictionary rxdictionary[]=
  47. {  "ACOS"   ,rxacs,
  48.    "ASIN"   ,rxasn,
  49.    "ATAN"   ,rxatn,
  50.    "COS"    ,rxcos,
  51.    "EXP"    ,rxexp,
  52.    "LN"     ,rxln,
  53.    "SIN"    ,rxsin,
  54.    "SQRT"   ,rxsqr,
  55.    "TAN"    ,rxtan,
  56.    "TOPOWER",rxpower,
  57.    0        ,0
  58. };
  59.  
  60. #endif /* sgi */
  61.  
  62. /* Two utility functions not provided by the REXX calculator:
  63.    getdouble(argc)   checks that argc==1 and then gets that argument from the
  64.                      stack, returning it as a double.
  65.    stackdouble(argc) stacks the given double on the calculator stack and
  66.                      formats it according to the current settings */
  67.  
  68. double getdouble(argc)
  69. int argc;
  70. {
  71.    char *num;
  72.    char c;
  73.    int len;
  74.    double ans;
  75.    if(argc!=1)die(Ecall);
  76.    num=delete(&len);
  77.    if(len<0)die(Ecall);
  78.    
  79.    /* now remove all spaces (which are legal in REXX) for sscanf */
  80.    while(len&&num[0]==' ')num++,len--;
  81.    if(num[0]=='+'||num[0]=='-')
  82.       while(len>1&&num[1]==' ')num[1]=num[0],num++,len--;
  83.    while(len&&num[len-1]==' ')len--;
  84.  
  85.    /* Now the number is converted to double with sscanf */
  86.    num[len]=0;
  87.    if(sscanf(num,"%lf%c",&ans,&c)!=1)die(Enum);
  88.    if(ans==HUGE_VAL)die(Eoflow);
  89.    return ans;
  90. }
  91.  
  92. void stackdouble(num)
  93. double num;
  94. {  /* the number is simply sprintf-ed, stacked, and formatted. */
  95.    void rxformat();    /* the REXX format() function */
  96.    char result[25];
  97.    sprintf(result,"%.18G",num);
  98.    stack(result,strlen(result));
  99.    rxformat(1);
  100. }
  101.  
  102. /* The following functions implement the REXX math functions using
  103.    floating-point math. Some functions perform range checking first. */
  104.    
  105. int rxsin(name,argc)
  106. char *name;
  107. int argc;
  108. {
  109.    stackdouble(sin(getdouble(argc)));
  110.    return 1;
  111. }
  112. int rxcos(name,argc)
  113. char *name;
  114. int argc;
  115. {
  116.    stackdouble(cos(getdouble(argc)));
  117.    return 1;
  118. }
  119. int rxtan(name,argc)
  120. char *name;
  121. int argc;
  122. {
  123.    double arg=getdouble(argc);
  124.    double s=sin(arg);
  125.    double c=cos(arg);
  126.    if(c==0.)return -Edivide;
  127.    stackdouble(s/c);
  128.    return 1;
  129. }
  130. int rxasn(name,argc)
  131. char *name;
  132. int argc;
  133. {
  134.    double arg=getdouble(argc);
  135.    if(arg>1.||arg<-1.)return -Ecall;
  136.    stackdouble(asin(arg));
  137.    return 1;
  138. }
  139. int rxacs(name,argc)
  140. char *name;
  141. int argc;
  142. {
  143.    double arg=getdouble(argc);
  144.    if(arg>1.||arg<-1.)return -Ecall;
  145.    stackdouble(acos(arg));
  146.    return 1;
  147. }
  148. int rxatn(name,argc)
  149. char *name;
  150. int argc;
  151. {
  152.    stackdouble(atan(getdouble(argc)));
  153.    return 1;
  154. }
  155. int rxexp(name,argc)
  156. char *name;
  157. int argc;
  158. {
  159.    double val=exp(getdouble(argc));
  160.    if(val==HUGE_VAL)return -Eoflow;
  161.    stackdouble(val);
  162.    return 1;
  163. }
  164. int rxln(name,argc)
  165. char *name;
  166. int argc;
  167. {
  168.    double arg=getdouble(argc);
  169.    if(arg<=0.)return -Ecall;
  170.    arg=log(arg);
  171.    if(arg==HUGE_VAL)return -Eoflow;
  172.    stackdouble(arg);
  173.    return 1;
  174. }
  175. int rxpower(name,argc)
  176. char *name;
  177. int argc;
  178. {
  179.    double arg1,arg2;
  180.    if(argc!=2)return -Ecall;
  181.    arg2=getdouble(1);
  182.    arg1=getdouble(1);
  183.    if(arg2==0.){        /* Firstly, anything to the power 0 is 1. */
  184.       stack("1",1);
  185.       return 1;
  186.    }
  187.    if(arg1==0.&&arg2>0.){  /* Also, 0 to any positive power is 0. */
  188.       stack("0",1);
  189.       return 1;
  190.    }                    /* Otherwise return x**y = exp(y log x) */
  191.    if(arg1<0.)return -Ecall;  
  192.    arg1=exp(arg2*log(arg1));
  193.    if(arg1==HUGE_VAL)return -Eoflow;
  194.    stackdouble(arg1);
  195.    return 1;
  196. }
  197.  
  198. /* Two utility functions for the square root program below.  They provide
  199.    add and subtract functions which are quicker than going via REXX's
  200.    calculator. */
  201.    
  202. static int add(num1,len1,num2,len2)  /* add num 2 to num 1. */
  203. char *num1,*num2;                    /* Return the carry */
  204. int len1,len2;
  205. {
  206.    int i;
  207.    char c=0;
  208.    num1+=len1;
  209.    num2+=len2;
  210.    for(i=len2;i--;){
  211.       *--num1+=*--num2+c-'0';
  212.       if(c=(*num1>'9'))*num1-=10;
  213.    }
  214.    for(i=len1-len2;c&&i--;)
  215.       if(c=((++*--num1)>'9'))*num1-=10;
  216.    if(c)*--num1='1';
  217.    return c;
  218. }
  219.  
  220. static void sub(num1,len1,num2,len2)  /* subtract num2 from num1. */
  221. char *num1,*num2;
  222. int len1,len2;
  223. {
  224.    int i;
  225.    char c=0;
  226.    num1+=len1;
  227.    num2+=len2;
  228.    for(i=len2;i--;){
  229.       *--num1-=*--num2+c-'0';
  230.       if(c=(*num1<'0'))*num1+=10;
  231.    }
  232.    for(i=len1-len2;c&&i--;)
  233.      if(c=((--*--num1)<'0'))*num1+=10;
  234. }
  235.  
  236. int rxsqr2(name,argc)  /* Quick square root evaluator, for "low" precision */
  237. char *name;
  238. int argc;
  239. {
  240.    double arg=getdouble(argc);
  241.    if(arg<0.)return -Ecall;
  242.    stackdouble(sqrt(arg));
  243.    return 1;
  244. }
  245.  
  246. /* The following function implements the REXX square root function.
  247.    If the number of digits required is low enough, the floating-point
  248.    method (above) is used.  Otherwise, a digit-by-digit method is used. */
  249.    
  250. int rxsqr(name,argc)
  251. char *name;
  252. int argc;
  253. {
  254.    int n,m,e,l,z;
  255.    char x;
  256.    char *arg,*ans,*residue,*string,*tmp1,*tmp2;
  257.    int rlen,t1len,t2len;
  258.    extern int precision;  /* Extract the precision and the workspace from */
  259.    extern char *workptr;  /* the REXX calculator */
  260.    extern int eworkptr,worklen;
  261.    if(precision<=16)return rxsqr2(name,argc);
  262.    /* now the digit by digit algorithm for greater precision than FP math */
  263.    eworkptr=1;  /* Clear the work space, but allow for one digit of left-
  264.                    overflow (to be put at workptr[0]) */
  265.    if((n=num(&m,&e,&z,&l))<0)return -Enum;  
  266.    if(m)return -Ecall;
  267.    delete(&m);
  268.    if(z){       /* the square root of zero is zero. */
  269.       stack("0",1);
  270.       return 1;
  271.    }
  272.    /* workptr[n] is now the first in a sequence of digits of a positive
  273.       number.  The decimal point belongs after the first digit.  The
  274.       following code moves the decimal point to after the second digit,
  275.       and at the same time makes the exponent even. */
  276.    if(e%2)e--;
  277.    else workptr[--n]='0', /* (add leading 0) */
  278.         l++;
  279.    /* Expand the workspace to hold all the temporary material */
  280.    mtest(workptr,worklen,n+l+6*precision+20,n+l+6*precision+20-worklen);
  281.    arg=workptr+n;
  282.    if(l%2)arg[l++]='0'; /* make the length even by adding a trailing 0 */
  283.    ans=arg+l;
  284.    residue=ans+precision+1;
  285.    rlen=0;
  286.    string=residue+precision+precision+2;
  287.    tmp1=string+precision+2;
  288.    tmp2=tmp1+precision+3;
  289.    string[0]='0';
  290.    /* Now we have five numbers: ans, residue, string, tmp1, tmp2, all of
  291.       ample length, in the workspace as well as the argument. */
  292.       
  293.    for(n=0;n<=precision;n++){
  294.       /* invariants: length of arg used so far = 2n
  295.                      length of ans = n
  296.              length of string = n+1
  297.              length of residue = rlen <= 2n       */
  298.              
  299.       if(n+n<l)residue[rlen++]=arg[n+n],  /* Bring down two more digits */
  300.                residue[rlen++]=arg[n+n+1];/* from the argument to residue */
  301.       else residue[rlen]=residue[rlen+1]='0',
  302.            rlen+=2;
  303.        
  304.       /* find least x s.t. (string||x)*x>residue. Find (string||x-1)*(x-1). */
  305.       /* Invariant wil be: tmp1=(string||x-1)*(x-1) and tmp2=(string||x)*x  */
  306.       memcpy(tmp2,string,n+1);
  307.       tmp2[n+1]='1';
  308.       memset(tmp1,'0',t1len=t2len=n+2);
  309.       x='1';
  310.       /* Remove some leading zeros from the residue */
  311.       while(rlen>t2len&&residue[0]=='0')residue++,rlen--;
  312.       /* Below says "while(tmp2<=residue)" */
  313.       while(t2len<rlen||(t2len==rlen&&memcmp(tmp2,residue,t2len)<=0)){
  314.          memcpy(tmp1,tmp2,t1len=t2len);
  315.      string[n+1]=x;
  316.      if(add(tmp2,t2len,string,n+2))tmp2--,t2len++;
  317.      x++;
  318.      if(add(tmp2,t2len,&x,1))tmp2--,t2len++;
  319.       }
  320.       if(t2len>n+2)tmp2++; /* Restore tmp2 pointer to its original value */
  321.       
  322.       /* Append x-1 to the answer; subtract (string||x-1)*(x-1) from residue;
  323.          append 2(x-1) to the string. */
  324.       ans[n]=--x;
  325.       sub(residue,rlen,tmp1,t1len);
  326.       string[n+1]=x+x-'0';
  327.       if(x>='5')string[n+1]-=10,string[n]++;
  328.    }
  329.    stacknum(ans,n,e/2,0); /* We have the answer! */
  330.    return 1;
  331. }
  332.  
  333.