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 & time it
times = 0.0_real64
call system_clock(cnt(1)) ! <-- 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) > 1d100) then
write(tout,'(a1)') 'n'
elseif(abs(trigs(j,i)) > 1.0) then
write(tout,'(f10.6)') trigs(j,i)
elseif(abs(trigs(j,i)) < 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
!> @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) = &
[ 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 30 iterations (should be good enough, need testing)
do j=1,n
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(n > 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 < beta) then
gamma = beta - mod(beta - alpha, TAU) + TAU
else
gamma = beta + mod(alpha - beta, TAU)
endif
end subroutine angleShift
end program SinCosTan
If you didn't have to significant figures precision requirements, I would use Bhaskara’s approximation for the Sine http://scholarworks.umt.edu/cgi/viewcontent.cgi?article=1313&context=tme — an incredible feat for the 7th century! The maximum error is −0.001631765 (downwards) and 0.0013436967 (upwards).
– sergiol – 2017-03-07T00:08:26.800I can't use scientific notation and have output of
1.000000 0.000000 n
at the same time. Also why are you so picky about howNaN
is output? – gggg – 2018-03-07T16:57:33.3731Are we supposed to time only the trig calculations, or include the io in the timing? – gggg – 2018-03-07T20:03:45.060
Does the output have to be on a single line? It looks so pretty when it's formated with an enter key... Also, is there a specific date at which a winner is picked? – Ephraim – 2014-06-09T08:35:41.867
@Ephraim what do you mean by formatted with an enter key? no, there isn't a specific date. I really need to test all these solutions, but I haven't made the test input yet ;( – None – 2014-06-09T13:04:36.673
@professorfish - see the output in my answer. Every
sin
,cos
andtan
is on a new line. Do I need to change it to output the answers onto a single line? – Ephraim – 2014-06-09T15:59:41.0972@Ephraim The output format really doesn't matter (this isn't code-golf) as long as it outputs the sin cos and tan for every angle and they are separate – None – 2014-06-09T16:02:50.067
Does the output have to be rounded to 7 significant figures? or can it just be correct to that level? – Οurous – 2014-06-12T09:15:02.047
@Ourous it has to be rounded – None – 2014-06-12T10:32:04.970