Računanje kvadratnog korena u C++ na veliki broj decimala (bilo: Kako na veliko u C++)

Računanje kvadratnog korena u C++ na veliki broj decimala (bilo: Kako na veliko u C++) 15.01.2005. u 14:49 · pre 265 meseci #1

Ljudi ne mogu izračunati kvadratni korijen iz 2 na 40 decimalnih mjesta u c++.
Pokušao sam sa algoritmom za ručno računanje al ne znam kako
da ga implementiram u c++.Naime potrebno mi je da nekako brojeve
predstavim veće od double. Evo kako je to izvedeno algoritmom:

Code:

1.Izaberi grubu približnu vrijednost B od  sqrt(A).
  B je predpostavljena neka gruba vrednost korena
  A je broj iz kojeg se traži kvadratni koren

2.Podeli A sa G pa saberi to sa G i to sve podeli sa 2, tj:

  C = ((A/G)+G)/2

3.Ukoliko je C dovoljno prezican, stani. Inače , ako B = C vrati se
  na korak 2.

Evo kako bi to bilo za kvadratni koren iz 2.
Da izračunamo (2), izaberimo da je B = 1.5.

C = (2/1.5 + 1.5)/2 = 1.41666666666
C = (2/1.41666666666+1.41666666666)/2=1.41421568628
C = (2/1.41421568628+1.41421568628)/2=1.41421356238
C = (2/1.41421356238+1.41421356238)/2=1.41421356238

U c++ to bi izgledalo ovako:

Code:

//#################################### INCLUDES ###############################
#ifndef  iostream.h
   #include <iostream.h>
#endif
//ifndef string.h
   //#include <string.h>
//#endif
//#ifndef conio.h
   //#include <conio.h>
//#endif
//#ifndef ctype.h
   //#include <ctype.h>
//#endif
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

//##################### PROTOTYPES ###################################

    double findSqrt ( double number , double approximation  ,  int placesOfAccuracy  );

//##########################################################

//*************************************** Main *****************************************
void main()
{

   cout << "Enter the number of which to take a root : ";

   double number;

   cin >> number;
   cout << endl << "Enter your estimate : " ;

   double estimate;

   cin >> estimate;

   cout << "Enter the digits of accuracy required : ";

   int accuracy;

   cin >> accuracy;

   double result = findSqrt ( number , estimate , accuracy );

   cout << "The result is : " << result;

   cout << endl << "goodbye.............." << endl;

}
/*******************************************************************************
Function Name: findSqrt
Parameters:
Description:
Returns:
Comments:

*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*/
double findSqrt (  double number , double approximation  ,  int placesOfAccuracy  )
{

   placesOfAccuracy  /= 2;

   while ( placesOfAccuracy > 0 )
   {

       approximation = ( ( number / approximation ) + approximation )  /  2 ;
       placesOfAccuracy --;
   }
   return approximation;
}

No meni treba na 40 decimala.
Na webu sam pronasao približno resenje al sa još nekim funkcijama koje bi trebalo
izbrisati a i program je totalno kompliciran. Nalazi se na:
http://www.people.memphis.edu/~wjermann/download.htm
Potrebna su ova tri filea:
number.h
number.cpp
numtest.cpp
Kad bi mogao to sve obrisat i staviti u jedan file ili slično jer siguran sam da veliki dio nije potreban. Ja sam skratio nummtest.cpp i izgleda ovako

Code:

#include "number.cpp"

main()

{ Number y;

_DP = 40;

cout << "sqrt(2) is:\n\n";

y = sqrt("2"); y.Numpnt();

return 0;

}

ali problem je što ne znam koko je do toga došlo.

AKO NEKO MOŽE DA MI POMOGNE BIO BI VEOMA ZAHVALAN.
UNAPRED HVALA!!!

Re: Kako na veliko u C++ 15.01.2005. u 19:28 · pre 265 meseci #2

Pretpostavimo da ne smes da koristis 3rd party biblioteke.
Zaboravi Njutnov iterativni obrazac. Deljenje se moze izbeci, ali ti ostaje realizacija mnozenja dva broja sa velikim brojem decimala.
Mozes upotrebiti sumu koja se dobija Tejlorovim razvojem (1+x)^0.5. Ova suma jako sporo konvergira. Postoje i sume koje brze konvergiraju. Mozes upotrebiti npr. ovu:

Code:
sqrt(2) = 7 / 5 * ( 1 + 1 / 100 + ((1 * 3) / (100*200)) + ((1 * 3 * 5) / (100 * 200 * 300)) + ...)

Preostaje realizacija potrebnih fixed point operacija. Brojeve predstavis kao nizove long varijabli. Npr.

Code:
long suma[15], sabirak[15];

Nulti element niza nek cuva celobrojnu vrednost. Prvi element niza neka sadrzi celobrojnu vrednost koja odgovara prvim 4 decimalama broja, itd... Znaci, decimalna tacka je izmedju nultog i prvog elementa.

Treba da realizujes sledece operacije:
-deljenje broja celobrojnom konstantom
-mnozenje broja celobrojnom konstantom
-sabiranje dva broja

Izracunavanje radis po algoritmu:

Code:
suma = 1 0000 0000 0000 0000 0000 0000 ...  = 1.00000... (prazan prostor odvaja elemente niza)
sabirak = 1 0000 0000 0000 0000 0000 0000...  = 1.00000... (takodje)

n = 1.. potreban broj koraka
   pomnozi sabirak sa (2 * n - 1)
   podeli sabirak sa (100 * n)
   dodaj sabirak sumi
pomnozi sumu sa 7
podeli sumu sa 5

Deljenje radis po algoritmu:

Code:
temp = 0;
i = 0 .. zadnji element niza
   temp = 10000L * temp + broj[i]
   broj[i] = temp / konstanta kojom delimo broj
   temp %= konstanta

Za sabiranje i mnozenje se malo pomuci :) Hint: Ne racuna se od 0 .. zadnji element niza, vec obrnuto
Prilikom izracunavanja moras paziti da rezultat ne predje maksimalne dozvoljene vrednosti za long.

E sad, da ti ja ne bih uradio ceo zadatak, malo se pomuci. Sta ti nije jasno pitaj.

Potpis

O_o

Re: Kako na veliko u C++ 16.01.2005. u 08:08 · pre 265 meseci #3

Ovo je, nadam se da ce se gospoda moderatori sloziti, tema za "The Art of Computer Programming" diskusionu grupu...

Potpis

Dejan Lekic
software engineer, MySQL/PgSQL DBA, sysadmin

Re: Računanje kvadratnog korena u C++ na veliki broj decimala (bilo: Kako na veliko u C++) 17.01.2005. u 23:11 · pre 265 meseci #4

Autor teme nije pokazao interesovanje da nesto nauci, vec je posegnuo za gotovim resenjima, te je nasao jedno tipa "the shortest c program for ...", koje na prvi pogled, a Boga mi i na drugi nije razumljivo pocetniku. Zasto? No pain no gain. Uz malo truda zadatak se moze lako resiti. Resenje koje sledi je nadam se razumljivo svakom ko poznaje osnove programiranja.

Code:

//sqrt(2) = 7 / 5 * ( 1 + 1 / 100 + ((1 * 3) / (100 * 200)) + ((1 * 3 * 5) / (100 * 200 * 300)) + ...)
#include <iostream>
using namespace std;
const int VEL = 3337; // izracunava se (VEL - 4) * 3 decimala - znaci, izracunavamo 9999 decimala
const unsigned long D = 1000UL;
unsigned long suma[VEL] = {1UL, 400UL}; // = 7 / 5 = 1.4000000...
unsigned long sabirak[VEL] = {1UL, 400UL}; // = 7 / 5 = 1.4000000...
int main(void){
    int k = 0;
    for(unsigned long n = 1; k < VEL; n++){ // prvi sabirak je vec dodat; dodajemo potreban broj sabiraka
        unsigned long temp = 0;
        for(int i = k; i < VEL; i++){ // sabirak = sabirak / (100 * n)
            temp = D * temp + sabirak[i];
            sabirak[i] = temp / (100UL * n);
            temp %= 100UL * n;
        }
        unsigned long temp2 = temp = 0;
        for(int i = VEL - 1; i >= k; i--){
            sabirak[i] = sabirak[i] * (2UL * n - 1UL) + temp; // sabirak = sabirak * (2 * n - 1)
            temp = sabirak[i] / D;
            sabirak[i] %= D;
            suma[i] += sabirak[i] + temp2; // suma = suma + sabirak
            temp2 = suma[i] / D;
            suma[i] %= D;
        }
        for( ; (sabirak[k] == 0) && (k < VEL); k++); // skracujemo izracunavanja u sledecim ciklusima
        if(k && (k < VEL)) k--;
    }
    cout << "sqrt(2) = 1."; // ispis
    for(int i = 1; i < VEL - 3; i++){
        cout.width(3);
        cout.fill('0');
        cout << suma[i];
    }
    return 0;
}

Potpis

O_o

Re: Računanje kvadratnog korena u C++ na veliki broj decimala (bilo: Kako na veliko u C++) 07.03.2005. u 22:01 · pre 263 meseca #5

Njutnovu iteraciju za računanje kvadratnog korena zaboravi jer je previše spora (deljenje u svakoj iteraciji, često se ono samo računa iteracijom). Probaj da kreneš od iteracije za recipročni koren, a finalni rezultat dobijaš kad pomnožiš vrednost iz iteracije sa x što unosi izvesnu grešku, ali najviše 1 cifra. VC tretira long double kao double, ali mislim da postoji mogućnost da se ovo isključi i da long double bude 80 bitova što je striktno vezano za x86. Sumiranje redova, naročito sporo konvergentnih, zaboravi jer će greška zaokruživanja upropastiti rezultat. Eventualno da probaš da primeniš neku transformaciju kao Aitkenovu/Kumerovu da malo ubrzaš celu stvar. Bitna stvar je da ne možeš dobiti veću preciznost od one koju ti nude fp brojevi, dakle ako imaš 16 tačnih cifara tvoj rezultat ne može imati više, samo može manje. Druga mogućnost je da upotrebiš FFT i da sve radiš pomoću toga (vidi Numerical Recipies, imaš tamo detaljno). Stvar radi prilično brzo, možeš koristiti Intel/AMD biblioteke za sam FFT. Treća varijanta, koju toplo NE preporučujem je da napišeš svoju klasu za rad sa fp brojevima u kojoj ćeš moći mirno da staviš za mantisu i eksponent koliko želiš. Samo cela ta skalamerija ima da ide prilično sporo, plus što ćeš morati da prostudiraš IEEE standard (zaokruživanje, izuzeci, normalizovani i denormalizovani brojevi), a ako poželiš sin, cos, exp, ln sa sličnom tačnošću cela stvar će postati noćna mora, numerika ipak nije baš tako prosta kao što izgleda.

Rezimirano, ako želiš iole veći broj decimala druga varijanta (FFT) je ono na šta treba da ideš sa tim da koren ražunaš kao (1/sqrt(x))*x Prvo probaj sa floatom, ako te to ne zadovolji idi na double, mada ćeš tada imati mnogo više tačnih cifara tako da ćeš mnoge stvari moći da uradiš naivno a da one ipak rade, eventualno ćeš to platiti u brzini, ali ne i u tačnosti.

Potpis

One word, one instruction, one cycle. RISC.

Re: Računanje kvadratnog korena u C++ na veliki broj decimala (bilo: Kako na veliko u C++) 17.02.2011. u 12:30 · pre 191 mesec #6

Ne, korenovanje nije sporo Njutnovim metodom tangente, ako se realizuje efikasno preko brzog množenja.

Pretpostavka je da je implementirano brzo množenje pomoću FFT. Funkciju $f_n(t)=\frac 1{\sqrt[n]{t}}$ možemo računati Njutnovim metodom tangente na sledeći način:

$x=\frac 1{\sqrt[n]{a}}\quad\Leftrightarrow\quad \frac 1{x^n}-a=0$ ,

$x_{k+1}=x_k\left(1+\frac{1-ax_k^n}n\right)$ .

Dakle, za $n=1$ dobijamo algoritam delenja preko množenja i sabiranja

$x_{k+1}=x_k(2-ax_k)$ ,

a za $n=2$ algoritam računanja recipročne vrednosti kvadratnog korena, takođe preko množenja i sabiranja

$x_{k+1}=x_k\left(1+\frac{1-ax_k^2}2\right)$ .

Ako je množenje realizovano pomoću FFT sa složenošću $O(n\log(n))$ , onda su ovakvo delenje i računanje recipročne vrednosti kvadratnog korena (i uopšte $t^{-1/n}$ ) računaju u vremenu $O(n\log^2(n))$ . No onda se kvadratni koren dobija još jednim množenjem po formuli $\sqrt t=t\,\frac 1{\sqrt{t}}$ , dakle sa ukupnim vremenom $O(n\log^2(n))$ .

Potpis

Nije bitno koji su zaključci izvučeni, već kako se do njih došlo.

Forum kao aplikacija na telefonu

Srodne teme