find as many digits of the square root of 2 as possible

Question

#include <iostream>
#include <cmath>

using namespace std;

int main()
{
 double a = sqrt(2);
 cout << a << endl;
}

hi this is the program to find sqrt of 2 it prints just 1.41421 in the output how to implement it in way such that it will print 200000 digits after decimal point

1.41421..........upto 200 000 digits

Is there any approach to print like this?

Isidraisidro

A double will give you about 16 significant digits. If you want 200000, you need an arbitrary precision library (GMP for example). – Spoon

en.wikipedia.org/wiki/Methods_of_computing_square_roots – Rosenda

I don't think using the built-in sqrt function will do what you want. I'm guessing this is a homework assignment. – Glycerinate

For floating-point, s/GMP/MPFR/. – Mangosteen

See #823234 – Saire

@Glycerinate you guesses right , anything that does not make sense is always a homework/assignment !! – Seko

Unfortunately, we don't do people's homework for them on this site. – Glycerinate

You'll need to re-implement the sqrt() function yourself, as the standard one only returns up to a long double, or you'll need a library to do it. If you just need the answer, go to sagenb.org, create an account, and enter in something like: N(sqrt(2),2000) – Casino

Assuming the assignment precludes the use of libraries like gmp, you might look into implementing a square root as a continued fraction expansion (see en.wikipedia.org/wiki/Methods_of_computing_square_roots for a description). You'll still need a multi-precision integer however to contain that many digits. – Gemmulation

Answer 1

Here is the code for your question which using GNU GMP library. The code will print 1000 digits of sqrt(2), increase the number in the lines with comment to satisfy your request.

#include <stdio.h>
#include <gmp.h>

int main(int argc, char *argv[])
{
  mpf_t res,two;
  mpf_set_default_prec(1000000); // Increase this number.
  mpf_init(res);
  mpf_init(two);
  mpf_set_str(two, "2", 10);
  mpf_sqrt (res, two); 
  gmp_printf("%.1000Ff\n\n", res); // increase this number.
  return 0;
}

Please compile it with the following command:

$gcc gmp.c  -lgmp -lm -O0 -g3

Theological

Oh my god using your solution I increased the default prec to 10 million and I could print 1 million digits correctly :) – Microspore

Answer 2

It can be shown that

sqrt(2) = (239/169)*1/sqrt(1-1/57122)

And 1/sqrt(1-1/57122) can be computed efficiently using the Taylor series expansion:

1/sqrt(1-x) = 1 + (1/2)x + (1.3)/(2.4)x^2 + (1.3.5)/(2.4.6)x^3 + ...

There's also a C program available that uses this method (I've slightly reformatted and corrected it):

/*
** Pascal Sebah : July 1999
**
** Subject:
**
**    A very easy program to compute sqrt(2) with many digits.
**    No optimisations, no tricks, just a basic program to learn how
**    to compute in multiprecision.
**
** Formula:
**
**    sqrt(2) = (239/169)*1/sqrt(1-1/57122)
**
** Data:
**
**    A big real (or multiprecision real) is defined in base B as:
**      X = x(0) + x(1)/B^1 + ... + x(n-1)/B^(n-1)
**      where 0<=x(i)<B
**
** Results: (PentiumII, 450Mhz)
**
**    1000   decimals :   0.02seconds
**    10000  decimals :   1.7s
**    100000 decimals : 176.0s
**
** With a little work it's possible to reduce those computation
** times by a factor of 3 and more.
*/

#include <stdio.h>
#include <stdlib.h>

long B = 10000; /* Working base */
long LB = 4;    /* Log10(base)  */

/*
** Set the big real x to the small integer Integer
*/
void SetToInteger(long n, long* x, long Integer)
{
  long i;
  for (i = 1; i < n; i++)
    x[i] = 0;
  x[0] = Integer;
}

/*
** Is the big real x equal to zero ?
*/
long IsZero(long n, long* x)
{
  long i;
  for (i = 0; i < n; i++)
    if (x[i])
      return 0;
  return 1;
}

/*
** Addition of big reals : x += y
**  Like school addition with carry management
*/
void Add(long n, long* x, long* y)
{
  long carry = 0, i;
  for (i = n - 1; i >= 0; i--)
  {
    x[i] += y[i] + carry;
    if (x[i] < B)
      carry = 0;
    else
    {
      carry = 1;
      x[i] -= B;
    }
  }
}

/*
** Multiplication of the big real x by the integer q
*/
void Mul(long n, long* x, long q)
{
  long carry = 0, xi, i;
  for (i = n - 1; i >= 0; i--)
  {
    xi = x[i] * q;
    xi += carry;
    if (xi >= B)
    {
      carry = xi / B;
      xi -= carry * B;
    }
    else
      carry = 0;
    x[i] = xi;
  }
}

/*
** Division of the big real x by the integer d
**  Like school division with carry management
*/
void Div(long n, long* x, long d)
{
  long carry = 0, xi, q, i;
  for (i = 0; i < n; i++)
  {
    xi    = x[i] + carry * B;
    q     = xi / d;
    carry = xi - q * d;
    x[i]  = q;
  }  
}

/*
** Print the big real x
*/
void Print(long n, long* x)
{
  long i;
  printf("%ld.", x[0]);
  for (i = 1; i < n; i++)
    printf("%04ld", x[i]);
  printf("\n");
}

/*
** Computation of the constant sqrt(2)
*/
int main(void)
{
  long NbDigits = 200000, size = 1 + NbDigits / LB;
  long* r2 = malloc(size * sizeof(long));
  long* uk = malloc(size * sizeof(long));
  long k = 1;
  /*
  ** Formula used:
  **    sqrt(2) = (239/169)*1/sqrt(1-1/57122)
  ** and
  **   1/sqrt(1-x) = 1+(1/2)x+(1.3)/(2.4)x^2+(1.3.5)/(2.4.6)x^3+...
  */
  SetToInteger(size, r2, 1); /* r2 = 1 */
  SetToInteger(size, uk, 1); /* uk = 1 */
  while (!IsZero(size, uk))
  {
    Div(size, uk, 57122); /* uk = u(k-1)/57122 * (2k-1)/(2k) */
    Div(size, uk, 2 * k);
    Mul(size, uk, 2 * k - 1);
    Add(size, r2, uk);    /* r2 = r2+uk */
    k++;
  }
  Mul(size, r2, 239);
  Div(size, r2, 169);  /* r2 = (239/169)*r2 */

  Print(size, r2);     /* Print out of sqrt(2) */

  free(r2);
  free(uk);

  return 0;
}

It takes about a minute to calculate 200,000 digits of sqrt(2).

Note, however, at 200,000 digits the last 11 digits produced are incorrect due to the accumulated rounding errors and you need to run it for 200,012 digits if you want 200,000 correct digits.

Capita

Cool! I guess (1393/985)*1/sqrt(1-1/1940450) would also work. – Lave

Is there other formulas for which convergence is faster ? – Whitesell

(275807/195025)*1/sqrt(1-1/76069501250) works but I don't find faster formula – Whitesell

Answer 3

Here is a solution that computes 1 million digits of sqrt(2) in less than a minute in the good old Prolog programming language. It is based on solving the pell equation, see also here:

p*p+1 = 2*q*q

The recurence relation p′=3p+4q and q′=2p+3q can be cast as a matrix multiplication. Namely we see that if we multiply the vector [p,q] by the matrix of the coefficients we get the vector [p',q']:

| p' |    | 3  4 |   | p |
|    | =  |      | * |   |
| q' |    | 2  3 |   | q |

For matrix A we can use a binary recursion so that we can calculate A^n in O(log n) operations. We will need big nums. I am using this experimental code here whereby the main program is then simply:

/**
  * pell(N, R):
  * Compute the N-the solution to p*p+1=2*q*q via matrices and return
  * p/q in decimal format in R.
  */
pell(N, R) :-
   X is [[3,4],[2,3]],
   Y is X^N,
   Z is Y*[[1],[1]],
   R is Z[1,1]*10^N//Z[2,1].

The following screenshot shows the timing and some results. I was using 10-times a million iterations. One can compare the result with this page here.

Whats missing is still a clear cut criteria and computation that says how many digits are stable. We will need some more time to do that.

Edit 20.12.2016:
We improved the code a little bit by an upper bound of the relative error, and further compute stable digits by nudging the result. The computation time for 1 million digits is now below 2 secs:

?- time(pell(653124, _, S)).
% Uptime 1,646 ms, GC Time 30 ms, Thread Cpu Time 1,640 ms
S = -1000000

Lave

Answer 4

The example you give is accurate in so far as the precision of double arithmetic goes, this is the highest precision most C++ compilers use. In general computers are not tooled to do higher precision calculation. If this is homework of some sort then I suspect you are required to figure out an algorithm for calculating - you need to keep you own array of digits in some way to keep all the precision that you need. If you have some real-world application you should definitely use a high precision library specifically made to do this sort of arithmetic (GMP is a good open-source possibility) - this is a complicated wheel that doesn't need reinventing.

Violetvioleta

Answer 5

This following Javascript function will take an integer and the digits of /precision required after the decimal, then return a string format of the square root.

// Input: a positive integer, the number of precise digits after the decimal point
// Output: a string representing the long float square root

function findSquareRoot(number, numDigits) {
    function get_power(x, y) {
      let result = 1n;
      for (let i = 0; i < y; i ++) {
        result = result * BigInt(x);
      }
      return result;
    }

    let a = 5n * BigInt(number);
    let b = 5n;
    const precision_digits = get_power(10, numDigits + 1);

    while (b < precision_digits) {
      if (a >= b) {
        a = a - b;
        b = b + 10n;
      } else {
        a = a * 100n;
        b = (b / 10n) * 100n + 5n;
      }
    }

    let decimal_pos = Math.floor(Math.log10(number))
    if (decimal_pos == 0) decimal_pos = 1
    let result = (b / 100n).toString()
    result = result.slice(0, decimal_pos) + '.' + result.slice(decimal_pos)    
    return result
  }

Offoffbroadway