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?

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

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.

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

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.