Fast Trig Calculation Fast Trigonometry Calculations

Your task is to create a program which can calculate the sine, cosine and tangent of an angle in degrees.

Rules


No built-in trigonometry functions (not even secant, cosecant and cotangent if your language has them).
You may use lookup tables, but their total size must not exceed 3000 members (for all of the three operations put together). Please make it read the tables from a file (e.g. trig.lookup) so they don't confuse the code.
No network access.
You must correctly round your output as explained below. Don't use floor or ceiling.
You may use any method to calculate the values, for example continued fractions, as long as it is correct to 7 significant figures.
Your code must be able to time itself. Exclude the file I/O operations from your time - so just time the function(s) that do(es) the trig and any rounding.
I must be able to run your code. Please post a link to a freely available compiler/interpreter and give the instructions needed to compile/run the code (e.g. what options to pass to GCC).
Standard loopholes apply.


Input Format


Read from a file called trig.in unless your language doesn't support file I/O.
Angles are between 0 and 360 inclusive.
Input will consist of angles to ten significant figures in decimal digits, separated by new lines. For example:



  90.00000000
  74.54390000
  175.5000000


Output Format


For each angle supplied, you must output its sine, cosine and tangent to 7 significant figures, separated by spaces, on a single line. Use "scientific notation", e.g. 1.745329E-5 for tan 0.001 or 1.000000E+0 for sin 90.
Denote infinity or NaN by n, for example the output for 90.00000000 should be 1.000000 0.000000 n.
If the input is three angles separated by newlines, your output should consist of three lines, each containing the sine, cosine and tangent.
You may not output anything else.
Output to a file called trig.out unless your language does not support file I/O.


Scoring


fastest-code. The challenge is to write a program which computes these three values as quickly as possible. Fastest time wins.
Everyone will receive the same test input of many angles.
Times will be recorded on my machine.
Your score is the average of three runs on the same input (you can't save anything in between runs obviously).
Compiling time not included. This challenge is more about the method used than the language. (If someone could point me to how I would exclude compilation time for languages such as Java, I'd be very grateful)
My machine is an Ubuntu 14.04 install. The processor's stats are on Pastebin (obtained by running cat /proc/cpuinfo).
I will edit your time into your answer when I've tested it.
Fortran 90

I employ the CORDIC method with a pre-tabulated array of 60 arctan values (see Wiki article for details on why that's necessary).

This code requires a file, trig.in, with all the values on newlines to be stored in the same folder as the Fortran executable. Compiling this is,

gfortran -O3 -o file file.f90


where file is whatever filename you give it (probably SinCosTan.f90 would be easiest, though it's not necessary to match program name and file name). If you have the Intel compiler, I'd recommend using

ifort -O3 -xHost -o file file.f90


as the -xHost (which doesn't exist for gfortran) provides higher level optimizations available to your processor.

My test runs were giving me about 10 microseconds per calculation when testing 1000 random angles using gfortran 4.4 (4.7 or 4.8 is available in the Ubuntu repos) and around 9.5 microseconds using ifort 12.1. Testing only 10 random angles will result in an indeterminable time using Fortran routines, as the timing routine is accurate to the millisecond and simple math says it should take 0.100 milliseconds to run all 10 numbers.



EDIT Apparently I was timing IO which (a) made the timing longer than necessary and (b) is contrary to bullet #6. I've updated the code to reflect this. I've also discovered that using a kind=8 integer with the intrinsic subroutine system_clock gives microsecond accuracy.

With this updated code, I am now computing each set of values of the trigonometric functions in about 0.3 microseconds (the significant digits in the end vary run-to-run, but it's consistently hovering near 0.31 us), a significant reduction from the previous iteration that timed IO.



program SinCosTan
   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), allocatable :: trigs(:,:), angles(:)
   real(real64) :: time(2), times, b
   character(len=12) :: tout
   integer :: i,j,ierr,amax
   integer(kind=8) :: cnt(2)

   open(unit=10,file='trig.out',status='replace')
   open(unit=12,file='CodeGolf/trig.in',status='old')
! check to see how many angles there are
   i=0
   do
      read(12,*,iostat=ierr) b
      if(ierr/=0) exit
      i=i+1
   enddo !- 
   print '(a,i0,a)',"There are ",i," angles"
   amax = i

! allocate array
   allocate(trigs(3,amax),angles(amax))

! rewind the file then read the angles into the array
   rewind(12)
   do i=1,amax
      read(12,*) angles(i)
   enddo !- i

! compute trig functions &amp; time it
   times = 0.0_real64
   call system_clock(cnt(1)) ! &lt;-- system_clock with an 8-bit INT can time to us
   do i=1,amax
      call CORDIC(angles(i), trigs(:,i), 40)
   enddo !- i
   call system_clock(cnt(2))
   times = times + (cnt(2) - cnt(1))

! write the angles to the file
   do i=1,amax
      do j=1,3
         if(trigs(j,i) &gt; 1d100) then
            write(tout,'(a1)') 'n'
         elseif(abs(trigs(j,i)) &gt; 1.0) then
            write(tout,'(f10.6)') trigs(j,i)
         elseif(abs(trigs(j,i)) &lt; 0.1) then
            write(tout,'(f10.8)') trigs(j,i)
         else
            write(tout,'(f9.7)') trigs(j,i)
         endif
         write(10,'(a)',advance='no') tout
      enddo !- j
      write(10,*)" "
   enddo !- i

   print *,"computation took",times/real(i,real64),"us per angle"
   close(10); close(12)
 contains
   !&gt; @brief compute sine/cosine/tangent
   subroutine CORDIC(a,t,n)
     real(real64), intent(in) :: a
     real(real64), intent(inout) :: t(3)
     integer, intent(in) :: n
! local variables
     real(real64), parameter :: deg2rad = 1.745329252e-2
     real(real64), parameter :: angles(60) = &amp;
       [ 7.8539816339744830962e-01_real64, 4.6364760900080611621e-01_real64, &amp;
         2.4497866312686415417e-01_real64, 1.2435499454676143503e-01_real64, &amp;
         6.2418809995957348474e-02_real64, 3.1239833430268276254e-02_real64, &amp;
         1.5623728620476830803e-02_real64, 7.8123410601011112965e-03_real64, &amp;
         3.9062301319669718276e-03_real64, 1.9531225164788186851e-03_real64, &amp;
         9.7656218955931943040e-04_real64, 4.8828121119489827547e-04_real64, &amp;
         2.4414062014936176402e-04_real64, 1.2207031189367020424e-04_real64, &amp;
         6.1035156174208775022e-05_real64, 3.0517578115526096862e-05_real64, &amp;
         1.5258789061315762107e-05_real64, 7.6293945311019702634e-06_real64, &amp;
         3.8146972656064962829e-06_real64, 1.9073486328101870354e-06_real64, &amp;
         9.5367431640596087942e-07_real64, 4.7683715820308885993e-07_real64, &amp;
         2.3841857910155798249e-07_real64, 1.1920928955078068531e-07_real64, &amp;
         5.9604644775390554414e-08_real64, 2.9802322387695303677e-08_real64, &amp;
         1.4901161193847655147e-08_real64, 7.4505805969238279871e-09_real64, &amp;
         3.7252902984619140453e-09_real64, 1.8626451492309570291e-09_real64, &amp;
         9.3132257461547851536e-10_real64, 4.6566128730773925778e-10_real64, &amp;
         2.3283064365386962890e-10_real64, 1.1641532182693481445e-10_real64, &amp;
         5.8207660913467407226e-11_real64, 2.9103830456733703613e-11_real64, &amp;
         1.4551915228366851807e-11_real64, 7.2759576141834259033e-12_real64, &amp;
         3.6379788070917129517e-12_real64, 1.8189894035458564758e-12_real64, &amp;
         9.0949470177292823792e-13_real64, 4.5474735088646411896e-13_real64, &amp;
         2.2737367544323205948e-13_real64, 1.1368683772161602974e-13_real64, &amp;
         5.6843418860808014870e-14_real64, 2.8421709430404007435e-14_real64, &amp;
         1.4210854715202003717e-14_real64, 7.1054273576010018587e-15_real64, &amp;
         3.5527136788005009294e-15_real64, 1.7763568394002504647e-15_real64, &amp;
         8.8817841970012523234e-16_real64, 4.4408920985006261617e-16_real64, &amp;
         2.2204460492503130808e-16_real64, 1.1102230246251565404e-16_real64, &amp;
         5.5511151231257827021e-17_real64, 2.7755575615628913511e-17_real64, &amp;
         1.3877787807814456755e-17_real64, 6.9388939039072283776e-18_real64, &amp;
         3.4694469519536141888e-18_real64, 1.7347234759768070944e-18_real64]
     real(real64), parameter :: kvalues(33) = &amp;
       [ 0.70710678118654752440e+00_real64, 0.63245553203367586640e+00_real64, &amp;
         0.61357199107789634961e+00_real64, 0.60883391251775242102e+00_real64, &amp;
         0.60764825625616820093e+00_real64, 0.60735177014129595905e+00_real64, &amp;
         0.60727764409352599905e+00_real64, 0.60725911229889273006e+00_real64, &amp;
         0.60725447933256232972e+00_real64, 0.60725332108987516334e+00_real64, &amp;
         0.60725303152913433540e+00_real64, 0.60725295913894481363e+00_real64, &amp;
         0.60725294104139716351e+00_real64, 0.60725293651701023413e+00_real64, &amp;
         0.60725293538591350073e+00_real64, 0.60725293510313931731e+00_real64, &amp;
         0.60725293503244577146e+00_real64, 0.60725293501477238499e+00_real64, &amp;
         0.60725293501035403837e+00_real64, 0.60725293500924945172e+00_real64, &amp;
         0.60725293500897330506e+00_real64, 0.60725293500890426839e+00_real64, &amp;
         0.60725293500888700922e+00_real64, 0.60725293500888269443e+00_real64, &amp;
         0.60725293500888161574e+00_real64, 0.60725293500888134606e+00_real64, &amp;
         0.60725293500888127864e+00_real64, 0.60725293500888126179e+00_real64, &amp;
         0.60725293500888125757e+00_real64, 0.60725293500888125652e+00_real64, &amp;
         0.60725293500888125626e+00_real64, 0.60725293500888125619e+00_real64, &amp;
         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 &lt; -half*PI) then
       theta = theta + PI
       sign_factor = -1.0_real64
    else if( half*PI &lt; 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 30 iterations (should be good enough, need testing)
    do j=1,n
       sigma = merge(-1.0_real64, +1.0_real64, theta &lt;  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 &gt; 60) then
          angle = angle * 0.5_real64
       else
          angle = angles(j+1)
       endif
    enddo !- j

    if(n &gt; 0) then
       c = c*Kvalues(min(n,33))
       s = s*Kvalues(min(n,33))
    endif
    c = c*sign_factor
    s = s*sign_factor

    t = [s, c, s/c]
   end subroutine CORDIC

   subroutine angleShift(alpha, beta, gamma)
     real(real64), intent(in) :: alpha, beta
     real(real64), intent(out) :: gamma
     if(alpha &lt; beta) then
        gamma = beta - mod(beta - alpha, TAU) + TAU
     else
        gamma = beta + mod(alpha - beta, TAU) 
     endif
   end subroutine angleShift

end program SinCosTan