Unable to do tail recursion in C++ using clang with -O3 option

Question

I am unable to tell the clang compiler to perform tail-recursive optimisation in my C++ code. I saw this post Which, if any, C++ compilers do tail-recursion optimization?, where the advice is to use the -O3 flag with clang and it would do it. But in spite of that, I am having stack overflow for large inputs to the below primality testing code.

#include <iostream>
#include <chrono>
#include <gmp.h>
#include <string>

bool mpz_is_divisible(mpz_t n, mpz_t d);

void mpz_sqr(mpz_t res, mpz_t x) {
   //computes square of x and puts the result into res

   mpz_pow_ui(res, x, 2);
}

//the below function causes stack overflow for large inputs

void smallest_divisor(mpz_t n, mpz_t div, mpz_t div_sqr, mpz_t smallest_div) {
 //finds the smallest number which divides n and puts the result into smallest_div

  mpz_sqr(div_sqr, div);
  
  if (mpz_cmp(div_sqr, n) > 0) {
     mpz_set(smallest_div, n);
     return;
  }

  if (mpz_is_divisible(n, div)) {
    mpz_set(smallest_div, div);
    return;
  }

  mpz_add_ui(div, div, 1);
  smallest_divisor(n, div, div_sqr, smallest_div); //<-- should do tail recursion optimisation?
}

bool mpz_prime_test_basic(mpz_t n) {
  //checks if n is prime

  mpz_t div, div_sqr, smallest_div;

  mpz_inits(div, div_sqr, smallest_div, NULL);
  mpz_set_ui(div, 2);

  smallest_divisor(n, div, div_sqr, smallest_div);

  if (mpz_cmp(smallest_div, n) == 0) {
    mpz_clear(div);
    mpz_clear(div_sqr);
    mpz_clear(smallest_div);
    return true;
  }
  mpz_clear(div);
  mpz_clear(div_sqr);
  mpz_clear(smallest_div);
  return false;
}

bool mpz_is_divisible(mpz_t n, mpz_t d) {
  //checks if n is divisible by d
  mpz_t rem;
  mpz_init(rem);
  mpz_tdiv_r(rem, n, d);
  int cmp = mpz_cmp_si(rem, 0); //checks if remainder is equal to 0
  if (cmp == 0) return true;
  return false;
}

int main() {
  std::string num_str;
  mpz_t num;
  mpz_init(num);
  std::cout << "Enter number to check" << '\n';
  std::cin >> num_str;
  mpz_set_str(num, num_str.c_str(), 10);

  bool is_prime = mpz_prime_test_basic(num);

  if (is_prime) {
    std::cout << num_str << " is a prime\n";
  } else {
    std::cout << num_str << " is not a prime\n";
  }

  return 0;
}

I am going through chapter 1 of SICP right now and there's an exercise (1.22) about checking primality using an O(√n) algorithm. Tail-recursion is enforced by the Scheme standard, hence there is no problem there for large inputs but I am doing the same problem in C++, and for large numbers, the smallest-divisor function consumes the stack. I am using the GMP library for big integer arithmetic.

Is it possible to enforce tail recursion in such a program?

Thanks

Answer 1

Something is odd about mpz_is_divisible. For one, you forget to free the memory held by rem. Even if you add the mpz_clear(rem) call, you don't get TCO in smallest_divisor. What you have to do is mark it __attribute__((noinline)) to get Clang to do the optimization (though GCC does it without it). Very strange, and also moot, since mpz_is_divisible is just a worse version of mpz_divisible_p.

void smallest_divisor(mpz_t n, mpz_t div, mpz_t div_sqr, mpz_t smallest_div) {    
    mpz_sqr(div_sqr, div);
    if (mpz_cmp(div_sqr, n) > 0) {
        mpz_set(smallest_div, n);
        return;
    }
    if (mpz_divisible_p(n, div)) {
        mpz_set(smallest_div, div);
        return;
    }
    mpz_add_ui(div, div, 1);
    smallest_divisor(n, div, div_sqr, smallest_div); // gets TCO
}

Godbolt (steals GMP from trunk GCC, may have to update some paths with the current date to make it work)

I've also cleaned up the code a lot to produce a new version. In particular, my FP intuition tells me that smallest_divisor is really not supposed to be its own function. It's a recursive worker that belongs inside of mpz_prime_test_basic. We can use this answer to write local recursive definitions. I couldn't find gmpxx.h on Godbolt, so I also wrote a clone of mpz_class. Then we can write

bool is_prime(mpz_t n) noexcept {
    mpz div(2L);
    fix{[](auto rec, mpz_t n, mpz_t div, mpz_t div_sqr) noexcept -> void {
        mpz_mul(div_sqr, div, div);
        if(mpz_cmp(div_sqr, n) > 0) mpz_set(div, n);
        else if(!mpz_divisible_p(n, div)) {
            mpz_add_ui(div, div, 1);
            rec(n, div, div_sqr);
        }
    }}(n, div, mpz());
    return mpz_cmp(div, n) == 0;
}

Since we're not actually changing the arguments between the calls, just mutating behind them, we can also just write

bool is_prime(mpz_t n) noexcept {
    mpz div(2L), div_sqr;
    fix{[&](auto rec) noexcept -> void {
        mpz_mul(div_sqr, div, div);
        if(mpz_cmp(div_sqr, n) > 0) mpz_set(div, n);
        else if(!mpz_divisible_p(n, div)) {
            mpz_add_ui(div, div, 1);
            rec();
        }
    }}();
    return mpz_cmp(div, n) == 0;
}

Godbolt (same caveat)

Answer 2

In principle, all major C++ compilers make tail call optimizations. And your code should benefit from it.

Having looked at the assembly obtained from your code on Godbolt with Clang 11, it appears that the -O3 flag indeed leads to tail code optimization. Your recursive function is compiled as:

smallest_divisor(test, test, test, test): # @smallest_divisor(test, test, test, test)
  push rax
.LBB1_1: # =>This Inner Loop Header: Depth=1  <------ TCO Loop 
  mov edi, 2
  call mpz_pow_ui(test, test, unsigned int)
  call mpz_cmp(test, test)
  test eax, eax
  jg .LBB1_4                                  <------ go to return sequence
  call mpz_init(test)
  call mpz_tdiv_r(test, test, test)
  xor edi, edi
  call mpz_cmp_si(test, int)
  test eax, eax
  je .LBB1_4
  mov edi, 1
  call mpz_add_ui(test, test, unsigned int)
  jmp .LBB1_1                                <--- no recursive call but loop
.LBB1_4:                                     <--- this is the end of recursion
  pop rax
  jmp mpz_set(test, test) # TAILCALL         <-- no ret, but another tail call :-)

I cannot test your code, since I do not have the library and do not have the time to invest in getting it and installing it.

From the symptoms you describe, I can see the following possible root causes:

• there's some UB in the library, leading to a segfault
• optimization is disable due to compiler settings. If you're using XCode for instance, you'll need to go to menu Product > Scheme > Edit Scheme... and then chose "Release" as build configuration.

Having tried the assembly in XCode (Product > Perform action > Assemble), the generated assembly code shows evidence of TCO:

__Z16smallest_divisor4testS_S_S_:       ## @_Z16smallest_divisor4testS_S_S_
Lfunc_begin1:
    .loc    67 29 0                 ## Test0105-1/main.cpp:29:0
    .cfi_startproc
## %bb.0:
    pushq   %rbp
    .cfi_def_cfa_offset 16
    .cfi_offset %rbp, -16
    movq    %rsp, %rbp
    .cfi_def_cfa_register %rbp
LBB1_1:                                 ## =>This Inner Loop Header: Depth=1
Ltmp3:
    .loc    67 24 4 prologue_end    ## Test0105-1/main.cpp:24:4
    movl    $2, %edi
    callq   __Z10mpz_pow_ui4testS_j
   ....
    movl    $1, %edi
    callq   __Z10mpz_add_ui4testS_j
Ltmp13:
    .loc    67 0 3 is_stmt 0        ## Test0105-1/main.cpp:0:3
    jmp LBB1_1                             <<<<<<<<<<<< Tail call optimization (loop)
LBB1_4:
    popq    %rbp
Ltmp14:
    jmp __Z7mpz_set4testS_      ## TAILCALL