Newsgroup: it.comp.lang.c++
Titolo: [subset - OT ?] analisi numerica...

Sono un po' arrugginito (si fa per dire. In realta' sono
ipermegastrasuperarrugginito: un ferrovecchio, praticamente da
buttare) in fatto di analisi numerica. Il pezzo di codice 
alla fine del posto l'ho estratto dalla libreria submit che sembra
contenere varie routine interessanti (almeno interessanti secondo me)
per calcoli numerici di vario tipo, nonche' generazione di alcune cose
che andavo cercando da tempo (tipo combinazioni, permutazioni, ecc).
Il problema e' che la libreria contiene le due funzioni seguenti
r8_epsilon e r8_agm.
Sul mio notebook (HP nx7400 con processore Dual Core2 ) il programma 
"si pianta" nel senso che il ciclo indicato come "ciclo incriminato"
non esce mai e non esce mai perche' in pratica si arriva ad un
punto in cui a2 e b2 valgono sempre gli stessi valori quando viene
eseguita if(fabs(a2-b2) <= tol)
che risulta essere sempre falsa.
Infatti a2, b2 assumono valori rispettivamente pari a
2.6040081905309407
e
2.6040081905309402
e a2-b2 vale sempre all'incirca 4.44089e-16
mentre la precisione restituita da e8_epsilon vale circa
1.0842e-19
e per farlo funzionare devo sostituire la
riga 
return 2.0 * r; 
in r8_epsilon 
con
return 10000.0 * r;
Qualcuno si ricorda qual e' un modo corretto di
scrivere il codice sotto riportato, senza che questo dipenda
dalla floating point unit della macchina su cui il
codice andra' a girare ? 
E se qualcuno mi dice che sto cercando di farmi aiutare
a risolvere i compiti di calcolo numerico, vorrei ricordare che
ho compiuto 40 anni lo scorso aprile e che l'esame
di calcolo numerico l'ho superato quai 20 anni fa; semmai, il problema
e' che in questi 20 anni di problemi legati alla precisione
finita del calcolatore ne ho affrontati pochi e l'ultima
volta (e forse unica volta dopo l'esame) e' stata circa 10 anni fa.
 
#include <stdio.h> 
#include <math.h> 
double r8_epsilon ( void )
//****************************************************************************80
//
//  Purpose:
//
//    R8_EPSILON returns the R8 round off unit.
//
//  Discussion:
//
//    R8_EPSILON is a number R which is a power of 2 with the property that,
//    to the precision of the computer's arithmetic,
//      1 < 1 + R
//    but 
//      1 = ( 1 + R / 2 )
//
//  Modified:
//
//    06 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double R8_EPSILON, the round-off unit.
//
{
  double r;
  r = 1.0;
  while ( 1.0 < ( double ) ( 1.0 + r )  )
  {
    r = r / 2.0;
  }
  //return ( 10000.0 * r );
   return ( 2.0 * r ); // se non
}
//****************************************************************************80
double r8_agm ( double a, double b )
//****************************************************************************80
//
//  Purpose:
//
//    R8_AGM finds the arithmetic-geometric mean of two R8's.
//
//  Discussion:
//
//    The AGM of (A,B) is produced by the following iteration:
//
//      (A,B) -> ( (A+B)/2, SQRT(A*B) ).
//
//    The sequence of successive values of (A,B) quickly converge to the
AGM.
//
//  Modified:
//
//    06 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double A, B, the numbers whose AGM is desired.  A and B should
//    both be non-negative.
//
//    Output, double R8_AGM, the AGM of the two numbers.
//
{
  double a1;
  double a2;
  double b1;
  double b2;
  double tol;
  if ( a < 0.0 )
  {
    return -1.0;
  }
  if ( b < 0.0 )
  {
    return -1.0;
  }
  if ( a == 0.0 || b == 0.0 )
  {
    return 0.0;
  }
  if ( a == b )
  {
    return a;
  }
  tol = r8_epsilon ( ) * ( a + b + 1.0 );
  a1 = a;
  b1 = b;
  for ( ; ; )
  {
    a2 = ( a1 + b1 ) / 2.0;
    b2 = sqrt ( a1 * b1 );
printf("%g\n",fabs(a2-b2));
    if ( fabs ( a2 - b2 ) <= tol ) //ciclo incriminato
    {
      return ( ( a2 + b2 ) / 2.0 );
    }
    a1 = a2;
    b1 = b2;
  }
}
int main()
{
  int a=1;
  int b= 5;
 
   r8_agm(a,b);
  return 0;
}

 1.0842e-19
Cio' e' strano: normalmente su architetture ia32 (ma credo sia cosi' 
anche nel caso tu utilizzi un sistema operativo a 64 bit) il tipo double 
e' 8 byte e 15 cifre, per cui non si spiega quel e-19 e il fatto che 
moltiplicando per 10000 si risolva mi sembra molto indicativo.
Mi viene da pensare che il tuo compilatore converta le variabili di 64 
bit in 80 bit (la lunghezza di parola dello stack della FPU), esegua i 
calcoli, ed effettui alla fine l'operazione contraria.
Su gcc, ad esempio, vedi le opzioni -ffloat-store (e successive).
Su VC++ non ricordo l'esatta opzione, ma viene battezzata "improve 
floating point precision" o giu' di li'.
  bye
	av.

 floating point precision" o giu' di li'.
 
VC++ 6.0 pare funzionare gia' con le opzioni di default. Ora dovrei
scoprire quali opzioni utilizzare con g++ per vedere
se riesco a farlo funzionare (la sola -ffloat-store pare non essere 
sufficiente)
Ciao e grazie