<math.h> es la biblioteca estándar para realizar operaciones matemáticas avanzadas sobre tipos de punto flotante. Proporciona funciones para trigonometría, logaritmos, potencias, raíces y diversos métodos de redondeo. Debido a que la precisión de un float, un double y un long double varía significativamente, las funciones están sobrecargadas mediante sufijos: sin para double, sinf para float y sinl para long double [disponible desde C99]. Esta distinción es crucial porque el uso de la versión incorrecta puede resultar en una pérdida de precisión innecesaria o en un consumo de ciclos de CPU impropio si trabajas en sistemas embebidos con una FPU (Floating Point Unit) limitada.
En lugar de implementar estas aproximaciones polinómicas complejas manualmente, delegamos el trabajo en la biblioteca para aprovechar las optimizaciones específicas de la arquitectura (como instrucciones FMA de hardware). Cuando realizas operaciones que pueden salirse de los límites representables, como la raíz cuadrada de un número negativo, el estándar define dos mecanismos de error: el uso de la variable global errno (configurado con EDOM para errores de dominio o ERANGE para desbordamientos) y la propagación de valores especiales como NaN (Not a Number) o Inf (Infinito). Si ignoras estos retornos o no verificas el estado de la máquina, los valores inválidos se propagarán silenciosamente por toda tu cadena de cálculos, convirtiendo el depurador en un infierno cuando el resultado final sea un NaN sin una causa aparente. Para facilitar el manejo de tipos genéricos, existe <tgmath.h> [C99], que selecciona la función correcta según el tipo de los operandos, evitando que tengas que escribir código diferente para cada precisión.
#define _POSIX_C_SOURCE 200809L
#include <stdio.h>
#include <math.h>
#include <errno.h>
int main(void) {
double x = 3.0;
double y = 4.0;
double a = 0.1;
double b = 0.2;
/* 1. Uso de hypot para evitar desbordamientos intermedios */
// En lugar de sqrt(x*x + y*y), hypot gestiona mejor la escala
double distancia = hypot(x, y);
printf("Distancia (hypot): %.2f\n", distancia);
/* 2. atan2 para determinación de cuadrantes */
// atan(y/x) falla si x=0 o si no se conoce el signo de x y y
double angulo = atan2(y, x);
printf("Ángulo (atan2): %.4f rad\n", angulo);
/* 3. Fused Multiply-Add (FMA) para precisión y velocidad */
// Calcula (a * b) + c con un solo redondeo al final
double producto_suma = fma(a, b, 0.0000000000000001);
printf("FMA (a*b + c): %.17f\n", producto_suma);
/* 4. log1p para alta precisión con valores cercanos a cero */
// log1p(x) calcula log(1 + x) de forma mucho más precisa que log(1 + x)
double pequeño = 1e-15;
printf("log(1+x) manual: %.17f\n", log(1.0 + pequeño));
printf("log1p(x) real: %.17f\n", log1p(pequeño));
/* 5. Manejo de errores y valores especiales */
errno = 0;
double error_val = sqrt(-1.0);
if (isnan(error_val)) {
printf("Resultado es NaN (Not a Number)\n");
}
if (errno == EDOM) {
printf("Error detectado: EDOM (Dominio inválido)\n");
}
return 0;
}
Desglose del ejemplo
En la primera parte, utilizamos hypot(x, y). Si intentáramos calcular la distancia usando sqrt(x*x + y*y), y x o y fueran valores extremadamente grandes, x*x podría causar un overflow (desbordamiento) hacia el infinito, haciendo que sqrt devuelva Inf incluso si el resultado final debería ser un número representable. hypot implementa una lógica interna que escala los valores para evitar este problema.
Para la trigonometría, atan2(y, x) es la versión robusta. A diferencia de atan(y/x), atan2 analiza los signos de ambos argumentos para determinar en qué cuadrante se encuentra el vector, y maneja correctamente el caso donde x es cero, evitando una división por cero.
En el bloque de fma, estamos utilizando la instrucción Fused Multiply-Add. En un cálculo estándar como (a * b) + c, el procesador suele redondear el resultado de a * b y luego vuelve a redondear al sumar c. fma realiza la operación con una precisión extendida interna y realiza un único redondeo al final, lo que reduce el error de redondeo acumulado, algo crítico en algoritmos de álgebra lineal o filtros digitales.
Un detalle de precisión profesional es log1p(x). Cuando x es muy pequeño (como 1e-15), la expresión 1.0 + x en punto flotante pierde la información de x debido al redondeo, porque el exponente de 1.0 es mucho mayor. log1p utiliza una serie de Taylor o algoritmos aproximados que mantienen la precisión de los bits significativos de x sin necesidad de sumarlo explícitamente a 1.0.
Finalmente, el manejo de errores muestra la dualidad del estándar. Al evaluar sqrt(-1.0), la función devuelve NaN. Sin embargo, la implementación también debe configurar errno a EDOM (Domain Error). Es vital usar isnan() para verificar el valor, ya que las comparaciones normales como if (valor == NaN) siempre fallan, incluso si valor es NaN.
El error frecuente
// ERROR: Comparación incorrecta de NaN y pérdida de precisión
double v = sqrt(-1.0);
if (v == 0) { /* Nunca entrará, v es NaN */ }
// ERROR: Uso de log(1+x) para valores minúsculos
double tiny = 1e-16;
double res = log(1.0 + tiny); // res será 0.0 por pérdida de precisión
El uso de log(1.0 + x) en lugar de log1p(x) para valores muy cercanos a cero es un error clásico de precisión numérica. El resultado es un 0.0 silencioso que puede invalidar modelos matemáticos enteros. Además, intentar comparar resultados con NaN usando == es un error de lógica; siempre se debe usar isnan() o isfinite().
N° 108