home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_100 / 188_01 / trns1.c < prev    next >
Text File  |  1987-09-30  |  11KB  |  314 lines
  1. /*
  2. HEADER:        ;
  3. TITLE:        C elementary transcendentals;
  4. VERSION:    1.1;
  5.  
  6. DESCRIPTION║ááá    Sourcσááá codσááá fo≥áá al∞áá standarΣááá ├ ì
  7.                transcendentals; revision of CUG 187
  8.  
  9.         Employ≤á ldexp(⌐á anΣá frexp(⌐á functions╗á iµ ì
  10.                suitablσá version≤á oµ thesσ arσá no⌠á provideΣ ì
  11.                b∙á ß giveε compiler¼á thσ version≤ provideΣ iε ì
  12.                sourcσá codσá wil∞ requirσá adaptatioεá t∩á thσ ì
  13.                doublσ floa⌠ format≤ oµ thσ compiler.
  14.  
  15. KEYWORDS:    Transcendentals, library;
  16. SYSTEM:        CP/M-80, V3.1;
  17. FILENAME:    TRNS1.C;
  18. WARNINGS║áá    frexp(⌐áá anΣá ldexp(⌐á arσáá implementatioε ì
  19.                dependent«á  Thσá compile≥á employeΣá doe≤á no⌠ ì
  20.                suppor⌠ááá minu≤áá (-⌐áá unar∙áá operator≤áá iε ì
  21.                initialize≥á lists¼á whicΦ arσ requireΣ b∙á thσ ì
  22.                code.
  23. AUTHORS:    Tim Prince;
  24. COMPILERS:    MIX C, v2.0.1; Eco C 3.45; Aztec C 1.06B
  25. */
  26. /* Copyright (C) 1986, Chandler Software,
  27. ** 4456 W. Maple Rd., Birmingham MI 48010-1923
  28. ** (313) 855-4587
  29. ** prepared for non-profit distribution by C Users' Group
  30. ** revised again 17.V.86 and 12.X.86
  31. ** all functions work in MIX Z80 C except sin(), atan()
  32. ** To fix sin() and atan(), CONVERT the .MIX file to ASCII,
  33. ** replace the minus signs on the table[] constants, adding 1
  34. ** to the leading hex numbers which give the number of ASCII
  35. ** characters in each constant.
  36. ** Eco C requires the negative bit to be set in the DEFW
  37. ** fields of the .MAC file.  See the normal source version for
  38. ** the initializer list constants which are negative.
  39. ** All functions work in Aztec C except frexp() and ldexp(),
  40. ** which are already present and working in m.lib.
  41. ** Eco-C version, uses Eco-C constants and dint() instead of
  42. ** floor() */
  43. #define ODD(i) ((i)&1) /* also like Pascal */
  44. #define ECOC 1 /* Ecosoft has some special definitions in S0 */
  45. #ifdef ECOC
  46. #define HALF _inv2
  47. #define SR5 _sqrt5
  48. #define INFINITY _FPMAX /*IEEE 0x7ff0000000000000*/
  49. #define FRND(x,y) dint(x+(y>=0)-HALF) /* no floor() */
  50. #define PI _pi
  51. #define HPI _pi_2
  52. #define ONE _FPONE
  53. extern double PI,HPI,HALF,SR5,INFINITY,ONE,dint();
  54. #else
  55. #define HALF .5
  56. #define SR5 .7071
  57. #define ONE 1
  58. #define INFINITY 1.7e38
  59. #define FRND(x,y) floor(x+.5)
  60. extern double floor();
  61. #define PI 3.1415926535897932385
  62. #define HPI 1.5707963267948966192
  63. #endif
  64. #define ROUND(x) (int)((x>=0)-HALF+x) /* Pascal ROUND */
  65. /* teach the preprocessor some algebra */
  66. #ifdef FULLDEFS
  67. #define cos(x) sin(HPI-(x)) 
  68. #define fabs(x) (x>=0?x:-(x))
  69. #define MAXLN 88 /* log(overflow theshold) */
  70. #define LN2 .69314718055994530942
  71. #define L2E 1.442695040888963407
  72. #define atan2(x,y) (y>0?atan((x)/(y)):y==0?x>=0?HPI:-HPI\
  73.   :x>=0?PI+atan((x)/(y)):atan((x)/(y))-PI)
  74. #define acos(x) (HPI-asin(x))
  75. #define log(x) log2(x)*LN2 /* double log2() */
  76. #define exp(x) exp2((x)*L2E) /* double exp2() */
  77. #define cosh(x) sqrt(sinh(x)*sinh(x)+ONE)/* x>=20?fabs(sinh(x)) */
  78. #define tanh(x) (x<-20?-1:x>=20?ONE:sinh(x)/cosh(x))
  79. #endif
  80. #define pow(x,y) exp2(log2(x)*(y))
  81. main(){ /* test accuracy of tan,sin,cos,atan,exp2,log2 */
  82.   double tan(),log2(),exp2(),asin();
  83.   double x,ang;
  84.   for(x=1e-36;x<1e33;x*=1e4){
  85.     ang=atan(x);
  86.     printf("x:%23.16e tan(ang):%23.16e\n ang:%23.16e\n asin(sin(ang)):%23.16e\n pow:%23.16e\n",
  87.       x,tan(ang),ang,asin(sin(ang)),pow(x,ONE));
  88.   }
  89. }
  90. /* dblfmt describes floating point format
  91. ** machine dependencies in case frexp() and ldexp() are not
  92. ** supplied (Aztec has their own version, which works) */
  93. struct dblfmt{
  94. /* MIX C like Microsoft format */
  95. #ifndef ECOC
  96.   char mant[7]; /* full reverse byte order */
  97. #endif
  98. /* IEEE int sign:1;unsigned expn:11;unsigned mant:4
  99. ** with DEC reverse byte pairing, can't cross 16-bit
  100. ** boundaries */
  101.   char expn;
  102. #ifdef ECOC
  103.   char mant[7]; /* forward byte order: Ecosoft */
  104. #endif
  105. };
  106. double frexp(x,scale)
  107.   int *scale;
  108.   double x;
  109. {
  110.   #define BIAS 128 /* IEEE 1023 */
  111.   *scale=((struct dblfmt *)&x)->expn-BIAS;
  112.   ((struct dblfmt *)&x)->expn=BIAS;
  113.   return(x);
  114. }
  115. double ldexp(x,scale)
  116.   double x;
  117.   char scale;
  118. {
  119.   ((struct dblfmt *)&x)->expn+=scale; /* return x*2**scale */
  120.   return(x);
  121. }
  122. /* use fast polynomial evaluation (possibly non-portable) 
  123. ** loop unrolling or odd & even summing or such strategies
  124. ** may be faster on fast floating point machines */
  125. double poly(x,n,table)
  126.   double x;
  127.   char n;
  128.   double *table;
  129. {
  130.   register double sum;
  131.   sum=*table++;
  132.   while(--n)
  133.     sum=sum*x+*table++;
  134.   return(sum);
  135. }
  136. double asin(x) double x;
  137. { double sqrt();
  138.   return(x>=ONE?HPI:(x<=-1?-HPI:atan(x/sqrt(ONE-x*x))));
  139. }
  140. double tan(x) double x;
  141. { double xs;
  142.   x-=FRND(x/PI,x)*PI; /* if your hardware prefers, change sign*/
  143.   xs=x*x;
  144. /* minimax relative error fit obtained by Ruzinski program 
  145. ** if sign was changed above, change here to get correct result*/
  146.   return(x*(33281902.4774370361+xs*(xs*(
  147.     131095.960756651362+xs*(xs-968.863630016633889))
  148.     -4572612.25618456766))/
  149.     (33281902.4774370360+xs*(xs*(915702.213319541242+xs*(
  150.     xs*44.4083430808097903-13491.8003635866138))
  151.     -15666579.7486635766)));
  152. }
  153. double log2(x)  double x;
  154. {
  155.   #define NTL 7
  156. /* Chebyshev fit to log2(x)*(1+x)/(1-x) using Taylor series
  157. ** expansion and then reverting to economized polynomial 
  158. ** this series has 2% less absolute error than a minimax
  159. ** relative error fit, without increasing relative error */
  160.   static double table[]={
  161.     .2429264,.261443335,.3206164184,
  162.     .412198401587,.5770780172495,.961796693924327,
  163.     2.8853900817779273};
  164.   int scale,incsc;
  165.   double poly();
  166. /* split x -> 2^scale*(reduced x near 1) */
  167.   incsc=frexp(x,&scale)<SR5;
  168.   x=ldexp(x,-(scale-=incsc));
  169. /* (x-.5)-.5 recommended for inaccurate arithmetic but it may
  170. **  not help; then you might as well use x=frexp(x,&scale)
  171. **  and x=(x-sqrt(.5))/(x+sqrt(.5)), poly()*x-.5+scale
  172. ** with sqrt(.5) reduced for best results (reduce by 2 ULPs
  173. ** from correctly rounded value on Harris Hnnn) */
  174.   x=(x-ONE)/(x+ONE);
  175.   return(poly(x*x,NTL,table)*x+scale);
  176. }
  177. double exp2(x)  double x;
  178. {
  179.   #define MAXEXP 127 /* IEEE 1023 */
  180.   #define MINLG2 -128
  181.   double xsq;
  182.   int scale;
  183.   if(x>=MAXEXP)return(INFINITY);
  184.   if(x<MINLG2)return(0);
  185. /* 2^x=2^ROUND(x)*2^(x-ROUND(x)) */
  186.   x-=scale=ROUND(x);
  187. /* approximation is ratio of polynomials differing only in
  188. ** sign of odd terms */
  189.   xsq=x*x;
  190. /* continued fraction approx for e^x
  191.   (x*(15120+xsq*(420+xsq))+30240+xsq*(3360+xsq*30))/... */
  192. /* minimax absolute error fit for 2^x, 18 significant digits */
  193.   x*=65556.325877407443+xsq*(874.804294426022+xsq);
  194.   xsq=189155.572484473087+xsq*(10097.515376265486+
  195.     xsq*43.302654542155);
  196.   return(ldexp((x+xsq)/(xsq-x),scale));
  197. }
  198. #ifdef FULLDEFS
  199. double sinh(x)  double x;
  200. {
  201.   #define NTSH 7
  202.   static double table[NTSH]={
  203. /* Chebyshev fit to sinh(x)/x has 2% less absolute error than
  204. ** minimax relative error fit */
  205.     1.632881e-10,2.50483893e-8,2.7557344615e-6,
  206.     1.984126975233e-4,.0083333333334816,.1666666666666574,
  207.     1.0000000000000001};
  208.   double poly(),exp2() /*,absx*/;
  209.   if(( /*absx=*/ fabs(x))<ONE)return(poly(x*x,NTSH,table)*x);
  210. /*  if(absx<MAXLN){ */
  211.     x=exp(x);
  212.     return(ldexp(x-ONE/x,-1));/*}
  213.   return(x>0?exp(x-LN2):-exp(LN2-x));*/
  214. /* good enough for unusual case */
  215. }
  216. #endif
  217. double sqrt(x) /* as if division hardware but no sqrt */
  218.   double x;
  219. { static float table[]={.59,.417,.295}; /* minimax 2 digits */
  220.   register double x1;
  221.   register int i;
  222.   register char iters=3; /* number of Newton iterations */
  223.   int scale;
  224.   if(x<=0)return(x); /* set range error ? */
  225.   x1=frexp(x,&scale); /* sqrt(2^x*y)=2^(x/2)sqrt(y) */
  226.   i=ODD(scale); /* separate fits for y<.5 & y>.5 */
  227.   x1=ldexp(x1*table[i]+table[i+1],(scale+1)>>1); /*crude sqrt
  228. ** if there is a crude firmware sqrt use that instead, then
  229. ** iters may be reduced */
  230.   do /* enough Newton iterations for full accuracy */
  231.     x1=ldexp(x/x1+x1,-1); /* ldexp might be faster than *.5 
  232. ** if division doesn't round up, add .5 ULP between last
  233. ** division and addition (tested on Honeywell DPS8) */
  234.   while(--iters);
  235.   return(x1);
  236. }
  237. /* #define ECOC 1 */
  238. #ifdef ECOC
  239. #definσ TMR│ _tsqr3
  240. extern double TMR3;
  241. #else
  242. #define TMR3 .268 
  243. #endif
  244. double sin(x)  double x;
  245. {
  246. #ifndef ECOC
  247.   #define NTS 8
  248. /* negative signs are ignored by some compilers! */
  249.   static double table[NTS]={
  250. /* put minus signs in output file if compiler can't 
  251. ** Chebyshev fit to sin(x)/x gives relative error from 25%
  252. ** less (in middle) to 20% more (at ends of interval) than
  253. ** minimax relative error fit, but coefficients are closer to
  254. ** Taylor series */
  255.     -.7374418e-12,.160481709e-9,-.25051882036e-7,
  256.     .2755731661057e-5,-.00019841269824875,
  257.     .00833333333328281,-.1666666666666607,.9999999999999999};
  258. #else
  259.   #define NTS 9
  260. /* note omitted - signs ! */
  261.   static double table[NTS]={.27215822e-14,
  262.     .7643027274e-12,.160589409072e-9,.25052106891169e-7,
  263.     .2755731921123654e-5,.00019841269841208719,
  264.     .00833333333333318652,.1666666666666666531,
  265.     .9999999999999999998};
  266. #endif
  267. /* the corresponding minimax relative error coefficients are:
  268. ** -.7370812e-12,.160478576e-9,-.25051871327e-7,
  269. ** .2755731642894e-5,-.00019841269823293,.00833333333327625,
  270. ** -.1666666666666597,.9999999999999999 (16 digits)
  271. ** .27205187e-14,-.7642921847e-12,.160589366594e-9,
  272. ** -.25052106802022e-7,.2755731921020009e-5,
  273. ** -.00019841269841202181,.00833333333333316634e-2,
  274. ** -.166666666666666507,.9999999999999999997 (18.5 digits) */
  275.   double pm;
  276. /* try to take care of x>maxint unless much too slow
  277. ** employ periodicity sin(x+n*PI)=(-1)^n*sin(x) */
  278.   x-=PI*(pm=FRND(x/PI,x)); /* here it may be better to change
  279. ** sign; now |x| < HPI; use series and fix sign */
  280.   return(poly(x*x,NTS,table)*(ODD((int)pm)?-x:x));
  281. }
  282. double atan(x) double x;
  283. {
  284.   #define TNPI6 .57735026918962576451
  285. /* on Harris Hnnn, better results are obtained by reducing
  286. ** TNPI6 by 1 ULP to compensate for inaccurate arithmetic */
  287.   #define NTA 9
  288.   static double table[]={
  289. /* if compiler doesn't support initializers, must edit output*/
  290. #ifndef ECOC
  291. /* Chebyshev fit to atan(x)/x gives 2% less to .4% more error
  292. ** than minimax relative error fit */
  293.    .04438671796,-.064831131188,.07679359305498,
  294.     -.090903705713302,.111110978788314,-.1428571410247392,
  295.     .1999999999872629,-.33333333333329920,
  296. #else /* signs omitted to satisfy Eco C compiler */
  297.     .04438671796,.064831131188,.07679359305498,
  298.     .090903705713302,.111110978788314,.1428571410247392,
  299.     .1999999999872629,.33333333333329920,
  300. #endif
  301.     .99999999999999998};
  302.   char invert,reduce,neg;
  303. /* atan(-x)=-atan(x) */
  304.   if(neg=(x<0))x=-x;
  305. /* tan(HPI-x)=1/tan(x) */
  306.   if(invert=(x>=ONE))x=ONE/x;
  307. /* tan(a-b)={tan(a)-tan(b)}/{tan(a)*tan(b)+1}; b=PI/6 */
  308.   if(reduce=(x>=TMR3))x=(x-TNPI6)/(x*TNPI6+ONE);
  309.   x*=poly(x*x,NTA,table);
  310.   if(reduce)x+=.52359877559829887308;
  311.   if(invert)x=HPI-x;
  312.   return(neg?-x:x);
  313. }
  314.