Advertisement
Guest User

Untitled

a guest
Jun 12th, 2014
395
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. program SinCosTan
  2.    !use iso_fortran_env, only: real64
  3.    implicit none
  4.    integer, parameter :: real64 = selected_real_kind(15,307)
  5.    real(real64), parameter :: PI  = 3.1415926535897932384626433832795028842
  6.    real(real64), parameter :: TAU = 6.2831853071795864769252867665590057684
  7.    real(real64), parameter :: half = 0.500000000000000000000_real64
  8.    real(real64) :: time(2), times
  9.    real(real64), allocatable:: ang(:), trigs(:,:)
  10.    character(len=12) :: tout
  11.    integer :: i,j,ierr,imax
  12.  
  13.    open(unit=10,file='trig.out',status='replace')
  14.    open(unit=12,file='trig.in',status='old')
  15.  
  16.    ! Read input values from file
  17.    i=0
  18.    do
  19.       read(12,*,iostat=ierr) ang
  20.       if(ierr/=0) then
  21.         imax = i
  22.         exit
  23.       end if
  24.       i=i+1
  25.    end do
  26.    
  27.    ! Now we know the number of input angles, allocate arrays
  28.    allocate(ang(imax))
  29.    allocate(trigs(imax,3))
  30.  
  31.    ! compute trig functions & time it
  32.    times = 0.0_real64
  33.    call cpu_time(time(1))
  34.    do i=1,imax
  35.       call findem(ang(i),trigs(i,:),40)
  36.    end do
  37.    call cpu_time(time(2))
  38.    times = times + time(2) - time(1)
  39.  
  40.    ! Output to trig.out
  41.    do i=1,imax
  42.      do j=1,3
  43.          if(trigs(i,j) > 1d100) then
  44.             write(tout,'(a1)') 'n'
  45.          elseif(abs(trigs(i,j)) > 1.0) then
  46.             write(tout,'(f10.6)') trigs(i,j)
  47.          elseif(abs(trigs(i,j)) < 0.1) then
  48.             write(tout,'(f10.8)') trigs(i,j)
  49.          else
  50.             write(tout,'(f9.7)') trigs(i,j)
  51.          end if
  52.          write(10,'(a)',advance='no') tout
  53.       end do
  54.       write(10,*)" "
  55.    end do
  56.  
  57.    ! Print timings and cleanup
  58.    print *,"computation took",times/real(imax,real64)*1e6,"us"
  59.    close(10); close(12)
  60.    deallocate(ang);deallocate(trigs)
  61.  
  62.  contains
  63.    !> @brief compute sine/cosine/tangent
  64.    subroutine findem(a,t,iters)
  65.      real(real64), intent(in) :: a
  66.      real(real64), intent(inout) :: t(3)
  67.      integer, intent(in) :: iters
  68. ! local variables
  69.      real(real64), parameter :: deg2rad = 1.745329252e-2
  70.      real(real64), parameter :: angles(60) = &
  71.        [ 7.8539816339744830962e-01_real64, 4.6364760900080611621e-01_real64, &
  72.          2.4497866312686415417e-01_real64, 1.2435499454676143503e-01_real64, &
  73.          6.2418809995957348474e-02_real64, 3.1239833430268276254e-02_real64, &
  74.          1.5623728620476830803e-02_real64, 7.8123410601011112965e-03_real64, &
  75.          3.9062301319669718276e-03_real64, 1.9531225164788186851e-03_real64, &
  76.          9.7656218955931943040e-04_real64, 4.8828121119489827547e-04_real64, &
  77.          2.4414062014936176402e-04_real64, 1.2207031189367020424e-04_real64, &
  78.          6.1035156174208775022e-05_real64, 3.0517578115526096862e-05_real64, &
  79.          1.5258789061315762107e-05_real64, 7.6293945311019702634e-06_real64, &
  80.          3.8146972656064962829e-06_real64, 1.9073486328101870354e-06_real64, &
  81.          9.5367431640596087942e-07_real64, 4.7683715820308885993e-07_real64, &
  82.          2.3841857910155798249e-07_real64, 1.1920928955078068531e-07_real64, &
  83.          5.9604644775390554414e-08_real64, 2.9802322387695303677e-08_real64, &
  84.          1.4901161193847655147e-08_real64, 7.4505805969238279871e-09_real64, &
  85.          3.7252902984619140453e-09_real64, 1.8626451492309570291e-09_real64, &
  86.          9.3132257461547851536e-10_real64, 4.6566128730773925778e-10_real64, &
  87.          2.3283064365386962890e-10_real64, 1.1641532182693481445e-10_real64, &
  88.          5.8207660913467407226e-11_real64, 2.9103830456733703613e-11_real64, &
  89.          1.4551915228366851807e-11_real64, 7.2759576141834259033e-12_real64, &
  90.          3.6379788070917129517e-12_real64, 1.8189894035458564758e-12_real64, &
  91.          9.0949470177292823792e-13_real64, 4.5474735088646411896e-13_real64, &
  92.          2.2737367544323205948e-13_real64, 1.1368683772161602974e-13_real64, &
  93.          5.6843418860808014870e-14_real64, 2.8421709430404007435e-14_real64, &
  94.          1.4210854715202003717e-14_real64, 7.1054273576010018587e-15_real64, &
  95.          3.5527136788005009294e-15_real64, 1.7763568394002504647e-15_real64, &
  96.          8.8817841970012523234e-16_real64, 4.4408920985006261617e-16_real64, &
  97.          2.2204460492503130808e-16_real64, 1.1102230246251565404e-16_real64, &
  98.          5.5511151231257827021e-17_real64, 2.7755575615628913511e-17_real64, &
  99.          1.3877787807814456755e-17_real64, 6.9388939039072283776e-18_real64, &
  100.          3.4694469519536141888e-18_real64, 1.7347234759768070944e-18_real64]
  101.      real(real64), parameter :: kvalues(33) = &
  102.        [ 0.70710678118654752440e+00_real64, 0.63245553203367586640e+00_real64, &
  103.          0.61357199107789634961e+00_real64, 0.60883391251775242102e+00_real64, &
  104.          0.60764825625616820093e+00_real64, 0.60735177014129595905e+00_real64, &
  105.          0.60727764409352599905e+00_real64, 0.60725911229889273006e+00_real64, &
  106.          0.60725447933256232972e+00_real64, 0.60725332108987516334e+00_real64, &
  107.          0.60725303152913433540e+00_real64, 0.60725295913894481363e+00_real64, &
  108.          0.60725294104139716351e+00_real64, 0.60725293651701023413e+00_real64, &
  109.          0.60725293538591350073e+00_real64, 0.60725293510313931731e+00_real64, &
  110.          0.60725293503244577146e+00_real64, 0.60725293501477238499e+00_real64, &
  111.          0.60725293501035403837e+00_real64, 0.60725293500924945172e+00_real64, &
  112.          0.60725293500897330506e+00_real64, 0.60725293500890426839e+00_real64, &
  113.          0.60725293500888700922e+00_real64, 0.60725293500888269443e+00_real64, &
  114.          0.60725293500888161574e+00_real64, 0.60725293500888134606e+00_real64, &
  115.          0.60725293500888127864e+00_real64, 0.60725293500888126179e+00_real64, &
  116.          0.60725293500888125757e+00_real64, 0.60725293500888125652e+00_real64, &
  117.          0.60725293500888125626e+00_real64, 0.60725293500888125619e+00_real64, &
  118.          0.60725293500888125617e+00_real64 ]
  119.     real(real64) :: beta, c, c2, factor, poweroftwo, s
  120.     real(real64) :: s2, sigma, sign_factor, theta, angle
  121.     integer :: j
  122.  
  123. ! scale to radians
  124.     beta = a*deg2rad
  125. ! ensure angle is shifted to appropriate range
  126.     call angleShift(beta, -PI, theta)
  127. ! check for signs
  128.     if( theta < -half*PI) then
  129.        theta = theta + PI
  130.        sign_factor = -1.0_real64
  131.     else if( half*PI < theta) then
  132.        theta = theta - PI
  133.        sign_factor = -1.0_real64
  134.     else
  135.        sign_factor = +1.0_real64
  136.     endif
  137.  
  138. ! set up some initializations...    
  139.     c = 1.0_real64
  140.     s = 0.0_real64
  141.     poweroftwo = 1.0_real64
  142.     angle = angles(1)
  143.  
  144. ! run for some iterations
  145.     do j=1,iters
  146.        sigma = merge(-1.0_real64, +1.0_real64, theta <  0.0_real64)
  147.        factor = sigma*poweroftwo
  148.  
  149.        c2 = c - factor*s
  150.        s2 = factor*c + s
  151.        c = c2
  152.        s = s2
  153. ! update remaining angle
  154.        theta = theta - sigma*angle
  155.  
  156.        poweroftwo = poweroftwo*0.5_real64
  157.        if(j+1 > 60) then
  158.           angle = angle * 0.5_real64
  159.        else
  160.           angle = angles(j+1)
  161.        endif
  162.     enddo !- j
  163.  
  164.     if(iters > 0) then
  165.        c = c*Kvalues(min(iters,33))
  166.        s = s*Kvalues(min(iters,33))
  167.     endif
  168.     c = c*sign_factor
  169.     s = s*sign_factor
  170.  
  171.     t = [s, c, s/c]
  172.    end subroutine findem
  173.  
  174.    subroutine angleShift(alpha, beta, gamma)
  175.      real(real64), intent(in) :: alpha, beta
  176.      real(real64), intent(out) :: gamma
  177.      if(alpha < beta) then
  178.         gamma = beta - mod(beta - alpha, TAU) + TAU
  179.      else
  180.         gamma = beta + mod(alpha - beta, TAU)
  181.      endif
  182.    end subroutine angleShift
  183.  
  184. end program SinCosTan
Advertisement
RAW Paste Data Copied
Advertisement