开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 操作系统开发 > 查看源码
EMATH.C
资源名称：ertos.rar　[点击查看]
上传用户：sunrenlu
上传日期：2022-06-13
资源大小：1419k
文件大小：11k
源码类别：

操作系统开发

开发平台：

DOS

EMATH.C：源码内容
  1. /*
  2.  * emath.c
  3.  *
  4.  * Copyright (c) 1987, 1999 Erick Engelke
  5.  */
  6. #include <string.h>
  7. #include <stdio.h>
  8. #include <stdlib.h>
  9. #include <ctype.h>
  11. /*------------------------------------------------------------------------*/
  12. #include <emath.h>
  14. /*------------------------------------------------------------------------*/
  15. EM EM_INT( EM x )
  16. {
  17.     int sign = 1;
  19.     if ( x < EM_ZERO ) {
  20.         sign = -1;
  21.         x = -x;
  22.     }
  23.     return( sign * (x & EM_INTEGERS) );
  24. }
  26. /*------------------------------------------------------------------------*/
  27. EM EM_ROUND( EM x )
  28. {
  29.     x += (x > 0) ? EM_POINT_5 : - EM_POINT_5 ;
  30.     return( x & EM_INTEGERS );
  31. }
  32. /*------------------------------------------------------------------------*/
  33. EM EM_SGN( EM x )
  34. {
  35.     if ( x < 0 ) return( EM_MINUS_ONE );
  36.     if ( x > 0 ) return( EM_ONE );
  37.     return( EM_ZERO );
  38. }
  39. /*------------------------------------------------------------------------*/
  40. EM EM_ABS( EM x )
  41. {
  42.     if ( x & EM_SIGN_BIT ) x *= -1 ;
  43.     return( x );
  44. }
  45. /*------------------------------------------------------------------------*/
  46. EM EM_MUL( EM a, EM b )
  47. {
  48.     EM ah, al, bh, bl, r;
  49.     long sign = 1;
  51.     if ( a < EM_ZERO ) {
  52.         sign = -1;
  53.         a = -a;
  54.     }
  55.     if ( b < EM_ZERO ) {
  56.         sign *= -1;
  57.         b = -b;
  58.     }
  61.     ah = a >> 16;
  62.     al = a & 65535;
  63.     bh = b >> 16;
  64.     bl = b & 65535;
  66.     r = (( ah * bh )  << 17)
  67.         + (  (  ( al * bh ) + (  ah * bl ) ) << 1 )
  68.       + ((( al>>1 ) * (bl>>1)) >> 13 );
  69.     return( r * sign );
  70. }
  71. /*------------------------------------------------------------------------*/
  72. EM EM_DIV( EM num , EM denom )
  73. {
  74.     int i;
  75.     EM r;
  76.     int sign = 1;
  78.     if ( num < 0 ) {
  79.         num = -num;
  80.         sign = -1;
  81.     }
  82.     if (denom <0) {
  83.         denom = -denom;
  84.         sign *= -1;
  85.     }
  87.     /* prescale to preserve accuracy */
  88.     i = EM_BDIGITS;
  90.     for ( i = EM_BDIGITS ; i ; i-- ) {
  91.         if ( ( num & EM_HIGH_BIT ) == 0 ) num *= 2; /* preserve sign */
  92.         else if ((denom & 1)==0) denom /= 2;
  93.         else break;
  94.     }
  96.     if ( denom == 0 ) return( EM_MAXINT );
  98.     r = num / denom;
  100.     /* finish scaling */
  101.     if ( i > 0 ) r <<= i;
  103.     return( r * sign );
  104. }
  105. /*------------------------------------------------------------------------*/
  106. /* Conversion Routines                                                    */
  107. /*------------------------------------------------------------------------*/
  108. /*------------------------------------------------------------------------*/
  109. EM EM_ASSIGN( long numerator, long denominator)
  110. {
  111.     EM scale = EM_ONE;
  113.     while ( !( scale & 1 ) && !( denominator &  1 )) {
  114.         scale >>= 1;
  115.         denominator >>= 1;
  116.     }
  117.     return( numerator * scale / denominator);
  118. }
  120. /*------------------------------------------------------------------------*/
  121. long emtol( EM em, long divisor )
  122. {
  123.     long basedivisor = EM_ONE;
  125.     /* try to minimize lost bits */
  126.     while ( ((basedivisor & 1 ) == 0 ) && ((divisor &1 )==0 ) ) {
  127.         basedivisor >>= 1;
  128.         divisor >>= 1;
  129.     }
  130.     return( (em * divisor) / basedivisor );
  131. }
  133. /*------------------------------------------------------------------------*/
  134. EM strtoem( char *s )
  135. {
  136.     EM x;
  137.     long j,k;
  138.     int i;
  139.     char *p;
  140.     int sign = 1;
  142.     /* skip leading spaces */
  143.     while ( isspace( *s )) s++;
  145.     /* fix sign, and more intermediate spaces */
  146.     if ( *s == '-' ) {
  147.         sign = -1;
  148.         s++;
  149.         while ( isspace( *s )) s++;
  150.     }
  152.     x = 0;
  153.     x = atol( s ) * EM_DENOMINATOR;
  154.     if ( ( p = strchr( s,'.' )) != NULL ) {
  155.         ++p;
  156.         for ( j = 1, i = 0 ; isdigit( p[i] ) ; ++i, j*=10 );
  157.         x += (atol( p ) * EM_DENOMINATOR + (j>>1)) / j;
  158.     }
  159.     if ( sign < 0 ) x = -x;
  160.     return( x );
  161. }
  162. /*------------------------------------------------------------------------*/
  163. char *emtostr( EM x, char *s )
  164. {
  165.     int sign = ( x < 0 ) ? -1 : 1;
  166.     char *p = s;
  167.     char ch;
  169.     if ( sign == -1  ) {
  170.         *p++ = '-';
  171.         x = -x;
  172.     }
  173.     sprintf( p, "%li.%0*li", x / EM_DENOMINATOR , EM_DEC_DIGITS,
  174.         (( x % EM_DENOMINATOR) * EM_DECIMAL_DENOMINATOR +
  175.             (EM_DECIMAL_DENOMINATOR>>1) )/ EM_DENOMINATOR );
  177.     /* kill trailing zeros */
  178.     p = strrchr( p, 0 );
  179.     while ( --p != s ) {
  180.         ch = *p;
  182.         /* don't need a trailing '.' */
  183.         if ( ch == '.' ) {
  184.             *p = 0;
  185.             break;
  186.         }
  187.         if ( ch == '0' ) *p = 0;
  188.         else break;
  189.     }
  190.     return( s );
  191. }
  193. /*------------------------------------------------------------------------*/
  194. /*------------------------------------------------------------------------*/
  195. /* Mathematic Routines                                                    */
  196. /*------------------------------------------------------------------------*/
  197. EM EM_SQR( EM in )
  198. {
  199.     EM newx, guess;
  200.     long scale;
  201.     int iteration;
  203.     if ( in <  EM_ZERO ) return( EM_MAXNEGINT );    /* invalid input */
  205.     /* reduce range */
  206.     for ( scale = 1 ; in > EM_ONE ; scale *= 10 )
  207.         in /= 100;
  209.     /* polynomial guess */
  210.     guess = EM_ADD( EM_ASSIGN(  1154,10000),
  211.             EM_MUL( EM_ASSIGN( 11544,10000), in ));
  213.     /* refine */
  214.     for ( iteration = 5 ; iteration ; --iteration ) {
  215.         guess = EM_MUL( EM_POINT_5, EM_ADD(  guess, EM_DIV( in , guess )));
  216.     }
  218.     return( guess * scale );
  219. }
  220. /*------------------------------------------------------------------------*/
  221. EM EM_LOG10( EM in )
  222. {
  223.     long sign = 1;
  224.     EM scale, x3, x4, y;
  225.     EM s, z, r;
  227.     if ( in <= EM_ZERO ) return( EM_MAXNEGINT );    /* invalid input */
  229.     if ( in < EM_ONE ) {
  230.         sign = -1;
  231.         in = EM_DIV( EM_ONE, in );
  232.     }
  234.     /* scale to useful */
  235.     scale = EM_ZERO;
  236.     while ( in >= EM_ONE ) {
  237.         in /= 10 ;
  238.         scale = EM_ADD( scale, EM_ONE );
  239.     }
  240.     in = EM_MUL( EM_ASSIGN( 31623, 10000 ), in );
  241.     in = EM_DIV( EM_SUB( in, EM_ONE ), EM_ADD( in, EM_ONE ));
  242.     y = EM_MUL( in, in );
  243.     x3 = EM_ADD(
  244.             EM_ASSIGN( 31878, 10000 ),
  245.             EM_MUL( y,
  246.                     EM_ADD(
  247.                         EM_ASSIGN( -26558,10000),
  248.                         EM_MUL( y, EM_ASSIGN(2669, 10000))
  249.                     )
  250.                  ) );
  252.     x4 = EM_ADD( EM_ASSIGN( 36701, 10000),
  253.                  EM_MUL( y,
  254.                     EM_ADD( EM_ASSIGN( -42810,10000), y) ));
  255.     y = EM_SUB(
  256.             EM_ADD( scale,
  257.                     EM_DIV(
  258.                         EM_MUL(in,x3),
  259.                         x4)),
  260.             EM_POINT_5) *sign;
  261.     return( y );
  262. }
  263. /*------------------------------------------------------------------------*/
  264. EM EM_POWER10( EM x )
  265. {
  266.     int inv = 0;
  267.     EM n, y, x3, x4;
  269.     if ( x < EM_ZERO ) {
  270.         inv = 1;
  271.         x = -x;
  272.     }
  273.     n = EM_INT( x );
  274.     x = EM_SUB( EM_SUB( x, n ), EM_POINT_5 );
  275.     y = EM_MUL( x, x );
  276.     x3 = EM_ADD( EM_ASSIGN( 695570,10000),
  277.                  EM_MUL( y, EM_ASSIGN( 87518, 10000)));
  278.     x4 = EM_ADD( EM_ASSIGN( 604164, 10000),
  279.                  EM_MUL( y,
  280.                     EM_ADD( EM_ASSIGN( 342951, 10000), y )));
  281. //    y = 2 * x3 * x  / ( x4 - x * x3 );
  283.     y = EM_DIV( 2 * EM_MUL( x3, x ),
  284.                 EM_SUB( x4,
  285.                         EM_MUL( x, x3 )));
  287.     y = EM_MUL( EM_ASSIGN(31623,10000) ,
  288.                 EM_ADD( EM_ONE, y ));
  290.     while ( n > 0 ) {
  291.         y *= 10;
  292.         n = EM_SUB( n, EM_ONE );
  293.     }
  294.     if ( inv ) y = EM_DIV( 1, y );
  295.     return( y );
  296. }
  297. /*------------------------------------------------------------------------*/
  298. EM EM_LN( EM x )
  299. {
  300. #define EM_LOG10_E 14231
  301.     return( EM_DIV( EM_LOG10( x ) , EM_LOG10_E ));
  302. }
  303. /*------------------------------------------------------------------------*/
  304. EM EM_EXP( EM x )
  305. {
  306.     return( EM_POWER10( EM_MUL( x , EM_LOG10_E )));
  307. }
  308. /*------------------------------------------------------------------------*/
  309. EM EM_POWER( EM x, EM y )
  310. {
  311.     return( EM_POWER10( EM_MUL( y, EM_LOG10( x ))));
  312. }
  313. /*------------------------------------------------------------------------*/
  314. EM EM_COS( EM rad )
  315. {
  316.     EM y;
  317.     int test;
  318.     int sign = 1;
  320.     if ( rad < 0 )
  321.         rad = -rad;
  323.     /* get into scale */
  324.     rad = EM_SUB( rad,
  325.                   EM_MUL( EM_INT( EM_DIV( rad, EM_TWO_PI )),
  326.                           EM_TWO_PI ));
  327.     /* stage of rad in range 0..2pi */
  328.     test = emtol( EM_DIV( rad, EM_90_DEG_RAD ), 1 );
  330.     switch (test) {
  331.         case 0 : break;
  332.         case 1 : rad = EM_ABS( EM_SUB( EM_PI, rad )) ; sign = -1 ; break;
  333.         case 2 : rad = EM_ABS( EM_SUB( EM_PI, rad )) ; sign = -1 ; break;
  334.         case 3 : rad = EM_ABS( EM_SUB( rad, EM_TWO_PI)); break;
  335.     }
  337.     rad = EM_MUL( rad, rad );
  338.     y = EM_ADD( EM_ONE,
  339.                 EM_MUL( rad,
  340.                         EM_ADD( EM_ASSIGN( -5, 10 ),
  341.                                 EM_MUL( rad,
  342.                                         EM_ADD( EM_ASSIGN( 417, 10000),
  343.                                                 EM_MUL( rad,
  344.                                                         EM_ASSIGN( -14, 10000)
  345.                                                        )
  346.                                         )
  347.                                  )
  348.                          )
  349.                  )
  350.                 );
  351.     return( y * sign );
  352. }
  354. /*------------------------------------------------------------------------*/
  355. EM EM_SIN( EM rad )
  356. {
  357.     return( EM_COS( rad - EM_PI/2 ));
  358. }
  359. /*------------------------------------------------------------------------*/
  360. EM EM_TAN( EM rad )
  361. {
  362.     EM x, y;
  364.     x = EM_SIN( rad );
  365.     y = EM_COS( rad );
  367.     return( EM_DIV( x, y ));
  368. }
  369. /*------------------------------------------------------------------------*/
  370. EM EM_ATAN( EM x )
  371. {
  372.     long sign = 1;
  373.     int flag = 0;
  374.     EM x3, x4, y;
  376.     if ( x < EM_ZERO ) {
  377.         sign = -1;
  378.         x = -x;
  379.     }
  380.     if ( x > EM_ONE ) {
  381.         x = EM_DIV( EM_ONE, x );
  382.         flag = 1;
  383.     }
  384.     y = EM_MUL( x , x );
  385.     x3 = EM_ADD( EM_ASSIGN( 8079, 100 ),
  386.                  EM_MUL( y,
  387.                          EM_ADD( EM_ASSIGN( 725814, 10000 ),
  388.                                  EM_MUL( y,
  389.                                          EM_ASSIGN( 111177, 10000 )))));
  390.     x4 = EM_ADD( EM_ASSIGN( 8079, 100 ),
  391.                  EM_MUL( y,
  392.                          EM_ADD( EM_ASSIGN( 995113, 10000 ),
  393.                                  EM_MUL( y,
  394.                                          EM_ADD( EM_ASSIGN( 281329, 10000),
  395.                                                  y )))));
  397.     y = EM_DIV( EM_MUL( x, x3) , x4 );
  398.     if ( flag )
  399.         y = EM_SUB( EM_ASSIGN(15708,10000), y );
  400.     return( y * sign );
  401. }
  403. /*------------------------------------------------------------------------*/
  404. EM EM_ASIN( EM x )
  405. {
  406.     EM c;
  408.     c = EM_DIV( x, EM_SQR( EM_SUB( EM_ONE, EM_MUL( x, x ))));
  409.     return( EM_ATAN( c ));
  410. }
  411. /*------------------------------------------------------------------------*/
  412. EM EM_ACOS( EM x )
  413. {
  414.     EM c;
  415.     c = EM_DIV( EM_SQR( EM_SUB( EM_ONE,
  416.                                 EM_MUL( x, x ))),
  417.                 x );
  418.     return( EM_ATAN( c ));
  419. }
  420. /*------------------------------------------------------------------------*/
  421. EM EM_DEG( EM rad )
  422. {
  423.     return( EM_MUL( rad , EM_ASSIGN( 572958, 10000 )));
  424. }
  425. /*------------------------------------------------------------------------*/
  426. EM EM_RAD( EM deg )
  427. {
  428.     return( EM_DIV( deg , EM_ASSIGN( 572958, 10000 )));
  429. }
  430.