home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_100 / 107_01 / clogs.c < prev    next >
Text File  |  1984-06-05  |  18KB  |  652 lines
  1.  
  2.  
  3.   /*  CLOGS.C       *****
  4. A group of programs in C, using the BDS-C floating point package
  5. as modified by LCC (FLOAT+44.*) and depending on the ability to
  6. insert null characters in a string available in BDS-C V 1.44.
  7. Four functions are handled:
  8.     log10, exp10, expe, pi
  9. In addition, there is a service function exprange() which returns
  10. a false (00) if the exponent of a floating point variable is reaching
  11. near the end of the usable range.
  12. These are discussed in detail in CLOGS.DOC
  13.  
  14. L. C. Calhoun
  15. 257 South Broadway
  16. Lebanon, Ohio 45036  513 day 433 7510 nite 932 4541
  17.  
  18.   */
  19.  
  20.  
  21. char *pi(result)
  22.  
  23.  /* result is a standard character array char result[5]; in calling
  24. program as used for floating point variables.  The return is a pointer
  25. to result, and the value of pi stored in the result array in floating
  26. point */
  27.  
  28. char *result;
  29. {
  30.  char *piconst, *fpasg();
  31.  
  32.  piconst = "\171\356\207\144\2";
  33.  return (fpasg(result,piconst) );
  34. }
  35.  
  36. char *expe(result,xin)
  37.  
  38.  /* computes the base of the natural log "e" raised to the "x'th"
  39.     power.  Checks are made for out of range values and result is
  40.     defaulted to 0, 1, or a large number as appropriate.  There are
  41.     no error flags.  The arguments are floating point character
  42.     arrays char result[5], x[5]; in calling program.  Return is
  43.     a pointer to result, and "e to the x" stored in result.
  44.  */
  45.  
  46. char *result, *xin;
  47. {
  48.     char *zero, *one, *large, *coef[7], *eghty6, *meghty6;
  49.     char intres[5], xint[5], x[5];
  50.     char  *fpmult(), *fpasg(), *fpdiv(), *fpchs();
  51.     int signx, index, bigexp, smallexp, zeroin;
  52.     int exprange();
  53.  
  54.     bigexp = smallexp = zeroin = 0;
  55.  
  56.     zero ="\0\0\0\0\0";
  57.     one = "\0\0\0\100\1";
  58.     large = "\0\0\0\100\175";
  59.     eghty6 = "\0\0\0\126\7";
  60.     meghty6 = "\0\0\0\252\7";
  61.  
  62.     coef[0] = "\0\0\0\100\1";
  63.     coef[1] = "\140\326\377\177\376";
  64.     coef[2] = "\130\373\3\100\374";
  65.     coef[3] = "\200\1\352\124\370";
  66.     coef[4] = "\351\253\362\131\364";
  67.     coef[5] = "\21\213\32\133\357";
  68.     coef[6] = "\371\330\260\134\354";
  69.  
  70.   /* preserve input datum */
  71.     fpasg(x,xin);
  72.  
  73.  /* check for sign */
  74.     if (x[3] < 128)   /* positive */
  75.        {signx = 1;
  76.          }
  77.     else
  78.          {signx = 0;
  79.           fpchs(x,x);
  80.          }
  81.  
  82.    /* check for zero and out of range of fp var */
  83.  
  84.    /* check for zero and very small numbers */
  85.     if ( ((x[4]>127) && (x[4]<226)) || ( (x[4]==0) &&
  86.         (x[3]==0) && (x[2]==0) && (x[1]==0) && (x[0]==0) ) )
  87.         zeroin = 1;
  88.  
  89.    /* check for very small exponent */
  90.     if ( fpcomp(xin,meghty6) < 0)
  91.         smallexp = 1;
  92.  
  93.   /*  check for very large exponent */
  94.     if ( fpcomp(x,eghty6) > 0 )
  95.         bigexp = 1;
  96.  
  97.   /*  now if small number or zero, result is one */
  98.   /*  now if big number and positive, result is large number */
  99.   /*  now if big number and negative, result is zero */
  100.  
  101.     if (zeroin) return (fpasg(result,one) );
  102.     if (smallexp) return (fpasg(result,zero) );
  103.     if (bigexp) return (fpasg(result,large) );
  104.  
  105.  /*  all exceptions taken care of, so evaluate rest  */
  106.         fpasg(result,one);
  107.         fpasg(xint,x);
  108.         index = 1;
  109.         while ( (index<=6) && exprange(xint) )
  110.            {
  111.            fpmult(intres,coef[index],xint);
  112.         fpadd(result,result,intres);
  113.         fpmult(xint,xint,x);
  114.         index++;
  115.            }
  116.   /* now do the square square */
  117.     fpmult(result,result,result);
  118.     fpmult(result,result,result);
  119.  
  120.   /* now treat sign appropriately */
  121.     if (signx) return (result);
  122.     else
  123.        {fpdiv(result,one,result);
  124.         return (result);
  125.        }
  126. }
  127.  
  128. char *exp10(result,xin)
  129.  
  130. /* similar to expe, except result returned is 10 raised to the x
  131. power    the antilogarithm to base 10  */
  132.  
  133. char *result, *xin;
  134. {
  135.     char *zero, *ten, *one, *large, *thty8;
  136.     char xint[5], *coef[7], intres[5], tenfac[5], x[5];
  137.     int index, bigexp, smallexp, signx, tenpower;
  138.     int exprange();
  139.  
  140.     bigexp = smallexp = 0;
  141.  
  142.     zero ="\0\0\0\0\0";
  143.     one = "\0\0\0\100\1";
  144.     large = "\0\0\0\100\175";
  145.     ten = "\0\0\0\120\4";
  146.     thty8 = "\0\0\0\114\6";
  147.  
  148.     coef[1] = "\65\264\256\111\1";
  149.     coef[2] = "\0\14\330\124\0";
  150.     coef[3] = "\0\46\354\100\377";
  151.     coef[4] = "\24\140\107\115\375";
  152.     coef[5] = "\242\304\361\155\372";
  153.     coef[6] = "\361\143\246\134\371";
  154.  
  155.  /* preserve input datum */
  156.     fpasg(x,xin);
  157.  
  158.  /* check for sign */
  159.     if (x[3] < 128)   /* positive */
  160.        signx = 1;
  161.     else     /* negative */
  162.        {signx = 0;
  163.         fpchs(x,x);
  164.        }
  165.  
  166.   /* check for very small or large numbers, check by exponent size */
  167.   /* check for zero or small */
  168.     if ( ((x[4]>127) && (x[4]<226)) || ( (x[4]==0) &&
  169.         (x[3]==0) && (x[2]==0) && (x[1]==0) && (x[0]==0) ) )
  170.         smallexp = 1;
  171.    /* check for big number */
  172.     if ( fpcomp(x,thty8) > 0)
  173.         bigexp = 1;
  174.  
  175.  /* if value is small or zero, return 1 as with expe */
  176.  /* if value is large and positive, return a large number */
  177.  /* if value is large and negative, return zero */
  178.  
  179.     if (smallexp) return (fpasg(result,one) );
  180.  
  181.     if(bigexp && signx) return (fpasg(result,large) );
  182.  
  183.     if(bigexp && !signx) return (fpasg(result,zero) );
  184.  
  185.  /* now reduce range of x to between zero and one */
  186.  
  187.     tenpower = ftoit(x);
  188.     itof(tenfac,tenpower);
  189.     fpsub(x,x,tenfac);
  190.     fpasg(tenfac,one);
  191.     while (tenpower)
  192.         {fpmult(tenfac,tenfac,ten);
  193.          tenpower--;
  194.         }
  195.  /* now evaluate series  */
  196.     fpasg(result,one);
  197.     fpasg(xint,x);
  198.     index = 1;
  199.  
  200.         while ( (index <= 6) && exprange(xint) )
  201.            {fpmult(intres,coef[index],xint);
  202.         fpadd(result,result,intres);
  203.         fpmult(xint,xint,x);
  204.         index += 1;
  205.            }
  206.  
  207.  /* now square result (note error in referenced article) */
  208.     fpmult(result,result,result);
  209.  
  210.  /* now check sign and make proper return */
  211.     fpmult(result,result,tenfac);  /* scale back up */
  212.  
  213.     if (signx) return (result);
  214.  
  215.     else return ( fpdiv(result,one,result) );
  216. }
  217.  
  218.  
  219. char *log10(result,sign,xin)
  220.  
  221.  /* computes briggsian logarithm of x which is a char[5]
  222.     floating point number.  Return is logarithm in result[5],
  223.     and sign pointed to by sign.  The logarithm is taken
  224.     of the magnitude, and sign information preserved
  225.     as required by sign.
  226.  */
  227. char *result, *xin;
  228. int *sign;
  229. {
  230.     char *zero, *ten, *one, *large;
  231.     char *sqrtten, x[5];
  232.     char gamma[5], gamnum[5], gamden[5], *coef[5];
  233.     char *half;
  234.     char intres[5], gamint[5];
  235.     int tenpower;
  236.     int index, bigexp, smallexp, signx;
  237.     int exprange();
  238.  
  239.     bigexp = smallexp = 0;
  240.  
  241.     zero ="\0\0\0\0\0";
  242.     one = "\0\0\0\100\1";
  243.     large = "\0\0\0\114\6";
  244.     ten = "\0\0\0\120\4";
  245.     half = "\0\0\0\100\0";
  246.     sqrtten = "\304\145\61\145\2";
  247.  
  248.     coef[0] = "\362\6\56\157\0";
  249.     coef[1] = "\30\346\21\112\377";
  250.     coef[2] = "\100\55\344\132\376";
  251.     coef[3] = "\106\73\244\140\375";
  252.     coef[4] = "\174\5\367\141\376";
  253.  
  254.  /*  preserve input variable */
  255.     fpasg(x,xin);
  256.  
  257.  /* check for sign */
  258.     if (x[3] < 128)   /* positive */
  259.        signx = 1;
  260.     else     /* negative */
  261.        {signx = -1;
  262.         fpchs(x,x);
  263.        }
  264.  
  265.   /* check for very small or large numbers, check by exponent size */
  266.   /* check for zero or small */
  267.     if ( ((x[4]>127) && (x[4]<209)) || ( (x[4]==0) &&
  268.         (x[3]==0) && (x[2]==0) && (x[1]==0) && (x[0]==0) ) )
  269.         smallexp = 1;
  270.    /* check for big number */
  271.     if ( (x[4]  >47) && (x[4] < 128) )
  272.         bigexp = 1;
  273.  
  274.  /* if very small, return - a large number
  275.     if very large, return + a large number  */
  276.  
  277.     if (smallexp)
  278.         {*sign = signx;
  279.          return (fpchs(result,large) );
  280.         }
  281.     if (bigexp)
  282.         {*sign = signx;
  283.          return (fpasg(result,large) );
  284.         }
  285.   /*  now bring into range 1 <= x < 10 */
  286.     tenpower = 0;
  287.     while ( fpcomp(x,ten) >= 0)
  288.            {tenpower++;
  289.         fpdiv(x,x,ten);
  290.            }
  291.     while ( fpcomp(x,one) < 0)
  292.            {tenpower--;
  293.         fpmult(x,x,ten);
  294.            }
  295.  
  296.   /* now if exactly one, no need to evaluate */
  297.     fpsub(gamnum,x,one);
  298.     if (((gamnum[4]>127)&&(gamnum[4]<209)) || ((gamnum[0]==0) &&
  299.        (gamnum[1]==0) && (gamnum[2] == 0) && (gamnum[3] == 0) ) )
  300.            {*sign = signx;
  301.         itof(result,tenpower);
  302.         return (result);
  303.            }
  304.  /* now compute gamma  for series  */
  305.  
  306.     fpsub(gamnum,x,sqrtten);
  307.   /* now check for size of numerator */
  308.     if (((gamnum[4]>127)&&(gamnum[4]<209)) || ((gamnum[0]==0) &&
  309.        (gamnum[1]==0) && (gamnum[2] == 0) && (gamnum[3] == 0) ) )
  310.            {itof(result,tenpower);
  311.         fpadd(result,result,half);
  312.         *sign = signx;
  313.         return (result);
  314.            }
  315.     fpadd(gamden,x,sqrtten);
  316.     fpdiv(gamma,gamnum,gamden);
  317.  
  318.  /* now set up for series (use gamnum as gamma squared) */
  319.     fpmult(gamnum,gamma,gamma);
  320.     fpasg(gamint,gamma);
  321.     index = 0;
  322.     fpasg(result,half);
  323.  
  324.  /* now do series evaluation */
  325.     while ( (index <= 4) && exprange(gamint) )
  326.            {fpmult(intres,coef[index],gamint);
  327.         fpadd(result,result,intres);
  328.         fpmult(gamint,gamint,gamnum);
  329.         index++;
  330.            }
  331.  
  332.  /* now do correction for range reduction */
  333.     if (tenpower != 0)
  334.           {itof(intres,tenpower);
  335.         fpadd(result,result,intres);
  336.           }
  337.  
  338.  /* now clean up and return */
  339.  
  340.     *sign = signx;
  341.     return (result);
  342. }
  343.  
  344. int exprange(x)
  345.  
  346.  /* The input argument is a floating point function from BDS C which
  347. consists of an array of 5 character data.  The function returns
  348. a 1 if the exponent is in the range of - 47 to + 47.  Outside this
  349. range a value of 0 (false) is returned.  This is a range of ten to
  350. plus or minus 14 power */
  351. char *x;
  352.  
  353. { if ( ((x[4]<128) && (x[4]>47) ) || ((x[4]>127) && (x[4] < 209) ) )
  354.     return (0);
  355.   else return (1);
  356. }
  357. r+s├┘%!É    n}╖╩7&!Æ    ^#Vr+s!═╥7&!₧    ^#Vr+sδ6 ├&├M&!₧    ^#Vr+s╒`inδßs├c&!₧    ^#Vr+s╒`inδßs├╕!₧    ~#fo6δ!Ü9∙δ┴╔├&├ì/┼!9∙DM═Én}╖╩⌐&!    ^#Vr+sδn&σ═|&╤├ç&┴╔├▒&├─*┼!9∙DM!    ~#foσ!    ~#foσ!σ═«&╤╤╤!    ~#fo┴╔┴╔┼!9∙DM!    n&0═∙═»┌ '!    
  358.  
  359.  /* TESTCLOG.C ****** testing program for
  360.     CLOGS.C testing pi, expe, exp10, and log10   functions
  361.     calls for entry of a number, multiplies it by pi,
  362.     obtains log10 of number and of pi times
  363.     obtains exp10 of two logs
  364.     obtains exp10 and expe of number and pi product
  365.   */
  366.  
  367.  main()
  368. {
  369.     char rlog[5], rlogpi[5], xin[5], x[5];
  370.     char rpi[5], rexp10[5], rexp10pi[5];
  371.     char rexpe[5], rexpepi[5], rexlog[5];
  372.     char rexlogpi[5];
  373.     char *log10(), *expe(), *pi(), *exp10();
  374.     char stringin[80], pival[5];
  375.      int *sign1, *sign2;
  376.  
  377. while(1)
  378.     {printf("\n Enter input value ");
  379.      gets(stringin);
  380.      atof(xin,stringin);
  381.      fpasg(x,xin);
  382.      pi(pival);
  383.      fpmult(rpi,pival,xin);
  384.      printf("\n Input %e pi %e times pi %e",x,pival,rpi);
  385.      expe(rexpe,x);
  386.      expe(rexpepi,rpi);
  387.      printf("\n expe %e expe of pi times %e",rexpe,rexpepi);
  388.      fpasg(x,xin);
  389.      exp10(rexp10,x);
  390.      exp10(rexp10pi,rpi);
  391.      printf("\n exp10 %e exp10 of pi times %e",rexp10,rexp10pi);
  392.      fpasg(x,xin);
  393.      log10(rlog,sign1,x);
  394.      log10(rlogpi,sign2,rpi);
  395.      printf("\n log10 %f sign %d log10 of pi %f sign %d",
  396.         rlog,*sign1,rlogpi,*sign2);
  397.      fpasg(x,xin);
  398.      exp10(rexlog,rlog);
  399.      exp10(rexlogpi,rlogpi);
  400.      if(*sign1 < 0) fpchs(rexlog,rexlog);
  401.      if(*sign2 < 0) fpchs(rexlogpi,rexlogpi);
  402.      printf("\n In %e PIVAL %e unlog %e unlogpi %e \n",
  403.         x, rpi, rexlog, rexlogpi);
  404.     }
  405. }
  406. ;
  407. {
  408.     char *zero, *one, *large, *coef[7], *eghty6, *meghty6;
  409.     char intres[5], xint[5], x[5];
  410.     char  *fpmult(), *fpasg(), *fpd
  411.  
  412.  
  413.  
  414.  
  415.  
  416.                       NOTES ON "LOG" FUNCTIONS FOR BDS - C
  417.  
  418.  
  419.                                   Introduction
  420.  
  421.  
  422.       These "log" functions were developed so that I could do some auxiliary
  423.  
  424. work on scaling and curve generation for a graphics package I'm now doing in
  425.  
  426. BDS-C.  They seem to work OK in my setup which is now:  
  427.  
  428.        o  Altair 8800b, 64K CPM 2.2 BDS-C vers 1.44 
  429.  
  430.        o  Tarbell SSSD 4 8 inch disk 
  431.  
  432.        o  Scion Microangelo 
  433.  
  434.        o  LSI ADM3A 
  435.  
  436. This package was developed by:  
  437.  
  438. L. C. Calhoun PE 
  439.  
  440. 257 South Broadway 
  441.  
  442. Lebanon, Ohio 45036 
  443.  
  444. <513> 932-4541/433-7510 
  445.  
  446.                        SPECIAL NOTE ON VERSION OF BDS - C
  447.  
  448.  
  449.       The CLOGS programs have been written to take advantage of the ability to
  450.  
  451. insert '\0' into string constants which make it possible to use string constant
  452.  
  453. as pseudo-static floating point constants.  
  454.  
  455.  
  456.  
  457.       These programs are written in BDS-C using the floating point package
  458.  
  459. modified to add truncation and magnitude functions.  This package is called
  460.  
  461. "FLOAT+44".  The following functions are mechanized:  
  462.  
  463.      char *expe(result,x) 
  464.  
  465.      char *result, *x;   /* usual [5] char arrays for fp */ 
  466.  
  467.           The program returns the base of natural logs "e" raised to the power
  468.  
  469.           given in x.  The program is limited to work within the bounds of the
  470.  
  471.  
  472.  
  473.  
  474.  
  475.                                         1
  476.  
  477.  
  478.  
  479.  
  480.  
  481.  
  482.  
  483.  
  484.  
  485.           floating point variable.  The function returns the pointer to the
  486.  
  487.           result.  Values outside floating point bounds for the result are set
  488.  
  489.           to either zero or a very large number on the order of 2e38.  No error
  490.  
  491.           flags exist or are set.  
  492.  
  493.      char *exp10(result,x) 
  494.  
  495.      char *result, *x;   /* as with expe */ 
  496.  
  497.           Identical to expe, except the base of briggs logs "10" is raised to
  498.  
  499.           the power indicated by x.  
  500.  
  501.      char *log10(result,sign,x) 
  502.  
  503.      char *result, *x;    /*as with expe */ 
  504.  
  505.      int *sign; 
  506.  
  507.      char *angle, *datum;    
  508.  
  509.           This returns the logarithm to the base 10 in result of the value in x.
  510.  
  511.           x is unchanged.  Logarithms are computed of the magnitude of x, and
  512.  
  513.           negative x values return a -1 in sign.  Positive x values produce a 1
  514.  
  515.           in sign.  Very large or small values are returned for out   of range
  516.  
  517.           data.  There is no over/underflow indication.  Also returns pointer to
  518.  
  519.           result.  
  520.  
  521.                                      Method
  522.  
  523.  
  524.       The methods used are outlined in "Functional Approximations" by Fred
  525.  
  526. Ruckdeschel; page 34 ff in BYTE for November 1978.  Note the corrections in the
  527.  
  528. January 1979 issue.  A number of references are given in that article, and are
  529.  
  530. recommended reading.  There is an error in Ruckdeschel' article, in Table 3b.
  531.  
  532. The term to the right of the equals should all be enclosed in brackets and
  533.  
  534. squared.  Refer to equation 4.2.47 in Ruckdeshel' reference 6.  The following
  535.  
  536. service function is used:  
  537.  
  538.      int exprange(x) 
  539.  
  540.  
  541.  
  542.  
  543.  
  544.                                         2
  545.  
  546.  
  547.  
  548.  
  549.  
  550.  
  551.  
  552.  
  553.  
  554.      char *x; 
  555.  
  556.           This is used in the series evaluations.  Input is a pointer to the
  557.  
  558.           BDS-C type floating point variable.  It returns a 1 (true) if -47 <=
  559.  
  560.           exponent <= 47.  Outside of the range of exponent a 0 (false) is
  561.  
  562.           returned.  The exponent is a power of 2, so the effective range is
  563.  
  564.           about 1.4e14 to 7.e-15 for exponent in decimals.  The function is used
  565.  
  566.           to avoid a series computation overflowing the exponent which wraps a
  567.  
  568.           small exponent into a large one and vice-versa.  
  569.  
  570.  
  571.  
  572.       A number of checks are made for very large and very small data, to protect
  573.  
  574. the evaluation from the underflow and overflow failures of the floating point
  575.  
  576. package.  I have used TESTCLOG to evaluate over a wide range of variables, and I
  577.  
  578. think! I got all the gotchas.  It looks as though the package has (except at
  579.  
  580. infinity) about an absolute accuracy of .00001.  You will note that I use the
  581.  
  582. properties of the floating point numbers to do magnitude and sign checks.  Lots
  583.  
  584. faster than using fpcomp().  There is another program included, COEFSTAT, which
  585.  
  586. I used to derive the five octal equivalents for the pseudo-static terms in the
  587.  
  588. series evaluations, etc.  This will only work with BDS C V 1.44 (and later, I
  589.  
  590. hope) which allows insertion of nulls ('\0') in string constants.  The string
  591.  
  592. constants are used as pseudo-static floating point constants..and work very
  593.  
  594. well.  
  595.  
  596.  
  597.  
  598.                         Components of the CLOGS Package
  599.        1. CLOGS.DOC             This documentation file 
  600.  
  601.        2. CLOGS.C               Source for trig package 
  602.  
  603.        3. COEFSTAT.C            Source for coeficient determing program 
  604.  
  605.        4. TESTCLOG.C            Source for trig function testing program 
  606.  
  607.  
  608.  
  609.  
  610.  
  611.  
  612.  
  613.                                         3
  614.  
  615.  
  616.  
  617. power);
  618.     fpsub(x,x,tenfac);
  619.     fpasg(tenfac,one);
  620.     while (tenpower)
  621.         Software contributions are received for inclusion into the
  622.   library with the understanding that the contributor is
  623.   authorized to make the material available to others for their
  624.   individual, non-commercial use.  The Users Group makes no
  625.   representations as to the utility of the material in the
  626.   library for any purpose.  Contributions should be submitted
  627.   on 8" single density diskettes in CP/M file form.  Please
  628.   cross reference any rewrites or bug-fixes to prior
  629.   distributions
  630.  
  631. National CP/M Users Group - Program Submission Form
  632.  
  633. Submission Date: 27 July 1981 
  634.  
  635. Files names:    FLOAT+44.*  Files pertaining to a modified
  636.                 floating point package for
  637.                 BDS-C.  These files upgrade
  638.                 my previous FLOATXT.*
  639.                 submission.  Files are
  640.         FLOAT+44.C  Source, includes fixed "z" in _spr
  641.         FLOAT+44.CRL compiled version with V 1.44
  642.         FLOAT+44.DOC Documentation from FLOAT.DOC, updated
  643.         NEWFLVAL.C   Source of test program for float pkg
  644.         NEWFLVAL.CRL compiled version with V 1.44 & FLOAT+44
  645.         NEWFLVAL.COM linked version with FLOAT+44 for test
  646.  
  647.         COEFSTAT.C   Source of program to derive octal string
  648.                 equivalent of fp constants
  649.         COEFSTAT.CRL compiled version with V 1.44
  650.  
  651.         CTRIG.C      Source of updated trignometric functions
  652.