Skip to content

Instantly share code, notes, and snippets.

@heisenbuug
Last active September 1, 2021 22:32
Show Gist options
  • Save heisenbuug/89ad02301d2c38f7723caf0e7232d6ab to your computer and use it in GitHub Desktop.
Save heisenbuug/89ad02301d2c38f7723caf0e7232d6ab to your computer and use it in GitHub Desktop.
digamma: boost vs mlpack

Implementation of digamma

template<typename T>
T digamma(T x)
{
  T result = 0;

  // Check for negative arguments and use reflection
  if (x <= -1)
  {
    // Reflect
    x = 1 - x;
    // Argument reduction for tan
    T reminder = x - floor(x);
    // Shift to negative if x > 0.5
    if (reminder > 0.5)
      reminder -= 1;

    // Check for evaluation at negative pole
    if (reminder == 0)
    {
      // throw proper exception
      throw std::runtime_error("Evaluation of fucntion at pole");
    }

    result = M_PI / tan(M_PI * reminder);
  }
  
  if (x == 0)
  { 
    // throw exception
    throw std::runtime_error("Evaluation of fucntion at pole");
  } 
  
  // Direct formula to calculate digamma
  result += std::log(x) - (1 / (2 * x));

  return result;
}

Output for x = 1 to x = 20

x mlpack boost
1 -0.5 -0.577216
2 0.443147 0.422784
3 0.931946 0.922784
4 1.26129 1.25612
5 1.50944 1.50612
6 1.70843 1.70612
7 1.87448 1.87278
8 2.01694 2.01564
9 2.14167 2.14064
10 2.25259 2.25175
11 2.35244 2.35175
12 2.44324 2.44266
13 2.52649 2.526
14 2.60334 2.60292
15 2.67472 2.67435
16 2.74134 2.74101
17 2.8038 2.80351
18 2.86259 2.86234
19 2.91812 2.91789
20 2.97073 2.97052

Ouput for x = 0.1 to 0.9

x mlpack boost
0.1 -10.3592 -10.4238
0.2 -5.23434 -5.28904
0.3 -3.45558 -3.50252
0.4 -2.52067 -2.56138
0.5 -1.92787 -1.96351
0.6 -1.50916 -1.54062
0.7 -1.19206 -1.22002
0.8 -0.939991 -0.965008
0.9 -0.732415 -0.754927

Output for x = -1.1 to -1.9

x mlpack boost
-1.1 10.1727 10.1542
-1.2 4.88521 4.86832
-1.3 2.89802 2.88254
-1.4 1.6879 1.67367
-1.5 0.716291 0.703156
-1.6 -0.257563 -0.269719
-1.7 -1.47444 -1.48572
-1.8 -3.47299 -3.48349
-1.9 -8.77654 -8.78635

Output for x = -0.1 to -0.9

x mlpack boost
-0.1 9.33908 9.24507
-0.2 4.15186 4.03499
-0.3 2.26237 2.11331
-0.4 1.15584 0.959381
-0.5 0.306853 0.03649
-0.6 -0.499624 -0.894718
-0.7 -1.44207 -2.07395
-0.8 -2.85944 -4.03904
-0.9 -6.19148 -9.31265

Output for x = -20.1 to -20.9

x mlpack boost
-20.1 12.6944 12.6942
-20.2 7.35443 7.35424
-20.3 5.31772 5.31753
-20.4 4.06078 4.06059
-20.5 3.04478 3.0446
-20.6 2.02875 2.02858
-20.7 0.77173 0.771553
-20.8 -1.26514 -1.26532
-20.9 -6.60553 -6.6057
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment