Run time errors while computing e^x

Question

Continuing from the question I posted before How to increase precision for function e^x

I made few changes in my code by taking advices given

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

long double exponential(long double x, long double n, long double p)
{
    long double i = 1;

    while (n > 0) // loop stops when becomes less than 0
        i = i + (x / p) * (exponential(x, n - 0.0001, p + 1));

    if (n < 0) // when n reaches 0 or becomes less than zero value of i will be returned
        return i;
}

int main()
{
    long double p, x, n;
    scanf("%Lf", &x);
    printf("math.h e^x =     %lf\n", exp(x));
    printf("calculated e^x = %Lf\n", exponential(x, 1, 1));
    return 0;
}

But I am not getting any output its just giving run time error(http://codepad.org/jIKoYGFC) and I don't know why . Please can some one help me why I am getting these errors

The function does not return anything if `n == 0`. Try removing `if (n < 0)`. Moreover, `n` is not modified in the `while(.)`loop. How can it stop? — Damien, Feb 23 '22 at 09:02

score 3 · Accepted Answer · answered Feb 23 '22 at 09:24

3

That loop is completely bogus. You're not writing an iterative function (where it would make more sense). Further you have an edge case of zero returning undefined content.

Though I do not recommend floating point for loop control, nor do I advise thousands of invocations into recursive calls, your code should be more like

long double exponential(long double x, long double n, long double p)
{
    long double i = 1;
    if (n > 0)
        i += (x / p) * (exponential(x, n - 0.0001, p + 1));
    return i;
}

which ultimately is just this:

long double exponential(long double x, long double n, long double p)
{
    return 1 + ((n > 0) ? (x / p) * exponential(x, n - 0.0001, p + 1) : 0);
}

Fixing that (either way):

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

long double exponential(long double x, long double n, long double p)
{
    return 1 + ((n > 0) ? (x / p) * exponential(x, n - 0.0001, p + 1) : 0);
}

int main()
{
    long double x = 5;
    printf("math.h e^x =     %Lf\n", expl(x));
    printf("calculated e^x = %Lf\n", exponential(x, 1, 1));
    return 0;
}

Output

math.h e^x =     148.413159
calculated e^x = 148.413159

answered Feb 23 '22 at 09:24

WhozCraig

65,258
11
75
141

Thanks for correcting my code , I assume you have seen my last question , can you just how can I increase the precision of the function since its deviating after x= 30 – user Feb 23 '22 at 09:53
Pretty sure you're capping the precision at e^30. Most x86 and familiar processes use 80-bit extended precision, which has a 1-bit integer part, and 63bit fractional part (and a 15 bit exponent, not that you'd care). E.g. you're going to get about 19 significant digits out of the deal. If you look at where things start falling off the tails, it's actually at x=29, i.e. 3931334297144.042074. At x=30 you're at 10686474581524.462147, and now bleeding fractional bits. – WhozCraig Feb 23 '22 at 10:15
@user You should be aware that getting code like this to work really well is **hard**. Floating-point work can be tricky enough to begin with — and summing an infinite series to compute something like exp(x) is trickier still. If you're not careful, the various x^N or N! terms blow up. If you're not careful, you lose more and more of your precision to noise. Successfully avoiding all the pitfalls, and achieving the holy grail of 1/2 ulp precision, can require a *lot* of careful study and work. I'm not suggesting that it's impossible or that you shouldn't try, but it's not simple! – Steve Summit Feb 23 '22 at 12:52

Gerhard · Answer 2 · 2022-03-17T08:50:21.800

Created a version using LibGMP out of curiosity to see how this will improve the outcome.

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

#define BITS 128

The line above mean mantissa is 128 bits and exponent of 64 bits. This is my GMP based function:

void gmp_exponential(mpf_t * i, double x, double n, int p)
{
    mpf_t ii;
    mpf_init(ii);    
    mpf_set_d(ii,1);
    
    if (n > 0){
      mpf_t a,b,c;

      gmp_exponential(&ii, x, n - 0.0001, p + 1);
      
      mpf_inits (a, b, c, NULL);  
      mpf_set_d(a,x);
      mpf_div_ui(b, a, p) ;
      mpf_mul(c, b, ii);
      mpf_add (*i,*i,c);
      
      mpf_clears(a,b,c,NULL);
    }
    mpf_clear(ii);
}

Reusing the function from WhozCraig for comparison:

long double exponential(long double x, long double n, long double p)
{
    return 1 + ((n > 0) ? (x / p) * exponential(x, n - 0.0001, p + 1) : 0);
}

And code to run it all:

int main()
{
  double x = 30.0;
  
  mpf_t i;
  mpf_init2 (i, BITS);    
  mpf_set_d(i,1);
  gmp_exponential(&i, x, 1, 1);
  
  printf     ("math.h e^x              =  %Lf\n", expl(x));
  gmp_printf ("calculated e^x with mpf =  %.*Ff\n", 6, i);
  printf     ("calculated e^x          =  %Lf\n", exponential(x, 1, 1));
  
  mpf_clear(i);
  return 0;
}

Output:

math.h e^x              =  10686474581524.462147
calculated e^x with mpf =  10686474581524.462143
calculated e^x          =  10686474581524.462149

... and the best approximation goes to ... standard library :) — pmg, Mar 17 '22 at 09:08
@pmg: Also checked this with `mpfr_exp` in lib MPRF. Roll your own results are just less accurate. — Gerhard, Mar 17 '22 at 14:19

Run time errors while computing e^x

2 Answers2