Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- program SinCosTan
- !use iso_fortran_env, only: real64
- implicit none
- integer, parameter :: real64 = selected_real_kind(15,307)
- real(real64), parameter :: PI = 3.1415926535897932384626433832795028842
- real(real64), parameter :: TAU = 6.2831853071795864769252867665590057684
- real(real64), parameter :: half = 0.500000000000000000000_real64
- real(real64) :: time(2), times
- real(real64), allocatable:: ang(:), trigs(:,:)
- character(len=12) :: tout
- integer :: i,j,ierr,imax
- open(unit=10,file='trig.out',status='replace')
- open(unit=12,file='trig.in',status='old')
- ! Read input values from file
- i=0
- do
- read(12,*,iostat=ierr) ang
- if(ierr/=0) then
- imax = i
- exit
- end if
- i=i+1
- end do
- ! Now we know the number of input angles, allocate arrays
- allocate(ang(imax))
- allocate(trigs(imax,3))
- ! compute trig functions & time it
- times = 0.0_real64
- call cpu_time(time(1))
- do i=1,imax
- call findem(ang(i),trigs(i,:),40)
- end do
- call cpu_time(time(2))
- times = times + time(2) - time(1)
- ! Output to trig.out
- do i=1,imax
- do j=1,3
- if(trigs(i,j) > 1d100) then
- write(tout,'(a1)') 'n'
- elseif(abs(trigs(i,j)) > 1.0) then
- write(tout,'(f10.6)') trigs(i,j)
- elseif(abs(trigs(i,j)) < 0.1) then
- write(tout,'(f10.8)') trigs(i,j)
- else
- write(tout,'(f9.7)') trigs(i,j)
- end if
- write(10,'(a)',advance='no') tout
- end do
- write(10,*)" "
- end do
- ! Print timings and cleanup
- print *,"computation took",times/real(imax,real64)*1e6,"us"
- close(10); close(12)
- deallocate(ang);deallocate(trigs)
- contains
- !> @brief compute sine/cosine/tangent
- subroutine findem(a,t,iters)
- real(real64), intent(in) :: a
- real(real64), intent(inout) :: t(3)
- integer, intent(in) :: iters
- ! local variables
- real(real64), parameter :: deg2rad = 1.745329252e-2
- real(real64), parameter :: angles(60) = &
- [ 7.8539816339744830962e-01_real64, 4.6364760900080611621e-01_real64, &
- 2.4497866312686415417e-01_real64, 1.2435499454676143503e-01_real64, &
- 6.2418809995957348474e-02_real64, 3.1239833430268276254e-02_real64, &
- 1.5623728620476830803e-02_real64, 7.8123410601011112965e-03_real64, &
- 3.9062301319669718276e-03_real64, 1.9531225164788186851e-03_real64, &
- 9.7656218955931943040e-04_real64, 4.8828121119489827547e-04_real64, &
- 2.4414062014936176402e-04_real64, 1.2207031189367020424e-04_real64, &
- 6.1035156174208775022e-05_real64, 3.0517578115526096862e-05_real64, &
- 1.5258789061315762107e-05_real64, 7.6293945311019702634e-06_real64, &
- 3.8146972656064962829e-06_real64, 1.9073486328101870354e-06_real64, &
- 9.5367431640596087942e-07_real64, 4.7683715820308885993e-07_real64, &
- 2.3841857910155798249e-07_real64, 1.1920928955078068531e-07_real64, &
- 5.9604644775390554414e-08_real64, 2.9802322387695303677e-08_real64, &
- 1.4901161193847655147e-08_real64, 7.4505805969238279871e-09_real64, &
- 3.7252902984619140453e-09_real64, 1.8626451492309570291e-09_real64, &
- 9.3132257461547851536e-10_real64, 4.6566128730773925778e-10_real64, &
- 2.3283064365386962890e-10_real64, 1.1641532182693481445e-10_real64, &
- 5.8207660913467407226e-11_real64, 2.9103830456733703613e-11_real64, &
- 1.4551915228366851807e-11_real64, 7.2759576141834259033e-12_real64, &
- 3.6379788070917129517e-12_real64, 1.8189894035458564758e-12_real64, &
- 9.0949470177292823792e-13_real64, 4.5474735088646411896e-13_real64, &
- 2.2737367544323205948e-13_real64, 1.1368683772161602974e-13_real64, &
- 5.6843418860808014870e-14_real64, 2.8421709430404007435e-14_real64, &
- 1.4210854715202003717e-14_real64, 7.1054273576010018587e-15_real64, &
- 3.5527136788005009294e-15_real64, 1.7763568394002504647e-15_real64, &
- 8.8817841970012523234e-16_real64, 4.4408920985006261617e-16_real64, &
- 2.2204460492503130808e-16_real64, 1.1102230246251565404e-16_real64, &
- 5.5511151231257827021e-17_real64, 2.7755575615628913511e-17_real64, &
- 1.3877787807814456755e-17_real64, 6.9388939039072283776e-18_real64, &
- 3.4694469519536141888e-18_real64, 1.7347234759768070944e-18_real64]
- real(real64), parameter :: kvalues(33) = &
- [ 0.70710678118654752440e+00_real64, 0.63245553203367586640e+00_real64, &
- 0.61357199107789634961e+00_real64, 0.60883391251775242102e+00_real64, &
- 0.60764825625616820093e+00_real64, 0.60735177014129595905e+00_real64, &
- 0.60727764409352599905e+00_real64, 0.60725911229889273006e+00_real64, &
- 0.60725447933256232972e+00_real64, 0.60725332108987516334e+00_real64, &
- 0.60725303152913433540e+00_real64, 0.60725295913894481363e+00_real64, &
- 0.60725294104139716351e+00_real64, 0.60725293651701023413e+00_real64, &
- 0.60725293538591350073e+00_real64, 0.60725293510313931731e+00_real64, &
- 0.60725293503244577146e+00_real64, 0.60725293501477238499e+00_real64, &
- 0.60725293501035403837e+00_real64, 0.60725293500924945172e+00_real64, &
- 0.60725293500897330506e+00_real64, 0.60725293500890426839e+00_real64, &
- 0.60725293500888700922e+00_real64, 0.60725293500888269443e+00_real64, &
- 0.60725293500888161574e+00_real64, 0.60725293500888134606e+00_real64, &
- 0.60725293500888127864e+00_real64, 0.60725293500888126179e+00_real64, &
- 0.60725293500888125757e+00_real64, 0.60725293500888125652e+00_real64, &
- 0.60725293500888125626e+00_real64, 0.60725293500888125619e+00_real64, &
- 0.60725293500888125617e+00_real64 ]
- real(real64) :: beta, c, c2, factor, poweroftwo, s
- real(real64) :: s2, sigma, sign_factor, theta, angle
- integer :: j
- ! scale to radians
- beta = a*deg2rad
- ! ensure angle is shifted to appropriate range
- call angleShift(beta, -PI, theta)
- ! check for signs
- if( theta < -half*PI) then
- theta = theta + PI
- sign_factor = -1.0_real64
- else if( half*PI < theta) then
- theta = theta - PI
- sign_factor = -1.0_real64
- else
- sign_factor = +1.0_real64
- endif
- ! set up some initializations...
- c = 1.0_real64
- s = 0.0_real64
- poweroftwo = 1.0_real64
- angle = angles(1)
- ! run for some iterations
- do j=1,iters
- sigma = merge(-1.0_real64, +1.0_real64, theta < 0.0_real64)
- factor = sigma*poweroftwo
- c2 = c - factor*s
- s2 = factor*c + s
- c = c2
- s = s2
- ! update remaining angle
- theta = theta - sigma*angle
- poweroftwo = poweroftwo*0.5_real64
- if(j+1 > 60) then
- angle = angle * 0.5_real64
- else
- angle = angles(j+1)
- endif
- enddo !- j
- if(iters > 0) then
- c = c*Kvalues(min(iters,33))
- s = s*Kvalues(min(iters,33))
- endif
- c = c*sign_factor
- s = s*sign_factor
- t = [s, c, s/c]
- end subroutine findem
- subroutine angleShift(alpha, beta, gamma)
- real(real64), intent(in) :: alpha, beta
- real(real64), intent(out) :: gamma
- if(alpha < beta) then
- gamma = beta - mod(beta - alpha, TAU) + TAU
- else
- gamma = beta + mod(alpha - beta, TAU)
- endif
- end subroutine angleShift
- end program SinCosTan
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement