home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Guide / c-cplusplus-interactive-guide.iso / c_ref / csource3 / 188_01 / translib.for < prev    next >
Encoding:
Text File  |  1987-10-01  |  3.9 KB  |  161 lines
  1. /*
  2. HEADER:        ;
  3. TITLE:        C elementary transcendentals;
  4. VERSION:    1.0;
  5.  
  6. DESCRIPTION:       Source    code    for   all   standard    C                            
  7.  
  8.                transcendentals
  9.  
  10.         Employs  ldexp()  and  frexp()  functions;  if                
  11.  
  12.                suitable  versions  of these are  not  provided                
  13.  
  14.                by  a given compiler,  the versions provided in                
  15.  
  16.                source  code  will require  adaptation  to  the                
  17.  
  18.                double float formats of the compiler.
  19.  
  20. KEYWORDS:    Transcendentals, library;
  21. SYSTEM:        CP/M-80, V3.1;
  22. FILENAME:    TRANS.C;
  23. WARNINGS:      frexp()   and  ldexp()  are   implementation                
  24.  
  25.                dependent.   The  compiler  employed  does  not 
  26.  
  27.                support    minus   (-)   unary   operators   in                
  28.  
  29.                initializer  lists,  which are required by  the                
  30.  
  31.                code.
  32. AUTHORS:    Tim Prince;
  33. COMPILERS:    MIX C, v2.0.1;
  34. */
  35.     program transt
  36.     real log
  37.     sum=0
  38.     do 1 i=1,1000
  39.       x=i*.01
  40. 1    sum=abs(tan(atan(x))/x-1)+sum
  41.     write(3,2)sum
  42. 2    format(' built-in tan(atan) error:',g15.7)
  43.     sum=0
  44.     do 3 i=1,1000
  45.       x=i*.01
  46. 3    sum=abs(exp(log(x))/x-1)+sum
  47.     write(3,4)sum
  48. 4    format(' built-in exp(log) error:',g15.7)
  49.     sum=0
  50.     do 5 i=1,1000
  51.       x=i*.01
  52. 5    sum=abs(xtan(xatan(x))/x-1)+sum
  53.     write(3,6)sum
  54. 6    format(' tan(atan) error:',g15.7)
  55.     sum=0
  56.     do 7 i=1,1000
  57.       x=i*.01
  58. 7    sum=abs(xexp(xlog(x))/x-1)+sum
  59.     write(3,8)sum
  60. 8    format(' exp(log) error:',g15.7)
  61.     end
  62.  
  63.       function xtan(x)
  64.     xtan=xsin(x)/xcos(x)
  65.     return
  66.       end
  67.  
  68.       function xcos(x)
  69.     xcos=xsin(1.57079633-x)
  70.     return
  71.       end
  72.  
  73.       function xsin(x)
  74.     real table(5)
  75.     data pi/3.141592653589793238/,
  76.      &      table/2.60198e-6,-.0001980746,.0083330258,
  77.      &        -.16666657,.99999999/
  78. c series to 28 bit accuracy
  79. c nearest multiple of pi
  80.     pim=anint(x/pi)
  81. c get remainder
  82.     xr=x-pi*pim
  83. c sin(x-pim*pi) = (-1)**pim * sin(x)
  84.     hpim=scal2(pim,-1)
  85. c if pim/2 not integer, pim is odd so change sign
  86.     if(aint(hpim).ne.hpim)xr=-xr
  87. c evaluate Horner polynomial, odd terms
  88.     xsin=poly(xr*xr,5,table)*xr
  89.     return
  90.       end
  91.  
  92.       function scal2(x,i)
  93. c    scal2=x*2**i
  94.     integer*1 ix(4),i
  95.     equivalence(xi,ix(1))
  96.     xi=x
  97.     ix(4)=i+ix(4)
  98.     scal2=xi
  99.     return
  100.       end
  101.  
  102.       function poly(x,n,table)
  103.     real table(n)
  104.     poly=table(1)
  105.     do 1 i=2,n
  106. 1    poly=poly*x+table(i)
  107.     return
  108.       end
  109.  
  110.       function xatan(x)
  111.     logical invert
  112.     real table(8)
  113.     data table/-.0046930,.0242506,-.0594827,.0991394,
  114.      &      -.1401932,.1996969,-.3333199,.9999999/
  115.     xr=abs(x)
  116.     invert=xr.gt.1
  117.     if(invert)xr=1/xr
  118.     xr=poly(xr*xr,8,table)*xr
  119. c tan(1/x)= tan(pi/2-x)
  120.     if(invert)xr=1.57079633-xr
  121.     xatan=sign(xr,x)
  122.     return
  123.       end
  124.  
  125.       function xlog(x)
  126.     xlog=alog2(x)*.693147181
  127.     return
  128.       end
  129.  
  130.       function alog2(x)
  131.     real table(3)
  132.     integer*1 ix(4),iexp
  133.     equivalence(xr,ix(1))
  134. c polynomial to 24 bit precision
  135.     data table/.5957779,.9615884,2.8853901/
  136.     xr=x
  137. c z'35' is leading 7 bits (minus hidden bit) of sqrt(.5)
  138. c shift so xr close to 1.
  139. c .true. value -1 on this system, bias z'80'
  140.     iexp=(ix(3).le.z'35')+ix(4)-z'80'
  141. c subtract 1 first so no accuracy lost
  142.     xr=scal2(xr,-iexp)-1
  143.     xr=xr/(xr+2)
  144. c polynomial in (xr-1)/(xr+1), .7 < xr < 1.4
  145.     alog2=poly(xr*xr,3,table)*xr+iexp
  146.     return
  147.       end
  148.  
  149.       function xexp(x)
  150.     xexp=exp2(x*1.442695040888963407)
  151.     return
  152.       end
  153.  
  154.       function exp2(x)
  155.     integer*1 expn
  156.     real table(7)
  157. c 26 bits precision -.5 < xr < .5 for 2**x
  158. c sinh (|x|<.7) use odd terms (x*1.442695)
  159.     data table/.000154,.00133904,.00961814,.05550358,
  160.      &      .2402265,.69314718,1./
  161.     if(x.ge.-128)go to 1
  162.       exp2=0
  163.       return
  164. 1    xr=amin1(x,126.9999)
  165.     expn=anint(xr)
  166.     xr=xr-expn
  167.     exp2=scal2(poly(xr,7,table),expn)
  168.     return
  169.       end
  170.