La literatura sobre computación de la función elementalsin con tablas se refiere a la fórmula:

sin(x) = sin(Cn) * cos(h) + cos(Cn) * sin(h)

dóndex = Cn + h, Cn es una constante para la cualsin(Cn) ycos(Cn) se han calculado previamente y están disponibles en una tabla y, si se sigue el método de Gal,Cn ha sido elegido para que ambossin(Cn) ycos(Cn) están estrechamente aproximados por números de punto flotante. La cantidadh esta cerca de0.0. Un ejemplo de referencia a esta fórmula es esteartículo (página 7).

No entiendo por qué esto tiene sentido:cos(h), sin embargo, se calcula, probablemente será incorrecto en al menos 0.5 ULP para algunos valores dehy ya que está cerca de1.0, esto parece tener un efecto drástico en la precisión del resultadosin(x) cuando se calcula de esta manera.

No entiendo por qué la fórmula a continuación no se usa en su lugar:

sin(x) = sin(Cn) + (sin(Cn) * (cos(h) - 1.0) + cos(Cn) * sin(h))

Entonces las dos cantidades(cos(h) - 1.0) ysin(h) se puede aproximar con polinomios que son fáciles de precisar ya que producen resultados cercanos a cero. Los valores parasin(Cn) * (cos(h) - 1.0), cos(Cn) * sin(h) y para su suma todavía es pequeña y su precisión absoluta se expresa en ULP de la pequeña cantidad que representa la suma, por lo que sumar esta cantidad asin(Cn) está casi correctamente redondeado.

¿Me estoy perdiendo algo que hace que la fórmula anterior, popular y más simple se comporte bien también? ¿Los escritores dan por sentado que los lectores entenderán que la primera fórmula se implementa realmente como la segunda fórmula?

EDITAR: Ejemplo

Una tabla de precisión simple para calcular precisión simplesinf() ycosf() podría contener el siguiente punto en precisión simple:

         f             |        cos f          |       sin f      
-----------------------+-----------------------+---------------------
0.017967 0x1.2660bcp-6 |    0x1.ffead8p-1      |    0x1.265caep-6
                       |    (actual value:)    |    (actual value:)
                       | ~0x1.ffead8000715dp-1 | ~0x1.265cae000e6f9p-6

Las siguientes funciones son funciones especializadas de precisión simple para usar alrededor0.017967:

float sinf_trad(float x)
{
  float h = x - 0x1.2660bcp-6f;

  return 0x1.265caep-6f * cos_0(h) + 0x1.ffead8p-1f * sin_0(h);
}

float sinf_new(float x)
{
  float h = x - 0x1.2660bcp-6f;

  return 0x1.265caep-6f + (0x1.265caep-6f * cosm1_0(h) + 0x1.ffead8p-1f * sin_0(h));
}

La prueba de estas funciones entre 0.01f y 0.025f parece mostrar que la nueva fórmula da resultados más precisos:

$ gcc -std=c99 test.c && ./a.out 
relative error, traditional: 2.169624e-07, new: 1.288049e-07
sum of squares of absolute error, traditional: 6.616202e-12, new: 2.522784e-12

Tomé varios atajos, así que mirael programa completo.