1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
|
program sin_02
use, intrinsic :: iso_fortran_env, only: dp => real64
implicit none
real(dp), parameter :: pi = 3.1415926535897932384626433832795_dp
real(dp) :: x, r1, r2, err
x = 10.5_dp
r1 = dsin(x)
r2 = -0.87969575997167_dp
err = abs(r1-r2)
print *, x, r1, r2
print *, err
if (err > 1e-15_dp) error stop
contains
real(dp) function abs(x) result(r)
real(dp), intent(in) :: x
if (x >= 0) then
r = x
else
r = -x
end if
end function
elemental integer function floor(x) result(r)
real(dp), intent(in) :: x
if (x >= 0) then
r = x
else
r = x-1
end if
end function
elemental real(dp) function modulo(x, y) result(r)
real(dp), intent(in) :: x, y
r = x-floor(x/y)*y
end function
elemental real(dp) function min(x, y) result(r)
real(dp), intent(in) :: x, y
if (x < y) then
r = x
else
r = y
end if
end function
elemental real(dp) function max(x, y) result(r)
real(dp), intent(in) :: x, y
if (x > y) then
r = x
else
r = y
end if
end function
elemental real(dp) function dsin(x) result(r)
real(dp), intent(in) :: x
real(dp) :: y
integer :: n
y = modulo(x, 2*pi)
y = min(y, pi - y)
y = max(y, -pi - y)
y = min(y, pi - y)
r = kernel_dsin(y)
end function
! Accurate on [-pi/2,pi/2] to about 1e-16
elemental real(dp) function kernel_dsin(x) result(res)
real(dp), intent(in) :: x
real(dp), parameter :: S1 = 0.9999999999999990771_dp
real(dp), parameter :: S2 = -0.16666666666664811048_dp
real(dp), parameter :: S3 = 8.333333333226519387e-3_dp
real(dp), parameter :: S4 = -1.9841269813888534497e-4_dp
real(dp), parameter :: S5 = 2.7557315514280769795e-6_dp
real(dp), parameter :: S6 = -2.5051823583393710429e-8_dp
real(dp), parameter :: S7 = 1.6046585911173017112e-10_dp
real(dp), parameter :: S8 = -7.3572396558796051923e-13_dp
real(dp) :: z
z = x*x
res = x * (S1+z*(S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*S8)))))))
end function
end program
|
