Идея: поскольку | sin'x | < 1, можно считать его значения по таблице с равномерным шагом, равным требуемой точности результата.

Реализация:

diff -ruN src-org/sunset.cpp src-more/sunset.cpp
--- src-org/sunset.cpp  2007-09-16 12:04:44.000000000 +0400
+++ src-more/sunset.cpp 2007-11-07 09:55:27.000000000 +0300
@@ -45,6 +45,8 @@
 #include <omp.h>
 #endif
 #include "sunset.h"
+#include <ippvm.h>
+#include <ippm.h>

 #define MIN(x,y)    (((x) < (y)) ? (x) : (y))
 #define MAX(x,y)    (((x) < (y)) ? (y) : (x))
@@ -115,6 +117,33 @@
     }
 }

+#define SIN_TAB_SZ 16384
+#define PI2 (2*3.141592653589793f)
+float g_sinTab[2*SIN_TAB_SZ+1];
+
+void fillSinTab()
+{
+  for(int i=0;i<2*SIN_TAB_SZ+1;++i)
+  {
+    g_sinTab[i]=sinf((i-SIN_TAB_SZ)*PI2/SIN_TAB_SZ);
+  }
+}
+
+inline float tab_sinf(float v)
+{
+  int idx=(((int)(v*(SIN_TAB_SZ/PI2))) % SIN_TAB_SZ)+SIN_TAB_SZ;
+  return g_sinTab[idx];
+}
+
+inline float lin_sinf(float v)
+{
+  float i=fmodf(v,PI2)*SIN_TAB_SZ/PI2+SIN_TAB_SZ;
+  size_t idx=i;
+  float d=i-idx;
+
+  return g_sinTab[idx]*(1-d)+g_sinTab[idx+1]*d;
+}
+
 /*-----------------------------------------------------------*/

 #define DECLARE_ALIGNED_PTR(type, pName) \
@@ -320,6 +349,7 @@
 /*---------------------------------------------------------------------------*/
     if(iCurFrame == 1)
     {
+        fillSinTab();
         iWaveMeshSize = iWaveHarmNum * iAngleHarmNum;

         if(iAllocated == 1)
@@ -730,30 +760,34 @@
 !!!!!!!!!!!!!!!!! Water surface modelling !!!!!!!!!!!!!!!!!!!!!!!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 */
+                pFlTmp = flArgSin[currentthread].aptr;
+                flDerivX = 0.;
+                flDerivY = 0.;
+#if 1
                 for(t = 0; t < NKMAX; t++)
                 {
-                    OT  = flOmega[t] * flTime;
+                    OT  = flOmegaTime[t];
                     KX1 = flK[t] * flDecartX[i][j];
                     KY1 = flK[t] * flDecartY[i][j];

                     for(l = 0; l < iAngleHarmNum; l++)
                     {
                         iSinIndex1 = t * iAngleHarmNum + l;
-                        flArgSin[currentthread].aptr[iSinIndex1] = OT -
+                        pFlTmp[iSinIndex1] = OT -
                             KX1 * flAzimuthCosFi[l] - KY1 * flAzimuthSinFi[l] +
-                            flRandomPhase[t*iAngleHarmNum + l];
+                            flRandomPhase[iSinIndex1];
                     } /* end for l */
                 } /* end for t */

-                pFlTmp = flArgSin[currentthread].aptr;
-
+//                ippsSin_32f_A11(pFlTmp,pFlTmp,iWaveMeshSize);
                 #pragma ivdep
+                #pragma vector
                 for(t=0; t<iWaveMeshSize; t++)
-                    pFlTmp[t] = (float)sinf(pFlTmp[t]);
+                {
+                  pFlTmp[t] = (float)tab_sinf(pFlTmp[t]);
+                }

                 /* initialize the values of derivation */
-                flDerivX = 0.0f;
-                flDerivY = 0.0f;

                 /* dot product to compute derivation */
                 for(t = 0; t < iWaveMeshSize; t++)
@@ -761,6 +795,7 @@
                     flDerivX += pFlTmp[t] * flAmplitudeX[t];
                     flDerivY += pFlTmp[t] * flAmplitudeY[t];
                 }
+#endif

                 /* Near horizont area correction */
                 flDerivX *= P2;

Результат: gcc – медленно, icc – догоняем ippsSinf. Основной тормоз – смешанная арифметика, перевод плавающего числа в целое.