È sicuro eseguire la cascata di hypot()?

Question

Vorrei calcolare la norma (lunghezza) di vettori a tre e quattro dimensioni. Sto usando numeri in virgola mobile a precisione doppia e voglio fare attenzione a evitare inutili overflow o underflow.

La libreria di matematica C fornisce hypot(x,y) per calcolare la norma di vettori bidimensionali, facendo attenzione a evitare il trabocco/trabocco nei calcoli intermedi.

La mia domanda: è sicuro utilizzare hypot(x, hypot(y, z)) e hypot(hypot(w, x), hypot(y, z)) per calcolare le lunghezze dei vettori a tre e quattro dimensioni, rispettivamente?

Answer 1

E 'sicuro, ma è uno spreco: è sufficiente calcolare sqrt() una volta, ma quando in cascata hypot(), si chiamerà sqrt() per ogni chiamata a hypot(). Normalmente potrei non essere preoccupato per le prestazioni, ma questo può anche degradare la precisione del risultato. Potresti scrivere il tuo:

double hypot3(double x, double y, double z) { 
    return sqrt(x*x + y*y + z*z); 
} 

ecc. Questo sarà più veloce e più accurato. Non credo che qualcuno sarebbe confuso quando vedono hypot3() nel tuo codice.

La libreria standard hypot() potrebbe disporre di trucchi per evitare l'overflow, ma potrebbe non esserne preoccupato. In genere, hypot() è più preciso di sqrt(x*x + y*y). Vedi e_hypot.c nel codice sorgente GLibC.

Answer 2

È sicuro (quasi) utilizzare hypot(x, hypot(y, z)) e hypot(hypot(w, x), hypot(y, z)) per calcolare le lunghezze dei vettori a tre e quattro dimensioni.

C non con forza specifica che hypot()deve lavoro per un double x, y che hanno una limitata double risposta. Ha weasel words di "senza eccessivo overflow o underflow".

Tuttavia, dato che hypot(x, y) funziona, un'implementazione ragionevole hypot() eseguirà hypot(hypot(w, x), hypot(y, z)) come necessario. C'è solo 1 incremento (alla fine bassa)/decremento (alla fine alta) dell'intervallo esponenziale binario perso quando con 4-D contro 2-D.

Per quanto riguarda velocità, precisione e intervallo, profilo codice contro sqrtl((long double) w*w + (long double) x*x + (long double) y*y + (long double) z*z) in alternativa, ma sembra necessario solo con determinati obiettivi di codifica.

Answer 3

Ho fatto alcuni esperimenti con questo genere di cose. In particolare ho esaminato un'implementazione semplice, un'implementazione usando gli hypot e (una traduzione C della versione di riferimento di) la funzione BLAS DNRM2.

Ho scoperto che per quanto riguarda l'over e il underflow, le implementazioni BLAS e hypot erano le stesse (nei miei test) e di gran lunga superiori all'implementazione normale. Per quanto riguarda il tempo, per i vettori di dimensioni elevate (centinaia), il BLAS era circa 6 volte più lento del piano, mentre l'ipotetico era 3 volte più lento del BLAS. Le differenze di tempo erano un po 'più piccole per le dimensioni più piccole.

Answer 4

Se il codice non è in grado di utilizzare hypot() né tipi di precisione più ampi, un metodo lento esamina gli esponenti utilizzando frexp() e ridimensiona gli argumnet @greggo.

#include <math.h> 

double nibot_norm(double w, double x, double y, double z) { 
    // Sort the values by some means 
    if (fabs(x) < fabs(w)) return nibot_norm(x, w, y, z); 
    if (fabs(y) < fabs(x)) return nibot_norm(w, y, x, z); 
    if (fabs(z) < fabs(y)) return nibot_norm(w, x, z, y); 
    if (z == 0.0) return 0.0; // all zero case 

    // Scale z to exponent half-way 1.0 to MAX_DOUBLE/4 
    // and w,x,y the same amount 
    int maxi; 
    frexp(DBL_MAX, &maxi); 
    int zi; 
    frexp(z, &zi); 
    int pow2scale = (maxi/2 - 2) - zi; 

    // NO precision loss expected so far. 
    // except w,x,y may become 0.0 if _far_ less than z 
    w = ldexp(w, pow2scale); 
    x = ldexp(x, pow2scale); 
    y = ldexp(y, pow2scale); 
    z = ldexp(z, pow2scale); 

    // All finite values in range of squaring except for values 
    // greatly insignificant to z (e.g. |z| > |x|*1e300) 
    double norm = sqrt(((w * w + x * x) + y * y) + z * z); 

    // Restore scale 
    return ldexp(norm, -pow2scale); 
} 

codice di prova

#include <float.h> 
#include <stdio.h> 
#ifndef DBL_TRUE_MIN 
#define DBL_TRUE_MIN DBL_MIN*DBL_EPSILON 
#endif 

void nibot_norm_test(double w, double x, double y, double z, double expect) { 
    static int dig = DBL_DECIMAL_DIG - 1; 
    printf("  w:%.*e x:%.*e y:%.*e z:%.*e\n", dig, w, dig, x, dig, y, dig, z); 
    double norm = nibot_norm(w, x, y, z); 
    printf("expect:%.*e\n", dig, expect); 
    printf("actual:%.*e\n", dig, norm); 
    if (expect != norm) puts("Different"); 
} 

int main(void) { 
    nibot_norm_test(0, 0, 0, 0, 0); 
    nibot_norm_test(10/7., 4/7., 2/7., 1/7., 11/7.); 
    nibot_norm_test(DBL_MAX, 0, 0, 0, DBL_MAX); 
    nibot_norm_test(DBL_MAX/2, DBL_MAX/2, DBL_MAX/2, DBL_MAX/2, DBL_MAX); 
    nibot_norm_test(DBL_TRUE_MIN, 0, 0, 0, DBL_TRUE_MIN); 
    nibot_norm_test(DBL_TRUE_MIN, DBL_TRUE_MIN, DBL_TRUE_MIN, 
     DBL_TRUE_MIN, DBL_TRUE_MIN * 2); 
    return 0; 
} 

Risultati

 w:0.00000000000000000e+00 x:0.00000000000000000e+00 y:0.00000000000000000e+00 z:0.00000000000000000e+00 
expect:0.00000000000000000e+00 
actual:0.00000000000000000e+00 
    w:1.42857142857142860e+00 x:5.71428571428571397e-01 y:2.85714285714285698e-01 z:1.42857142857142849e-01 
expect:1.57142857142857140e+00 
actual:1.57142857142857140e+00 
    w:1.79769313486231571e+308 x:0.00000000000000000e+00 y:0.00000000000000000e+00 z:0.00000000000000000e+00 
expect:1.79769313486231571e+308 
actual:1.79769313486231571e+308 
    w:8.98846567431157854e+307 x:8.98846567431157854e+307 y:8.98846567431157854e+307 z:8.98846567431157854e+307 
expect:1.79769313486231571e+308 
actual:1.79769313486231571e+308 
    w:4.94065645841246544e-324 x:0.00000000000000000e+00 y:0.00000000000000000e+00 z:0.00000000000000000e+00 
expect:4.94065645841246544e-324 
actual:4.94065645841246544e-324 
    w:4.94065645841246544e-324 x:4.94065645841246544e-324 y:4.94065645841246544e-324 z:4.94065645841246544e-324 
expect:9.88131291682493088e-324 
actual:9.88131291682493088e-324