Newsgroup: it.comp.lang.c++
Titolo: Radici quadrate

Buon giorno a tutto il ng
scrivo in quanto ho  un problema che mi affigge da qualche giorno:
sto cercando di implementare un metodo per calcolare la radice quadrata
di un numero reale (nel caso specifico la radice di 2 che è
irrazionale). Mi sembra ovvio che non voglio usare la funzione sqrt(2)
di cmath.h in quanto questo va contro i miei intenti. Avevo in mente di
provare con il seguente algoritmo ricorsivo:
a_{n+1} = \frac{1}{2} \cdot (a_n + \frac{x}{a_n}) dove per n che tende
all'infinito si ha proprio la radice quadrata di x.
I problemi sono questi:
(i) il codice (che posterò in fondo) funziona ma appena arriva a
\sqrt(2) questa ricomincia a crescere senza che io possa sapere il
valore esatto (o in caso di 2 una approssimazione). Volevo implementare
un metodo che controlli il valore minimo della successione e mi
restituisca quel numero.
(ii) il calcolo è troppo impegnativo e per ottenere la radice ci mette
troppo. Ho provato anche a definire una funzione inline ma niente...
per ottenere la radice ci metto sempre almeno un paio di minuti. Come
ovviare a questo problema?
Ringrazio in anticipo per eventuale risposte
Saluti Simon
--------------------------------------------
#include <iostream>
using namespace std;
int main()
{
	double MaxValue = 1000000000;
	cout << "Programma per calcolare la radice quadrata di un numero
reale: \n Immetti il numero da calcolare: " << endl;
	double a;
	cin >> a;
	double a_n1;
	double a_n = 1;
	for (int i = 0; i <= MaxValue; i++)
	{
		a_n1 = (a_n + (a / a_n))/2;
		cout << a_n1 << endl;
		a_n = a_n + 0.0000001;
	}
	return 0;
}
--------------------------------
ps. Premetto di non essere molto ferrato nel c++ quindi scusatemi per
l'eventuale codice orrido.... (è il migliore che mi è venuto in mente
e comunque sto imparando perciò siate clementi con il sottoscritto :)
).

Provo a risponderti, anche se a occhio il tuo è, più che un problema di 
programmazione, un problema di calcolo numerico. Chiedo scusa ai 
frequentatori del newsgoup per l'indentazione, ma OE non mi ha inserito i 
'>' nel testo originale, non so il perchè (O_o)
<POST ORIGINALE>
a_{n+1} = \frac{1}{2} \cdot (a_n + \frac{x}{a_n}) dove per n che tende
all'infinito si ha proprio la radice quadrata di x.
I problemi sono questi:
(i) il codice (che posterò in fondo) funziona ma appena arriva a
\sqrt(2) questa ricomincia a crescere senza che io possa sapere il
valore esatto (o in caso di 2 una approssimazione). Volevo implementare
un metodo che controlli il valore minimo della successione e mi
restituisca quel numero.
(ii) il calcolo è troppo impegnativo e per ottenere la radice ci mette
troppo. Ho provato anche a definire una funzione inline ma niente...
per ottenere la radice ci metto sempre almeno un paio di minuti. Come
ovviare a questo problema?
Ringrazio in anticipo per eventuale risposte
</POST ORIGINALE>
Quella formula ricorsiva è il metodo di newton per la ricerca degli zeri di 
una funzione.
Perchè iterando tu sia sicuro di arrivare a un risultato esatto, devi 
cautelarti rispetto al valore di partenza.
Questo non è il newsgroup adatto per queste questioni, it.scienza.matematica 
andrebbe sicuramente meglio. In breve, il tuo problema sta principalmente in 
questa riga:
double a_n = 1;
per cui fissi il valore iniziale indipendentemente dal radicando. Prova a 
mettere a_n = a (proprio il radicando), e vedi se funziona =P
La condizione di uscita poi non è il massimo, puoi fare di meglio. Provo a 
riscriverti il ciclo:
#include <cmath>   // per abs(), valore assoluto
#include <limits>  // per epsilon(), ovvero il tuo epsilon macchina
// in 'a'  hai il radicando, come nel tuo caso
a_n1 = a;
do {
a_n = a_n1;
a_n1 = (a_n + (a / a_n))/2;
}while(abs(a_n1 - a_n) > numeric_limits<double>::epsilon());
cout << "La radice è: " << a_n1 << endl;
Spero di esserti stato utile =P !
'sera lietissima, Marco 

 'sera lietissima, Marco
Saluti Simon

 Mi spiegeresti questa riga? (Non conosco limits.h :( )
Guarda, ad essere sincero non lo conosco nemmeno io troppo...nel senso che 
non so ahimè con precisione che cosa sia numeric_limits<>. Detto alla buona, 
direi che è una classe template che ti fornisce, con una certa interfaccia, 
i limiti numerici del tuo compilatore/architettura. Ad esempio, ti può 
fornire il massimo intero rappresentabile con un tipo int, che vale appunto 
"numeric_limits<int>::max()", oppure ti fornisce il numero di cifre decimali 
di un long come "numeric_limits<long>::digits10()".
In questo caso, è interessante sapere quando non ha più senso iterare, cosa 
che succede ad esempio quando due termini consecutivi della successione sono 
indistinguibili nella precisione della tua macchina. Questa precisione la 
ottieni con "numeric_limits<double>::epsilon()", che si può dire sia il 
passo di discrtizzazione dei double. Quando quindi il valore assoluto della 
differenza tra due termini consecutivi della successione è pari all'epsilon 
macchina ti puoi fermare.
Mi sono capito =P ?
'sera lietissima, Marco 

Markobre ha scritto:

Penso di si adesso cerco con gooooooogle qualcosa in più. Grazie sei
stato gentilissimo (anche perchè praticamento ho risolto tutti i
problemi che si risolvono come limiti di successioni :p )


 meglio.
Aveva due problemi : il messo a giorno (si dice così in
italiano?) della variabile, e la condizione di fine. Corrigendo
questi *due* problemi (come fà il tuo codice), e
un'inizialisazione con a o con 1 si comportono stesso.
Più la prima approssimazione è bene, meno volte passaremmo nel
ciclo. Prova:
    int exp ;
    double a_n1 = frexp( a, &exp ) ;
    a_n1 = ldexp( a_n1, exp / 2 ) ;
da essempio. (Le funzioni frexp e ldexp manipulano il campo
dell'esposente direttamente -- qui, ho dividato l'esponante per
2.)
Quando ho implementato una biblioteca standard C, 20 anni fa, hò
estrato i quattro bits superiori della mantissa e il bit
inferiore dell'esposante (quello che si perde nella divisione
intiere), e mi son serviti di loro come indizio in una tavola
per ottener il valore iniziale. Che era abbastanza giusto per
che tre cicli bastano, sempre.
(Se no il tuo codice è giusto.)
--
James Kanze                                           GABI Software
Conseils en informatique orientée objet/
                   Beratung in objektorientierter Datenverarbeitung
9 place Sémard, 78210 St.-Cyr-l'École, France, 
+33 (0)1 30 23 00 34

 italiano?) 
"aggiornamento"
-- 
MSN: sirpsycho_78 -at- hotmail dot com
"Miss, simme juti a fernì int' 'o cess, e mo' ca rimango 'i sulo, te 
mann 'a fanculo."

 2.)
Bella l'idea di partire dividendo per due l'esponente ! Così certo parti 
parecchio vicino alla soluzione, e quindi si itera meno. 

Grazie a tutti adesso sono riuscito a implementare il codice e ci mette
un battito di ali :) .
Saluti Simon

 di un long come "numeric_limits<long>::digits10()".
Giusto. E' un buon modo di astrarre da cose altrimenti aberranti.
-- 
Matteo Settenvini
FSF Associated Member
Email : 
mat...@member.fsf.org
-----BEGIN GEEK CODE BLOCK-----
Version: 3.12
GCS d--(-) s+:- a-- C++ UL+++
P?>++ L+++>$ E+>+++ W+++ N++ o?
w--- O- M++ PS++ PE- Y+>++
PGP+++ t+ 5 X- R tv-- b+++ DI+
D++ G++ e h+ r-- y?
------END GEEK CODE BLOCK------

 
Se vuoi migliorarlo ancora, guardando alle condizioni da rispettare del
metodo di Newton, se ben ricordo puoi usare il metodo di bisezione per
iniziare a convergere sotto una certa soglia (specie con funzioni con
più zeri), e poi applicare Newton iterativamente.
Reminescenze vaghe...
Ciao,
