digamma: boost vs mlpack

digamma_comp.md

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