The tan function
The Fortran bli is in bli-tan.zip.
The starting and ending points for the tan could not be ±p/2 owing to the inablility of Fortran to produce the correct results. The following even required a slight modification to BLI to keep the upper part of the range exactly right. It was also necessary to change the **8 factor to **2 to properly see the singular regions at the end.
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION XI(5001),YI(5001),ERR(5001)
EXTERNAL BTAN
NP=5000
XI(1)=-ATAN(1D20)
XI(2)=ATAN(1D20)
PRINT*,XI(1),XI(2)
PRINT*,' TAN(XI(1))',TAN(XI(1))
PRINT*,'
TAN(XI(2))',TAN(XI(2))
NDO=2
CALL BLI(XI,YI,ERR,NP,NDO,BTAN)
OPEN(1,FILE='BLI.OUT')
WRITE(1,'(2G14.6)')(XI(I),YI(I),I=1,NP)
STOP
END
Bli modification
IF(NDO.EQ.2)THEN
NDO=NP/2
H=(XI(2)-XI(1))/(NDO-1)
XI(NDO)=XI(2) ßSo
that it would not be recalculated in the integer multiplication with h.
DO I=1,NDO-1
XI(I)=XI(1)+(I-1)*H
FI(I)=FEXT(XI(I))
ENDDO
FI(NDO)=FEXT(XI(NDO)) ß So
that the value would be as calculated.
ENDIF
Figure 1 Detail of tan(x) for x near p/2.
Near p/2, the sine is 1-(p/2-x)2, while the cosine is (p/2-x), thus the tan is 1/(p/2-x). The files below are in ttan.zip.
ttan.for
(4,1,36871) WATFOR-77 V3.0 PC/DOS
<beginning
of file>
IMPLICIT REAL*8 (A-H,O-Z)
DATA PI/3.141592653589793D0/
DO I=1,10
EPS=10D0**(-I-10)
PRINT*,EPS,TAN(PI/2-EPS)
ENDDO
END
<end of
file>
RUN
1.000000000000000D-011 9.999937915559000D+010
1.000000000000000D-012 9.998498645389550D+011
1.000000000000000D-013 1.000186738474650D+013
1.000000000000000D-014 9.947017564685660D+013
1.000000000000000D-015 8.536207345538900D+014
1.000000000000000D-016 1.632455227761910D+016
1.000000000000000D-017 1.632455227761910D+016
1.000000000000000D-018 1.632455227761910D+016
1.000000000000000D-019 1.632455227761910D+016
1.000000000000000D-020 1.632455227761910D+016
F:\Newton>WFL386
TTAN
Open Watcom
F77/32 Compile and Link Utility Version 1.3
wfc386 ttan.for
ttan.for: 7
statements, 101 bytes, 2 extensions, 0 warnings, 0 errors
creating a
Windows NT character-mode executable
F:\Newton>TTAN
9.9999999999999990D-012 9.9999379426637130D+010
1.0000000000000000D-012 9.9984989163587170D+011
1.0000000000000000D-013 1.0001870096265090D+013
1.0000000000000000D-014 9.9470443833547950D+013
1.0000000000000000D-015 8.5364048560630980D+014
1.0000000000000000D-016 1.6331778728383840D+016
1.0000000000000000D-017 1.6331778728383840D+016
1.0000000000000000D-018 1.6331778728383840D+016
1.0000000000000000D-019 1.6331778728383840D+016
9.9999999999999990D-021 1.6331778728383840D+016
F:\Newton\ttan\Debug>ttan
msdev-ttan.zip – compiled with Fortran power
station è that IEEE standard has this limit
1.000000000000001E-011 9.999937942663713E+010
1.000000000000001E-012 9.998498916358717E+011
1.000000000000001E-013 1.000187009626509E+013
1.000000000000001E-014 9.947044383354795E+013
1.000000000000002E-015 8.536404856063098E+014
1.000000000000002E-016 1.633177872838384E+016
1.000000000000002E-017 1.633177872838384E+016
1.000000000000002E-018 1.633177872838384E+016
1.000000000000002E-019 1.633177872838384E+016
1.000000000000002E-020 1.633177872838384E+016
cttan.wpj in ttan.zip
#include
<stdlib.h>
#include
<stdio.h>
#include
<math.h>
void main(){
double pi=3.141592653589793,eps;
int i;
for(i=1;i<11;i++){
eps=pow(10.,-i-10);
printf("%lg %lg
\n",eps,tan(pi/2-eps));}}
F:\Newton>cttan
1e-011
9.99994e+010
1e-012
9.9985e+011
1e-013
1.00019e+013
1e-014
9.94704e+013
1e-015
8.5364e+014
1e-016
1.63318e+016
1e-017
1.63318e+016
1e-018
1.63318e+016
1e-019
1.63318e+016
1e-020
1.63318e+016
Borland C++
#include
<stdlib.h>
#include
<stdio.h>
#include
<math.h>
void main(){
double pi=3.141592653589793,eps;
int i;
for(i=1;i<11;i++){
eps=pow(10.,-i-10);
printf("%lg %lg
\n",eps,tan(pi/2-eps));}}
F:\Newton>ttan
1e-11
9.99994e+10
1e-12
9.9985e+11
1e-13
1.00019e+13
1e-14
9.94702e+13
1e-15
8.53621e+14
1e-16
1.63246e+16
1e-17
1.63246e+16
1e-18
1.63246e+16
1e-19
1.63246e+16
1e-20
1.63246e+16
Conclusion
The tan is reliable to p/2 – 1e-11 with 5 digit accuracy and not reliable at all closer than 1e-20, but does produce a positive infinity rather than a negative one for positive epsilon. Note that there is an accuracy gain in going to the argument of the tan. That is the argument producing 9.99994e10 will be p/2 – 1e-11 to 16 digits.