Había escrito un código hace un tiempo que intentaba calcular sin usar las funciones de la biblioteca. Ayer, estaba revisando el código anterior, e intenté hacerlo lo más rápido posible (y correcto). Aquí está mi intento hasta ahora:
const double ee = exp(1);
double series_ln_taylor(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now;
int i, flag = 1;
if ( n <= 0 ) return 1e-300;
if ( n * ee < 1 )
n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
for ( term = 1; term < n ; term *= ee, lgVal++ );
n /= term;
/* log(1 - x) = -x - x**2/2 - x**3/3... */
n = 1 - n;
now = term = n;
for ( i = 1 ; ; ){
lgVal -= now;
term *= n;
now = term / ++i;
if ( now < 1e-17 ) break;
}
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
Aquí estoy tratando de encontrar para que esté justo por encima de n, y luego agrego el valor del logaritmo de , que es menor que 1. En este punto, La expansión Taylor de se puede utilizar sin preocuparse.e a n log(1-x)
Recientemente he crecido el interés en el análisis numérico, y es por eso que no puedo evitar hacer la pregunta, ¿cuánto más rápido se puede ejecutar este segmento de código en la práctica, sin dejar de ser lo suficientemente correcto? ¿Necesito cambiar a otros métodos, por ejemplo, usando una fracción continua, como esta ?
La función proporcionada con la biblioteca estándar de C es casi 5.1 veces más rápida que esta implementación.
ACTUALIZACIÓN 1 : Usando la serie de arctanos hiperbólicos mencionados en Wikipedia , el cálculo parece ser casi 2.2 veces más lento que la función de registro de biblioteca estándar de C. Sin embargo, no he comprobado exhaustivamente el rendimiento, y para grandes cantidades, mi implementación actual parece REALMENTE lenta. Si puedo administrar, quiero verificar tanto la implementación de errores como el tiempo promedio para una amplia gama de números. Aquí está mi segundo esfuerzo.
double series_ln_arctanh(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now, sm;
int i, flag = 1;
if ( n <= 0 ) return 1e-300;
if ( n * ee < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
for ( term = 1; term < n ; term *= ee, lgVal++ );
n /= term;
/* log(x) = 2 arctanh((x-1)/(x+1)) */
n = (1 - n)/(n + 1);
now = term = n;
n *= n;
sm = 0;
for ( i = 3 ; ; i += 2 ){
sm += now;
term *= n;
now = term / i;
if ( now < 1e-17 ) break;
}
lgVal -= 2*sm;
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
Cualquier sugerencia o crítica es apreciada.
ACTUALIZACIÓN 2: según las sugerencias que se hacen a continuación, he agregado algunos cambios incrementales aquí, que es aproximadamente 2.5 veces más lento que la implementación estándar de la biblioteca. Sin embargo, esta vez lo he probado solo para enteros , para números más grandes el tiempo de ejecución aumentaría. Por ahora. Todavía no conozco técnicas para generar números dobles aleatorios , por lo que aún no está completamente comparado. Para hacer que el código sea más robusto, he agregado correcciones para los casos de esquina. El error promedio para las pruebas que hice es de alrededor de . ≤ 1 e 308 4 e - 15
double series_ln_better(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now, sm;
int i, flag = 1;
if ( n == 0 ) return -1./0.; /* -inf */
if ( n < 0 ) return 0./0.; /* NaN*/
if ( n < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
/* the cutoff iteration is 650, as over e**650, term multiplication would
overflow. For larger numbers, the loop dominates the arctanh approximation
loop (with having 13-15 iterations on average for tested numbers so far */
for ( term = 1; term < n && lgVal < 650 ; term *= ee, lgVal++ );
if ( lgVal == 650 ){
n /= term;
for ( term = 1 ; term < n ; term *= ee, lgVal++ );
}
n /= term;
/* log(x) = 2 arctanh((x-1)/(x+1)) */
n = (1 - n)/(n + 1);
now = term = n;
n *= n;
sm = 0;
/* limiting the iteration for worst case scenario, maximum 24 iteration */
for ( i = 3 ; i < 50 ; i += 2 ){
sm += now;
term *= n;
now = term / i;
if ( now < 1e-17 ) break;
}
lgVal -= 2*sm;
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
fuente
frexpLa respuesta de Kirill ya tocó una gran cantidad de temas relevantes. Me gustaría ampliar algunos de ellos basados en la experiencia práctica de diseño de la biblioteca de matemáticas. Una nota por adelantado: los diseñadores de la biblioteca matemática tienden a usar todas las optimizaciones algorítmicas publicadas, así como muchas optimizaciones específicas de la máquina, no todas las cuales se publicarán. El código frecuentemente se escribirá en lenguaje ensamblador, en lugar de usar código compilado. Por lo tanto, es poco probable que una implementación sencilla y compilada logre más del 75% del rendimiento de una implementación de biblioteca matemática de alta calidad existente, suponiendo conjuntos de características idénticas (precisión, manejo de casos especiales, informes de errores, soporte de modo de redondeo).
La precisión generalmente se evalúa en comparación con una referencia (de terceros) de mayor precisión. Las funciones de precisión simple de argumento único se pueden probar fácilmente exhaustivamente, otras funciones requieren pruebas con vectores de prueba aleatorios (dirigidos). Claramente, uno no puede calcular resultados de referencia infinitamente precisos, pero la investigación sobre el Dilema del Creador de tablas sugiere que para muchas funciones simples es suficiente calcular una referencia con una precisión de aproximadamente tres veces la precisión objetivo. Ver por ejemplo:
Vincent Lefèvre, Jean-Michel Muller, "Los peores casos para el redondeo correcto de las funciones elementales en doble precisión". En Actas, 15º Simposio IEEE sobre aritmética informática , 2001,111-118). (preimpresión en línea)
En términos de rendimiento, uno debe distinguir entre la optimización de la latencia (importante cuando se observa el tiempo de ejecución de las operaciones dependientes) y la optimización del rendimiento (relevante cuando se considera el tiempo de ejecución de las operaciones independientes). En los últimos veinte años, la proliferación de técnicas de paralelización de hardware como el paralelismo a nivel de instrucción (por ejemplo, procesadores superescalares, fuera de orden), el paralelismo a nivel de datos (por ejemplo, instrucciones SIMD) y el paralelismo a nivel de hilo (por ejemplo, hiper-threading, procesadores de múltiples núcleos) ha llevado a un énfasis en el rendimiento computacional como la métrica más relevante.
La operación fusionada de adición múltiple ( FMA ), introducida por primera vez por IBM hace 25 años, y ahora disponible en todas las arquitecturas de procesadores principales, es un componente fundamental de las implementaciones modernas de la biblioteca matemática. Proporciona reducción de errores de redondeo, brinda protección limitada contra la cancelación sustractiva y simplifica enormemente la aritmética doble-doble .
C99log()C99fma()fuente