title: Inconsistent value returned from C++ function vs Python function for skew normal distribution
body: <p>I have a function to estimate the alpha parameter from the skew normal distribution in C++ and Python. The Python function is written using NumPy and the C++ function uses the STL. My issue is that my C++ implementation is giving me incorrect results. The two functions are essentially identical but the Python version gives me correct results whereas the C++ does not - I have investigated this in some detail and I cannot come to a conclusion as to what's causing the error, any help would be great. </p>
  
  <p><strong>Python Function</strong></p>
  
  <pre><code>import numpy as np
  def convert_to_alpha(skew):
      a = np.pi/2 
      skew_ = abs(skew)
      numerator = np.power(skew_, (2/3)) 
      b = (4-np.pi)/2
      b = np.power(b, (2/3))
      denom = numerator + b
      delta = np.sqrt(a * (numerator/denom))
      a = delta/np.sqrt((1-np.power(delta, 2)))
      return a * np.sign(skew)
  </code></pre>
  
  <p><strong>C++ Function</strong></p>
  
  <pre><code>double convert_to_alpha(double skew)
  {
      double pi = 3.141592653589793;
      double a = pi / 2;
      double skew_ = std::abs(skew);
      double numerator = std::pow(skew_, (2 / 3));
      double b = (4 - pi) / 2;
      b = std::pow(b, (2 / 3));
      double denom = numerator + b;
      double delta = std::sqrt(a * (numerator / denom));
      double alpha = delta / std::sqrt((1 - std::pow(delta, 2)));
      if (skew == 0) { return 0; }
      else if (std::signbit(skew) == 1) { return -1 * alpha; } 
      else return alpha; 
  }
  </code></pre>
  
  <p>The Python function returns the values I would expect whereas the C++ function does not, as examples for input 0.99 I'd expect 27.85xxxx or for input 0.5 I'd expect 2.17xxxx which is exactly what I get from the Python implementation, C++ gives me 1.91306. </p>
  
  <p>Also, strangely - regardless of the input, the C++ implementation seems to return 1.91306. </p>
  
  <p><strong>Driver code for C++</strong></p>
  
  <pre><code>#include &lt;cmath&gt;
  #include &lt;math.h&gt;
  #include &lt;iostream&gt;
  
  int main()
  {
      double convert_to_alpha(double skew);
      std::cout &lt;&lt; "skew: " &lt;&lt; convert_to_alpha(0.99);
      return 0;
  }
  double convert_to_alpha(double skew)
  {
      double pi = 3.141592653589793;
      double a = pi / 2;
      double skew_ = std::abs(skew);
      double numerator = std::pow(skew_, (2 / 3));
      double b = (4 - pi) / 2;
      b = std::pow(b, (2 / 3));
      double denom = numerator + b;
      double delta = std::sqrt(a * (numerator / denom));
      double alpha = delta / std::sqrt((1 - std::pow(delta, 2)));
      if (skew == 0) { return 0; } // if skew is 0 return 0 
      else if (std::signbit(skew) == 1) { return -1 * alpha; } // if skew is negative return -alpha
      else return alpha; // if skew is positive return alpha
  }
  </code></pre>
  
  <p>I'd expect the results to be very similar, definitely not as different as they are currently. I have not encountered an issue like this before so any help figuring out what's causing the inconsistency with the C++ implementation would be very helpful.</p>