Convert a recursive function to a tail recursive function

Question

In continuum to this question:

I have a function in c++ that calls itself over and over again. This is it:

#include <iostream>
#include <math.h>
#include <stdio.h>
#include <cmath>
using namespace std;
double g( double a, double x){
    if (x>=a) return (g(a,x-1)+g(a,x-a));
    else if (x<a) return 1;
    return 0; //Never Reached
}
int main(){
    cout << (unsigned long)g(sqrt(90),90) <<endl; // outputs 7564511
    cout << (unsigned long)g(sqrt(10000019),10000019)<<endl; // Signal: SIGSEGV (Segmentation fault)
}

I would like to know how this function can be converted into some sort of iteration or tail loop (or anything that stops the segfault), but more importantly I need to know how to actually do it yourself.

---

NOTE: I apologize in advance if this is a trivial question.

NOTE2: There are similar questions (like this or this), but none of the ones I found address the fact that my function calls itself twice each iteration.

Answer 1

As many person have already stated, this cannot be directly converted to a tail recursive or an iterative function as floating point function arguments make it difficult to build the result iteratively. However, with a bit of thought the function can be calculated quite efficiently.

First, because all the recursions are summed and end with 1, the function basically calculates the number of paths to the recursion end. For example, for g(5,2), one path would be g(2,5)-> g(2,3) -> g(2,1) (returns 1 here) and an other path g(5,2)-> g(4,2) -> g(3,2) -> g(2,2) -> g(0,2). So to calculate g we need to just calculate the number of the possible paths.

Let's start with the path where we subtract always a from x. Clearly we have only exactly one this kind of path. Next, consider the case where we once we choose the path where we subtract 1 and other times we subtract a, we have floor((x-a)/a) positions to select the 1. Hence, there is floor((x-a)/a) possible paths in this case. In a next iteration, we want to select step 1 twice. There are n*(n-1)/2 combination, where n=floor((x-1-a)/a) and n*(n-1)/2 is the binomial coefficient \binom{n,2} . Next step with three ones there are \binom{n,3} combinations where n is now=floor((x-2-a)/a) etc.

If you pre-calculate the binomial coefficient the algorithm is O(x), hence it probably could also calculate g(sqrt(10000019),10000019). However, the largest native c++ integer type (unsigned long long) overflows already around g(sqrt(500),500). You can use long double to get an approximate value for slightly larger input. Alternatively you can use boost Multiprecision Library to get more digits, but I assume you will run out of memory before you get the answer to g(sqrt(10000019),10000019).

The source code with an overflow check to calculate g()

#include <iostream>
#include <vector>
#include <limits>
#include <algorithm>
#include <cstdlib>

unsigned long long binomial(unsigned int n, unsigned int m) {
  if (n - m < m) m = n - m;
  std::vector<unsigned long long>bc(m+1, 0);
  bc[0] = 1;
  for (unsigned int i = 0;i <= n;++i) {
    for (unsigned int j = std::min(i, m);j > 0;j--) {
      if (std::numeric_limits<unsigned long long>::max()-bc[j-1] < bc[j]) {
        std::cout << "Error: too large binomial coefficient(" << n << "," << m << ")" << std::endl;
        exit(1);
      }
      bc[j] += bc[j - 1];
    }   
  }
  return bc[m];
}

unsigned long long g(double a, double x) {
  unsigned long long r = 1;
  int n = 0;
  for (int i = static_cast<int>(x);i >= a;--i) {
    ++n;
    int m = static_cast<int>((i - a) / a);
    unsigned long long b = binomial(m + n, n);
    if (std::numeric_limits<unsigned long long>::max() - b < r) {
      std::cout << "Error: too large sum " << b << "+" << r << std::endl;
      exit(1);
    }
    r += b;
  }
  return r;
}

int main(){
  std::cout << g(sqrt(90), 90) << std::endl;
  std::cout << g(sqrt(10000019), 10000019) << std::endl;
}

Answer 2

I guess just for the reference, here's implementation w/o recursion. I am not sure it will ever finish with the given inputs though:

#include <stack>
double gnr(double a, double x) {
   std::stack<double> stack;
   double result = 0;
   stack.push(x);
   while (!stack.empty()) {
      x = stack.top();
      stack.pop();
      //cout << stack.size() << " " << x << endl;
      if (x < a) {
         result++;
      } else {
         stack.push(x - 1);
         stack.push(x - a);
      }
   }
   return result;
}

Answer 3

Memorization can be an effective way to limit the number of recursive calls required for your computation. Try evaluating the function for some simple inputs, like g(2, 8) and you'll see that you end up evaluating the function for the same values over and over. By caching the result for each set of inputs the first time you calculate it, you can short circuit the recursion and dramatically reduce the size of the problem.

To make the function iterative rather than recursive, one strategy you can use is to try to turn the definition around and iterate from the bottom up. Consider the Fibonacci function:

fib(n) = fib(n-1) + fib(n-2)

To calculate fib(n) iteratively, we start from the base cases fib(1) + fib(0) and iterate up to fib(n). That lets you accumulate the value as you go, rather than having to remember where you were (over and over and over again) while you calculate intermediate values. So an iterative definition of fib() looks like:

fib(n) {
    a = 1;
    b = 0;
    fib = 0;
    i = 1;
    while (i < n) {
        fib = a + b;
        b = a;
        a = fib;
        i++;
    }
    return fib;
}

You should be able to do something similar with your g() function. I don't have time to play with it, but I bet that if you try evaluating it by hand for a few a, x pairs you'll notice a pattern that'll let you rewrite the function in an iterative form the way I did above for fib().