0

I have the following code and it returns: 2.7182818284590455348848081484902650118 where only 2.718281828459045 is right and the rest is wrong.

I'm using the following reference: 2.7182818284590452353602874713526624977 http://www.vaxasoftware.com/doc_edu/mat/nume15000.pdf

My code is:

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

#define EPSILON 1.0e-15

int main() {
    unsigned long long fact = 1;
    double e = 2.0, e0;
    int n = 2;
    do {
        e0 = e;
        fact *= n++;
        e += 1.0 / fact;
    }
    while (fabs(e - e0) >= EPSILON);
    printf("e = %.37f\n", e);
    return 0;
}

I appreciate if anyone can help me. Thank you so much!

emacs drives me nuts
  • 2,785
  • 13
  • 23
Brayan Angarita
  • 307
  • 2
  • 12
  • 3
    You will need to write extended-precision software or get and use such software, such as [GMP](https://gmplib.org). The `double` format in most C implementations has a precision around one part in 2^52, which is about 10^15.65. For 37 digits, you need precision of one part in 10^37 just for the result, and more for intermediate calculations. Also, you will need to change `EPSILON`, as of course terminating the calculation when the change in the sum reaches 10^−15 will be insufficient to calculate e to one part in 10^37. – Eric Postpischil Feb 12 '20 at 01:44
  • Thank you for your comment ... Yes sir, I'm using GMP, but I really don't know how I can extend the functionality – Brayan Angarita Feb 12 '20 at 01:45
  • 6
    You say you are using GMP, but the code you show in the question is not using GMP, except for including the header. Include the header does not change how C’s built-in types like `double` work. You need to call routines from the GMP library. To do that, start by reading the GMP documentation. – Eric Postpischil Feb 12 '20 at 01:46
  • Oh yes, you're right. I will read the documentation. I'm starting with GMP and C. Thanks for your help. – Brayan Angarita Feb 12 '20 at 01:52
  • 1
    I suspect that a lot of the problem is that it's doing additions (at `e += 1.0 / fact;`) in "largest to smallest" order. To improve precision you'd want the opposite, which isn't pretty (e.g. maybe store the values of `1/fact` in an array, then have a second loop that does "e = sum of values in array"?). – Brendan Feb 12 '20 at 04:11

3 Answers3

1

How can calculate 37 decimals for Euler in C?

2128 is about 3.4*1038. By forming a 128 bit integer numerator/denominator we can achieve the desired e.

Re-writing the loop as below and running up to a large enough n until num*i + 1 is about to overflow (n==33), we can arrive at/near the desired result.

// Algorithm
wide_floating_point e() {
  unsigned n = 33;
  my_uint128 num = 1;
  my_uint128 den = 1;
  for (unsigned i=1; i <= n; i++) {
    num = num*i + 1;
    den = den*i;
  }
  return (wide_floating_point) num/fact;
}

Now OP, may not have a 128 unsigned integer type, nor able to perform a 128/128 bit floating point divide. Yet OP only needs to create 2 functions (implementations not shown - but are essentially grade school math):

typedef struct {
  // or unsigned char decimal_digit[40] or 
  // however you want to save a big integer
  uint64_t hi,lo;
} my_uint128; 

my_uint128 times_add_128(my_uint128 x, unsigned m, bool a);
void print_quotient_128(my_uint128 num, my_uint128 dem);

If we test the algoithm and use __int128 and use a long double (with 80 precision) we get close to the goal.

long double e(int n) {
  unsigned __int128 fact = 1;
  unsigned __int128 num = 1;
  for (int i=1; i<=n; i++) {
    fact *= i;
    if (num >  (((unsigned __int128)-1)-1)/i) { printf("%d!\n", i); exit(0); }
    num = num*i + 1;
  }
  return 1.0L*num/fact;
}

int main() {
  for (int i=1; i<34; i++)  {
    printf("%d %.25Lf\n", i, e(i));
  }
}

   2,718281828459045235 3602874713526624977
33 2.718281828459045235 4281681

I leave the 2 functions times_add_128(), print_quotient_128() for OP to code:

chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256
1

My code is:

#include "gmp.h"

Just including this header does not use multi-precision arithmetic by magic; you'll have to use appropriate functions as documented in the GNU MP manual!

For example, you would use

{
  mpf_t x;
  mpf_init2 (x, 100);
  ...

To initialize floating-point x for precision of 100 bits, and for computation you'd use functions from here like mpf_add to add 2 such variables.

emacs drives me nuts
  • 2,785
  • 13
  • 23
-1

I decrease your value in.

printf("e = %.37f\n", e);

to

printf("e = %.15f\n", e);

answer is:

2,718281828459046

Liaquat Hussain
  • 1
  • 3