Newsgroup: it.comp.lang.c++
Titolo: Problemi elevamento a potenza

Ciao.
 Oggi mi sono scontrato con il seguente problema:
Calcolando queste due espressioni
exp(j*3/0.45) - exp(j*3/0.7)^(0.7/0.45)
exp(j*3/0.45) - exp(j*3/7)^(7/0.45)
la prima risulta essere diversa da 0 mentre la seconda fa 0
Cosa ancora piu' strana, lo stesso risultato lo otteno anche in Matlab
e Fortran.
Alla fine ho fatto un loop per verificare dove avevo problemi
Di seguito riporto il codice (C++):
----------------------------------------start-------------------
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <complex>
using namespace std;
main()
{
   complex<double>  a,b1,b2,b3,b4,c,j;
   j = complex<double>(0,1);
   a = exp(j*3./0.45);
   for(double x=-1.; x<=1. ; x+=0.0001)
   {
       b1 = pow(exp(j*3./x ),  x/0.45);
       c = a-b1;
       printf("%10.5lf %10.5lf %10.5lf\n", x,c.real(), c.imag());
   }
}
-----------------------------------stop------------------
Stesso risultato ottenuto su
cygwin, linux debian, Microsoft c++.
Possibile che sia tutto sbagliato?
Forse sto sbagliando io!!!!!!
Graditi suggerimenti e consigli.

 la prima risulta essere diversa da 0 mentre la seconda fa 0
Non ho letto bene il tuo problema, ma una cosa te la dico "a prescindere":
Cosa intendi per "diversa da 0" e per "fa 0" ?
Con i numeri in virgola mobile l'"ugualianga" � un concetto che nei casi 
pratici ha spesso bisogno di essere "elasticizzato": di solito anzich� 
espressioni si controlla se i due valori sono abbastanza "vicini".
Es.
double a, b;
if(a==b) //ugualianza perfetta. pu� soffrire di problemi di 
calcolo/rappresenzazione
if(abs(a-b)<=epsilon) //"vicinanza". nei casi pratici si usa questa
Vedere ad esempio
http://www.cygnus-software.com/papers/comparingfloats/Comparing%20floating%20point%20numbers.htm
HTH
Beppe

 Beppe
Provo ad essere piu' chiaro.
Riminiscenze delle scuole medie:
" La potenza di una potenze e' una potenza che ha per base la stessa
base e per esponente il prodotto degli esponenti "
In modo piu' stringato (concedimi il simbolo ^ per indicare
l'elevamento a potenza)
(x^a)^b = x^(a*b)
Detto cio' , nel codice che ho riportato ho fatto la potenza di una
potenza.
I due termini sono sostanzialmente uguali (scritti in forme diverse).
So perfettamente che lo 0 con i floating point e' un concetto relativo
e se avessi ottenuto come risultato
10^-7 non avrei scritto questo messaggio.
Per comidita' riporto alcune uscite che produce il codice
 -0.91200   -0.04941    0.15990
 -0.91190   -0.04911    0.15853
 -0.91180   -0.04881    0.15717   questi sono numeri sufficientemente
 -0.91170   -0.04850    0.15581   diversi da zero da farmi sorgere il
 -0.91160   -0.04820    0.15444   sospetto che ci sia qualcosa che non
va?
 -0.91150   -0.04789    0.15308
 -0.91140   -0.04758    0.15172
 -0.91130   -0.04727    0.15036
   0.99670    0.00000   -0.00000
   0.99680    0.00000    0.00000
   0.99690    0.00000   -0.00000
   0.99700    0.00000    0.00000
   0.99710    0.00000    0.00000
   0.99720    0.00000   -0.00000
il primo gruppo (secondo me) non funziona il secondo si.
Spero di essere stato piu' chiaro

 
...
Attento che questo � calcolo numerico e non centra na mazza con le 
scuole medie :)
Ad esempio uno dei possibili problemi nel tuo codice (non � detto sia 
l'unico e quello risolutivo) potrebbe essere il fatto che usi spesso 
valori interi. Purtroppo in c++ l'aritmetica intera (ma neanche quella 
fp) � quella delle scuole medie :)
Ad esempio:
exp(j*3/7)
porva a mettere un .0 dopo tutti gli interi
Ad Es.
exp(j*3.0/7.0)

exp(j*3.0/7.0)
double x=-1.; x<=1. ; x+=0.0001)
   {
       b1 = pow(exp(j*3./x ),  x/0.45);
       c = a-b1;
Si vedono abbastanza punti decimali ?
la x e' double;
x parte da -1. ed arriva ad 1. e si incrementa di 0.0001
poi c'e' 3./x ed x/0.45
Altri suggerimenti?

> Altri suggerimenti?
Avendo dovuto risolvere in passato problemi molto complessi mi sono 
riscontrato pure io in questo genere di effetti collateriali e ti 
spiego il motivo.
Semplicemente devi tenere conto che i problemi di arrotondamento dovuti 
alla limitazione di 8 byte vengono amplificati se inseriti in 
particolari contesti.
Ti faccio un esempio pratico.
Siamo tutti daccordo che due valori tipo
double A=0.00000000000001;
double B=0.000000000000000000000000000000001;
informaticamente parlando sono praticamente identici.
Ma se tu scrivi la formula
double C = A / B
Non e' che ottieni circa 1 ... ottieni 1E+19. Che e' un numero 
decisamente diverso.
Nella tua formula c'e' la potenza di una potenza e quindi per sua 
natura amplifica gli errori ma non solo, ci sono divisioni di numeri 
reali che producono risultati sicuramente con piu' di 15 cifre decimali 
previste dal tipo double (internamente M$ ne archivia 17, provare per 
credere).
Capisci che questi errori di arrotondamento inevitabilmente ci sono.
Quindi scrivere
j * 3 / 0.7
produrra un certo risultato con un certo errore.
Scrivere
j * 3 / 7
produrra un altro risultato certo ma sopratutto con un errore di 
approssimazione diverso, dell'ordine di 1 cifra decimale.
Ecco quindi che l'errore di + o - una cifra decimale sulla 15esima 
cifra te lo ritrovi giustamente sul risultato finale, peccato che non 
e' piu' sulla 15esima cifra ma sulla prima...
Spero di essere stato chiaro.
Comunque non ci sono soluzioni... io ho pravo a sbattermi ad utilizzare 
librerie per aumentare la precisione (dalla banalissima doubledouble a 
quelle che ti garantiscono numeri di decimali arbitrari) ma non puoi 
comunque farci nulla perche' come ti dicevo dipende dalla formula.

 j * 3 / 7
Prima di tutto ti ringrazio per la lunga mail nella quale mi hai
esaurientemente e professionalmente esposto la tua opinione.
Non ho certezze per dirti che hai torto o ragione ma sono comunque
perplesso.
La questione della precisione non mi convince tanto.
Tanto per provare ho fatto un altro tentativo che riporto:
------------------------------start-------------------
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <complex>
using namespace std;
main()
{
   complex<long double>   a;
   complex<long double>   b1,b2,b3,b4;
   complex<long double>   j;
   long double x = .1L;
   j = complex<long double>(0,1);
   a = exp(j);
   b1 = exp(j/x);
   b2 = pow(b1,  x);
   b3 = pow(exp(j/x ),  x);
   b4 = a - b3;
   printf("x                              = +0.1L\n");
   printf("a  = exp(j)                 = %+30.25lle %+30.25lle
\n",a.real(), a.imag());
   printf("b1 = exp(j/x)             = %+30.25lle %+30.25lle
\n",b1.real(), b1.imag());
   printf("b2 = pow(b2, x)        = %+30.25lle %+30.25lle \n",b2.real
(), b2.imag());
   printf("b3 = pow(exp(j/x), x) = %+30.25lle %+30.25lle \n",b3.real
(), b3.imag());
   printf("b4 = a - b3               = %+30.25lle %+30.25lle
\n",b4.real(), b4.imag());
}
----------------------------stop---------------
Il tipo long double risponde al sizeof con 12 il che equivale a
parlare di 96 bit.
Il risultato e' il seguente:
x                             = +0.1L
a  = exp(j)                = +5.4030230586813971741400736e-01
+8.4147098480789650666459081e-01
b1 = exp(j/x)             = -8.3907152907645245226058231e-01
-5.4402111088936981343191709e-01
b2 = pow(b2, x)        = +9.6724905798807620156166964e-01
-2.5382919418612824645313519e-01
b3 = pow(exp(j/x), x) = +9.6724905798807620156166964e-01
-2.5382919418612824645313519e-01
b4 = a - b3               = -4.2694675211993648414766228e-01
+1.0953001789940247530635159e+00
Non e' che il problema e' nell'implementazione della funzione pow?
Comunque come ho gia' detto grazie in ogni caso.

> Non e' che il problema e' nell'implementazione della funzione pow?
Lo escluderei a priori perche' gli stessi problemi (come da te 
confermato) li puoi riscontare in altri linguaggi.


Non e' che tutti si appoggiano alla medesima libreria che ha qualche
problema (libc)?

 problema (libc)?
Non ne ho idea ma come ti dicevo ho avuto modo di risontrarlo anche in 
altre situazioni ed ha comunque una sua spiegazione logica.
Ad esempio
double A = 7.0;
double B = 0.45;
double C = A/B;
double D = B/A;
double E = 1.0 / D;
double F = C - E; // (C==E) -> qundi F deve essere 0
double G = log(1E-15);
double H = log(1E-15 + F );
double I = G - H; //Essendo F=0, G==H e quindi I deve essere 0
E invece I quanto vale?
Durante la strada si sono verificati arrotondamenti che la log ha messo 
in evidenza. La stessa cosa succede con la pow.
N.B. Se poi invece di 7/0.45 usi altri valori otterrai errori diversi.

 la prima risulta essere diversa da 0 mentre la seconda fa 0
Come ti dicevo sull'altro NG, non � garantito che (x^y)^z = x^(yz) tra i 
numeri complessi. Questo dipende dal fatto che l'esponenziale complesso 
non � una funzione invertibile. La stessa cosa succede con i numeri 
negativi per i quali non � vero che (x^2)^(1/2) = x.
E.