Szybkie obliczenia trygonometryczne

Twoim zadaniem jest stworzenie programu, który może obliczyć sinus, cosinus i styczną kąta w stopniach.

Zasady

• Brak wbudowanych funkcji trygonometrii (nawet siecznych, cosecant i cotangent, jeśli ma je Twój język).
• Możesz użyć tabel odnośników, ale ich całkowity rozmiar nie może przekraczać 3000 elementów (dla wszystkich trzech razem wziętych operacji). Proszę odczytaj tabele z pliku (np. trig.lookup), Aby nie myliły kodu.
• Brak dostępu do sieci.
• Musisz poprawnie zaokrąglić wynik, jak wyjaśniono poniżej. Nie używaj podłogi ani sufitu.
• Możesz użyć dowolnej metody do obliczenia wartości, na przykład ułamków ciągłych , o ile jest to poprawne do 7 cyfr znaczących.
• Twój kod musi być w stanie sam się zmierzyć. Wyklucz operacje we / wy pliku z twojego czasu - więc po prostu określ czas funkcji, które wykonują wyzwalanie i dowolne zaokrąglanie.
• Muszę być w stanie uruchomić Twój kod. Proszę zamieścić link do darmowego kompilatora / interpretera i podać instrukcje potrzebne do skompilowania / uruchomienia kodu (np. Jakie opcje przekazać do GCC).
• Obowiązują standardowe luki .

Format wejściowy

• Czytaj z pliku o nazwie, trig.inchyba że Twój język nie obsługuje plików I / O.
• Kąty zawierają się w przedziale od 0 do 360 włącznie.
• Dane wejściowe będą się składały z kątów do dziesięciu cyfr znaczących w liczbach dziesiętnych, oddzielonych nowymi wierszami. Na przykład:

90,00000000
74,54390000
175,5000000

Format wyjściowy

• Dla każdego podanego kąta należy wyprowadzić jego sinus, cosinus i styczną do 7 znaczących liczb, oddzielonych spacjami, w jednym wierszu. Użyj „notacji naukowej”, np. 1.745329E-5Dla tan 0.001lub 1.000000E+0dla sin 90.
• Oznacz nieskończoność lub NaN n, na przykład wyjściem dla 90.00000000powinno być 1.000000 0.000000 n.
• Jeśli wejście ma trzy kąty oddzielone znakiem nowej linii, wynik powinien składać się z trzech linii, z których każda zawiera sinus, cosinus i styczną.
• Nie możesz wyprowadzać niczego innego.
• Wyjście do pliku o nazwie, trig.outchyba że twój język nie obsługuje I / O pliku.

Punktacja

• najszybszy kod . Wyzwanie polega na napisaniu programu, który oblicza te trzy wartości tak szybko, jak to możliwe. Najszybszy czas wygrywa.
• Każdy otrzyma ten sam testowy sygnał wejściowy pod wieloma kątami.
• Czasy zostaną zapisane na moim komputerze.
• Twój wynik to średnia z trzech przebiegów na tym samym wejściu (oczywiście nie możesz zapisać niczego pomiędzy przebiegami).
• Czas kompilacji nie jest wliczony. To wyzwanie dotyczy bardziej użytej metody niż języka. (Gdyby ktoś mógł mi wskazać, jak wykluczyłbym czas kompilacji dla języków takich jak Java, byłbym bardzo wdzięczny)
• Moja maszyna to instalacja Ubuntu 14.04. Statystyki procesora są na Pastebin (uzyskane przez uruchomienie cat /proc/cpuinfo).
• Zmienię twój czas na twoją odpowiedź, kiedy ją przetestuję.

Fortran 90

I zatrudnić CORDIC metodę z pre-tabelaryczne tablicy 60 wartości arctan (patrz artykuł Wiki Szczegółowe informacje o tym, dlaczego to jest konieczne).

Ten kod wymaga pliku, trig.in w którym wszystkie wartości nowego wiersza będą przechowywane w tym samym folderze, co plik wykonywalny Fortran. Kompilowanie tego jest

gfortran -O3 -o file file.f90

gdzie filejest SinCosTan.f90podana nazwa pliku (prawdopodobnie byłoby to najłatwiejsze, chociaż nie jest konieczne dopasowanie nazwy programu i nazwy pliku). Jeśli masz kompilator Intel, zalecamy użycie

ifort -O3 -xHost -o file file.f90

jak -xHost (który nie istnieje dla gfortran) zapewnia optymalizacje wyższego poziomu dostępne dla twojego procesora.

Moje testy dały mi około 10 mikrosekund na obliczenie podczas testowania 1000 losowych kątów za pomocą gfortran 4.4 (4.7 lub 4.8 jest dostępne w repozytoriach Ubuntu) i około 9,5 mikrosekund przy użyciu ifort 12.1. Testowanie tylko 10 losowych kątów da nieokreślony czas przy użyciu procedur Fortrana, ponieważ procedura pomiaru czasu jest dokładna do milisekundy, a prosta matematyka mówi, że uruchomienie wszystkich 10 liczb powinno zająć 0,100 milisekund.

EDIT Najwyraźniej mierzyłem czas IO, który (a) spowodował, że czas był dłuższy niż to konieczne i (b) jest sprzeczny z punktorem 6. Zaktualizowałem kod, aby to odzwierciedlić. Odkryłem również, że użycie kind=8liczby całkowitej z wewnętrznym podprogramem system_clockzapewnia dokładność w mikrosekundach.

Dzięki temu zaktualizowanemu kodowi obliczam teraz każdy zestaw wartości funkcji trygonometrycznych w około 0,3 mikrosekundy (znaczące cyfry na końcu zmieniają się między biegami, ale stale unosi się w pobliżu 0,31 us), co znacznie zmniejsza w porównaniu z poprzednim iteracja, która mierzyła czas 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

Python 2.7.x lub Java (wybierz)

Bezpłatny interpreter języka Python można pobrać stąd .
Bezpłatny interpreter Java można pobrać stąd .

Program może pobierać dane wejściowe zarówno z pliku o nazwie trig.in znajdującej się w tym samym katalogu, co plik programu. Dane wejściowe są oddzielone znakami nowej linii.

Pierwotnie robiłem to w Pythonie, ponieważ - cóż, uwielbiam Pythona. Ale ponieważ chcę również wygrać, przepisałem go później w Javie ...

Wersja Python: Mam około 21 µs na uruchomienie na moim komputerze. Dostałem około 32µs, gdy uruchomiłem go na IDEone .

Wersja Java: Dostaję około 0,4 µs na uruchomienie na moim komputerze i 1,8 µs na IDEone .

Dane techniczne komputera:

• Aktualizacja systemu Windows 8.1 1 64-bit z procesorem Intel Core i7-3632QM - 2,2 GHz

Test:

• Czas na metę”jest łączny czas potrzebny do obliczenia sin, cosi tanwszystkie kąty wejściowych.

• Dane wejściowe testu użyte dla obu są następujące:

  90,00000000  
  74,54390000  
  175,5000000  
  3600000.000  
  

O Kodeksie:
Założeniem tego programu było oszacowanie sini cosużycie wielomianów Taylora z 14 terminami, co, jak obliczyłem, było konieczne, aby oszacować błąd na mniej niż 1e-8. Stwierdziłem jednak, że obliczenie było szybsze sinniż cos, więc postanowiłem zamiast tego obliczyć cosza pomocącos=sqrt(1-sin^2)

Wersja Python:

import math
import timeit
import sys
import os
from functools import partial

#Global Variabls
pi2 = 6.28318530718
piover2 = 1.57079632679

#Global Lists
factorialList = [1.0,
                 -6.0,
                 120.0,
                 -5040.0,
                 362880.0,
                 -39916800.0,
                 6227020800.0,
                 -1307674368000.0,
                 355687428096000.0,
                 -121645100408832000.0,
                 51090942171709440000.0,
                 -25852016738884976640000.0,
                 15511210043330985984000000.0,
                 -10888869450418352160768000000.0,
                 8841761993739701954543616000000.0]

#simplifies angles and converts them to radians
def torad(x):  
    rev = float(x)/360.0
    if (rev>1) or (rev<0):
        return (rev - math.floor(rev))*pi2
    return rev*pi2

def sinyield(x):
    squared = x*x
    for n in factorialList:
        yield x/n
        x*=squared

def tanfastdivide(sin, cos):
    if (cos == 0):
        return "infinity"  
    return sin/cos

#start calculating sin cos and tan
def calcyield(outList):
    for angle in outList[0]:
        angle = torad(angle)
        sin = round(math.fsum(sinyield(angle)), 7)
        cos=math.copysign(math.sqrt(1-(sin*sin)),(piover2-angle))
        yield sin
        yield cos
        yield tanfastdivide(sin, cos) #tan

def calculations(IOList):
    calcyieldgen = calcyield(IOList)
    for angle in IOList[0]:
        IOList[1].append(next(calcyieldgen))
        IOList[2].append(next(calcyieldgen))
        IOList[3].append(next(calcyieldgen))
    return IOList

#Begin input from file
def ImportFile():
    try:
        infile = open("trig.in", "r")
    except:
        infile = sys.stdin
    inList = [[], [], [], []]

    #take input from file
    for line in infile:
        inList[0].extend([float(line)])
    return inList

#Begin output to file
def OutputResults(outList):
    try:
        outfile = open("trig.out", "w")
        PrintOutput(outList, outfile)    
    except:
        print 'Failed to write to file. Printing to stdout instead...'
    finally:
        PrintOutput(outList, sys.stdout)

def PrintOutput(outList, outfile):
    #outList[0][i] holds original angles
    #outList[1][i] holds sin values
    #outList[2][i] holds cos values
    #outList[3][i] holds tan values
    outfile.write('-----------------------------------------------------\n')
    outfile.write('                    TRIG RESULTS                     \n')
    outfile.write('-----------------------------------------------------\n')
    for i in range(len(outList[0])):
        if (i):
            outfile.write('\n')
        outfile.write("For angle: ")
        outfile.write(str(outList[0][i]))
        outfile.write('\n    ')
        outfile.write("Sin: ")
        outfile.write(str('%.7E' % float(outList[1][i])))
        outfile.write('\n    ')
        outfile.write("Cos: ")
        outfile.write(str('%.7E' % float(outList[2][i])))
        outfile.write('\n    ')
        outfile.write("Tan: ")
        outfile.write(str('%.7E' % float(outList[3][i])))


#Run the Program first
inList = ImportFile()
OutputResults(calculations(inList))

#Begin Runtime estimation
def timeTest(IOList):
    for output in calcyield(IOList):
        pass
def baselined(inList):
    for x in inList:
        pass

totime = timeit.Timer(partial(timeTest, inList))
baseline = timeit.Timer(partial(baselined, inList))
print '\n-----------------------------------------------------'
print '                    TIME RESULTS:                    '
print '-----------------------------------------------------'
OverheadwithCalcTime =  min(totime.repeat(repeat=10, number=10000))
Overhead = min(baseline.repeat(repeat=1, number=10000))
estimatedCalcTime = (OverheadwithCalcTime - Overhead)
estimatedTimePerAngle = estimatedCalcTime/len(inList)
estimatedTimePerCalc = estimatedTimePerAngle/3
print ' Estimated CalcTime+Overhead:.....', '%.10f' % (OverheadwithCalcTime*100), 'µsec'
print ' Estimated Overhead Time:..........', '%.10f' % (Overhead*100), 'µsec'
print ''
print ' Estimated CalcTime/Run:..........', '%.10f' % (estimatedCalcTime*100), 'µsec'
print ' Estimated CalcTime/Angle:.........', '%.10f' % (estimatedTimePerAngle*100), 'µsec'
print ' Estimated CalcTime/Cacluation:....', '%.10f' % (estimatedTimePerCalc*100), 'µsec'
print '-----------------------------------------------------'
print "                   COOL, IT WORKED!                  "
print '-----------------------------------------------------'

Wersja Java:

import java.io.FileNotFoundException;
import java.io.File;
import java.io.PrintStream;
import java.text.DecimalFormat;
import java.text.NumberFormat;
import java.util.ArrayList;
import java.util.Scanner;
import java.lang.Math;

class Trig {
   /**
    *Global Variables
    **/
    static final double pi2 = 6.28318530718;
    public long totalTime = 0L;
    static final double piover2 = 1.57079632679;
    static final double plusinfty = Double.POSITIVE_INFINITY;
    static final double minusinfty = Double.NEGATIVE_INFINITY;
    static final double factoriallist[] =
                             new double[]{
                         -6.0,120.0,-5040.0,362880.0,-39916800.0,
                         6227020800.0,-1307674368000.0,355687428096000.0,
                        -121645100408832000.0,51090942171709440000.0,
                        -25852016738884976640000.0,
                         15511210043330985984000000.0,
                        -10888869450418352160768000000.0,
                         8841761993739701954543616000000.0
                         };
//Begin Program
    public static void main(String[] args) {
        Trig mytrig = new Trig();
        double[] input = mytrig.getInput();
        double[][] output = mytrig.calculate(input);
        mytrig.OutputResults(output);
        Trig testIt = new Trig();
        testIt.timeIt(input);
    }

//Begin Timing
    public void timeIt(double[] input) {
        double overhead = 0L;
        totalTime = 0L;

        for (int i = 0; i < 1000001; i++) {
            calculate(input);
            if (i == 0) {
                overhead = totalTime;
                totalTime = 0L;
            }
        }
        double averageTime = ((Double.valueOf(totalTime-overhead))/1000000.0);
        double perAngleTime = averageTime/input.length;
        double perOpperationTime = perAngleTime/3;
        NumberFormat formatter = new DecimalFormat("0000.0000");
        System.out.println("\n-----------------------------------------------------");
        System.out.println("                    TIME RESULTS:                    ");
        System.out.println("-----------------------------------------------------");
        System.out.println("Average Total  Runtime:.......... " + formatter.format(averageTime) + " nsec");
        System.out.println("                                = " + formatter.format(averageTime/1000) + " usec\n");
        System.out.println("Average Runtime Per Angle:....... " + formatter.format(perAngleTime) + " nsec");
        System.out.println("                                = " + formatter.format(perAngleTime/1000) + " usec\n");
        System.out.println("Average Runtime Per Opperation:.. " + formatter.format(perOpperationTime) + " nsec");
        System.out.println("                                = " + formatter.format(perOpperationTime/1000) + " usec");
    }

//Begin Input
    double[] getInput() {
        Scanner in;
        ArrayList<Double> input = new ArrayList<Double>();
        try {
            in = new Scanner(new File("trig.in"));
        } catch (FileNotFoundException e) {
            new FileNotFoundException("Failed to read from file. Reading from stdin instead...").printStackTrace();
            in= new Scanner(System.in);
        }
        while (in.hasNextLine()) {
            Double toadd = Double.parseDouble(in.nextLine());
            input.add(toadd);   
        }
        in.close();
        double[] returnable = new double[input.size()];
        for(int i = 0; i < input.size(); i++) {returnable[i] = input.get(i);}
        return returnable;
    }

//Begin OutputStream Choice
    void OutputResults(double[][] outList) {
        PrintStream out;
        try {
            out = new PrintStream("trig.out");
            PrintOutput(outList, out);
            PrintOutput(outList, System.out);
        } catch (FileNotFoundException e) {
            new FileNotFoundException("Failed to write to file. Printing to stdout instead...").printStackTrace();
            PrintOutput(outList, System.out);
        }
    }

//Begin Output
    static void PrintOutput(double[][] outList, PrintStream out) {
        /**
         *outList[0][i] holds original angles
         *outList[1][i] holds sin values
         *outList[2][i] holds cos values
         *outList[3][i] holds tan values
         */
        NumberFormat formatter = new DecimalFormat("0.0000000E0");
        out.println("-----------------------------------------------------");
        out.println("                    TRIG RESULTS                     ");
        out.println("-----------------------------------------------------");
        for (int i=0; i<outList[0].length; i++) {
            out.println("For Angle: " + outList[0][i]);

            out.println("      sin: " + formatter.format(outList[1][i]));
            out.println("      cos: " + formatter.format(outList[2][i]));
            if (Double.valueOf(outList[3][i]).isInfinite() || Double.valueOf(outList[3][i]).isNaN()) {
                out.println("      tan: " + outList[3][i]);
            }
            else out.println("      tan: " + formatter.format(outList[3][i]));
        }
        if (out != System.out) out.close();
    }

//Begin Calculations
    double[][] calculate(double[] IOList) {
        double[][] outList = new double[4][IOList.length];
        double sin;
        double cos;
        double tan;
        double rads;
        int i = 0;
        long calctime = 0L;
        long startTime;
        long endTime;
        for (double input : IOList) {
            startTime = System.nanoTime();
            rads = toRad(input);
            sin=sin(rads);
            cos = ((piover2-rads)>=0) ? Math.sqrt((1.0-(sin*sin))) : -Math.sqrt((1.0-(sin*sin)));
            tan = (cos!=0.0d) ? sin/cos : ((sin>0) ? plusinfty : minusinfty);
            endTime = System.nanoTime();
            calctime = calctime + endTime - startTime;
            outList[0][i] = input;
            outList[1][i] = sin;
            outList[2][i] = cos;
            outList[3][i] = tan;
            i++;
        }
        totalTime = totalTime + calctime;
        return outList;
    }

//Convert Degrees to Radians
    double toRad(double deg) {
        double rev = deg/360.0;
        return (rev>1 || rev<0) ? Math.abs(rev - ((int)rev))*pi2 : rev*pi2;
    }

//Get sin
    double sin(double angle) {
        double sqr = angle*angle;
        double value = angle;
        for (double fct : factoriallist) {
            value += ((angle*=sqr)/fct);
        }
        return ((long)((value + Math.copySign(0.0000000005d, value))*10000000.0))/10000000.0;
    }   
}

Oktawa (lub Matlab) i C.

Trochę skomplikowany proces kompilacji, ale rodzaj nowatorskiego podejścia, a wyniki były zachęcające.

Podejście polega na generowaniu przybliżonych kwadratowych wielomianów dla każdego stopnia. Tak więc stopień = [0, 1), stopień = [1, 2), ..., stopień = [359, 360) będzie miał inny wielomian.

Oktawa - część budowlana

Oktawa jest publicznie dostępna - Google download octave.

Określa to najlepiej dopasowany wielomian kwadratowy do każdego stopnia.

Zapisz jako build-fast-trig.m:

format long;
for d = 0:359
    x = (d-1):0.1:(d+1);
    y = sin(x / 360 * 2 * pi);
    polyfit(x, y, 2)
endfor

C - część budynku

Konwertuje to podwójnie w formacie tekstowym na natywny format binarny w twoim systemie.

Zapisz jako build-fast-trig.c:

#include <stdio.h>

int main()
{
    double d[3];

    while (scanf("%lf %lf %lf", d, d + 1, d + 2) == 3)
        fwrite(d, sizeof(double), 3, stdout);

    return 0;
}

Skompilować:

gcc -o build-fast-trig build-fast-trig.c

Generowanie pliku współczynników

Biegać:

octave build-fast-trig.m | grep '^ ' | ./build-fast-trig > qcoeffs.dat

Teraz mamy qcoeffs.dat jako plik danych do użycia dla rzeczywistego programu.

C - część szybkiego wyzwalania

Zapisz jako fast-trig.c:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

#define INPUT    "qcoeffs.dat"

#define DEGREES    360

typedef struct {double a, b, c;} QCOEFFS;

double normalize(double d)
{
    if (d < 0.0)
        d += ceil(d / -(double)DEGREES) * (double)DEGREES;

    if (d >= (double)DEGREES)
        d -= floor(d / (double)DEGREES) * (double)DEGREES;

    return d;
}

int main()
{
    FILE *f;
    time_t tm;
    double d;
    QCOEFFS qc[DEGREES];

    if (!(f = fopen(INPUT, "rb")) || fread(qc, sizeof(QCOEFFS), DEGREES, f) < DEGREES)
    {
        fprintf(stderr, "Problem with %s - aborting.", INPUT);
        return EXIT_FAILURE;
    }
    fclose(f);

    tm = -clock();

    while (scanf("%lf", &d) > 0)
    {
        int i;
        double e, f;

        /* sin */
        d = normalize(d);
        i = (int)d;
        e = (qc[i].a * d + qc[i].b) * d + qc[i].c;

        /* cos */
        d = normalize((double)DEGREES / 4.0 - d);
        i = (int)d;
        f = (qc[i].a * d + qc[i].b) * d + qc[i].c;

        /* tan = sin / cos */

        /* output - format closest to specs, AFAICT */
        if (d != 0.0 && d != 180.0)
            printf("%.6e %.6e %.6e\n", e, f, e / f);
        else
            printf("%.6e %.6e n\n", e, f);
    }

    tm += clock();

    fprintf(stderr, "time: %.3fs\n", (double)tm/(double)CLOCKS_PER_SEC);    

    return EXIT_SUCCESS;
}

Skompilować:

gcc -o fast-trig fast-trig.c -lm

Biegać:

./fast-trig < trig.in > trig.out

Będzie odczytywał trig.in, zapisywał się trig.outi drukował, aby pocieszać upływający czas z milisekundową precyzją.

W zależności od zastosowanych metod testowania może się nie powieść na niektórych danych wejściowych, np .:

$ ./fast-trig 
0
-6.194924e-19 1.000000e+00 -6.194924e-19

Prawidłowe wyjście powinno być 0.000000e+00 1.000000e+00 0.000000e+00. Jeśli wyniki zostaną sprawdzone przy użyciu ciągów, dane wejściowe zakończą się niepowodzeniem, jeśli zostaną zweryfikowane przy użyciu błędu bezwzględnego, np fabs(actual - result) < 1e-06. Dane wejściowe przejdą.

Maksymalny błąd bezwzględny dla sini coswynosił ≤ 3e-07. Ponieważ tan, ponieważ wynik nie jest ograniczony do ± 1 i można podzielić stosunkowo dużą liczbę przez względnie małą liczbę, błąd bezwzględny może być większy. Od -1 ≤ tan (x) ≤ +1 maksymalny błąd bezwzględny wynosił ≤ 4e-07. Dla tan (x)> 1 i tan (x) <-1 maksymalny względny błąd , np. fabs((actual - result) / actual)Zwykle wynosił <1e-06, dopóki nie osiągniesz obszaru (90 ± 5) lub (270 ± 5) stopni, a następnie błąd się pogarsza.

W badaniach, średni czas był za jednego wejścia (1,053 ± 0,007 ms), który na moim komputerze było około 0,070 mikrosekundy szybciej niż rodzimy sini cos, tanjest zdefiniowana w ten sam sposób.

Kobra

class Trig
    const mod as float = 0.0174532925199433f #0.0174532925199432957692369076848861271344287188854172f
    var time as System.Diagnostics.Stopwatch = System.Diagnostics.Stopwatch()
    var file as String[] = File.readAllLines('trig.in')
    var sin_out as float[] = float[](1)
    var cos_out as float[] = float[](1)
    var tan_out as float[] = float[](1)
    def main
        .compute(@[1f])
        .time.reset
        input = float[](.file.length)
        for num, line in .file.numbered, input[num] = float.parse(line)
        .compute(input)
        for num in .file.length, .file[num] = (.sin_out[num].toString('0.000000E+0') + ' ' + .cos_out[num].toString('0.000000E+0') + ' ' + .tan_out[num].toString('0.000000E+0'))
        File.writeAllLines('trig.out', .file)
        print .time.elapsed
    def compute(list as float[])
        .sin_out = float[](list.length)
        .cos_out = float[](list.length)
        .tan_out = float[](list.length)
        .time.start
        for index in list.length
            degrees as float = list[index]
            #add `degrees %= 360` for numbers greater than 360
            rad as float = sin as float = degrees * .mod
            two as float = rad * rad
            sin -= (rad *= two) / 6
            sin += (rad *= two) / 120
            sin -= (rad *= two) / 5040
            sin += (rad *= two) / 362880
            sin -= (rad *= two) / 39916800
            sin += (rad *= two) / 6227020800
            sin -= (rad *= two) / 1307674368000
            sin += (rad *= two) / 355687428096000
            sin -= (rad *= two) / 121645100408832000
            sin += (rad *= two) / 51090942171709440000f
            sin -= (rad *= two) / 25852016738884976640000f
            sin += (rad *= two) / 15511210043330985984000000f
            sin -= (rad *= two) / 10888869450418352160768000000f
            sin += (rad *= two) / 8841761993739701954543616000000f
            cos as float = (1 - (sin * sin)).sqrt * ((degrees - 180).abs - 90).sign
            if cos.isNaN, cos = 0
            .tan_out[index] = Math.round((sin / cos) * 10000000) / 10000000
            .sin_out[index] = Math.round(sin * 10000000) / 10000000
            .cos_out[index] = Math.round(cos * 10000000) / 10000000
        .time.stop

Skompiluj to z cobra filename -turbo

Testy: AMD FX6300 @ 5.1GHz

• Test 360 * 10000 zastosowany w odpowiedzi C działa w ciągu 365 ms (w porównaniu do 190 ms)

• 4-wejściowy test używany w odpowiedziach na Python i Java przebiega w 0,32µs (vs 30µs, 3µs)

• Test 1000 losowych kątów zastosowany w odpowiedzi Fortrana przebiega przy 100ns na kąt (vs 10µs)

do

Oto moja próba. Działa to tak:

Zbuduj tabelę wszystkich wartości sin (x) od 0 do 450 stopni. Równolegle są to wszystkie wartości cos (x) od -90 do 360 stopni. Przy 2926 elementach jest wystarczająco dużo miejsca na wartość co 1 / 6,5 stopnia. Jednostka programu wynosi zatem 1 / 6,5 stopnia, a na ćwierć obrotu jest 585 jednostek.

Przelicz stopnie wejściowe na jednostki programu (pomnóż przez 6.5==110.1 binary. ) Znajdź najbliższe wartości sin i cos z tabeli. następnie zamień pozostałą część danych wejściowych (dx) na radiany.

zastosuj wzór sin(x+dx) == sin x +(d(sin x)/dx)*dx. zwróć uwagę, że(d(sin x)/dx)==cos x, tylko, jeśli użyjemy radianów.

niestety nie jest to wystarczająco dokładne, więc wymagany jest inny termin, oparty na następnej pochodnej. d2(sin x)/dx2 == -sin x.Należy to pomnożyć przezdx*dx/2 (nie jestem pewien, skąd bierze się współczynnik 2, ale działa).

Postępuj zgodnie z analogiczną procedurą dla cos x, a następnie oblicztan x == sin x / cos x .

Kod

Jest tu około 17 operacji zmiennoprzecinkowych. Można to nieco poprawić. Program zawiera dane wyjściowe do budowania tabel i testowania przy użyciu natywnych funkcji wyzwalania, ale algorytm nie. Dodam taktowanie i edycję, aby spełnić wymagania we / wy później (mam nadzieję, że w ten weekend.) Pasuje do wyjściowej funkcji natywnej, z wyjątkiem bardzo małych wartości sin x i cos x, które powinny być lepsze niż wyjściowa funkcja natywna za pomocą trochę poprawiania.

<#include <math.h>                                                 //only for table building and testing
int a;
double t[2926],                                                    //a table for sin x from 0 to 360+90=450deg every 1/6.5 degrees
x,                                                                 //input
s,c,                                                               //first guess sin and cos (from table)
sn,cs,                                                             //output sin and cos
pi1170=3.1415926535897932384626433832795/1170,                     // there are 1170 units of 1/6.5 degrees in a half circle 
pi180=3.1415926535897932384626433832795/180;                       // pi/180 only used for testing

main(){
  for (a=0;a<=2925;a++)t[a]=sin(a*pi1170);                         //build table. 

  scanf("%lf",&x);                                                 //input 
  printf("%le %le %le\n",sin(x*pi180),cos(x*pi180),tan(x*pi180));  //native output for testing purposes

  x*=6.5;                                                          //convert from deg to program units. 6.5=110.1 in binary, a fairly round number. 
  a=x+0.5;                                                         //a=round(x) add 0.5 to round, not truncate. Assigning to int, this is probably faster than the round function.
  x=(x-a)*pi1170;                                                  //(x-a)=dx in program units. Divide to get radians. 

  s=t[a];                                                          //get sin(a) from table
  c=t[a+585];                                                      //cos(a)=sin(a+90degrees)=sin(a+585units)
  sn=s+c*x-s*x*x/2;                                                //sin(x+dx)=sin(x)+cos(dx)-sin(dx^2/2)
  cs=c-s*x-c*x*x/2;                                                //cos(x+dx)=cos(x)-sin(dx)-cos(dx^2/2)
  printf("%le %le %le\n",sn,cs,sn/cs);                             //print sin,cos and tan=sin/cos
}