add a function to compute log(n) for integer n

This commit adds the function ln() to compute log(n) for integer n. It uses a
lookup table for n < 100 to speed things up.

4 files changed, 146 insertions, 2 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 6ddd576..76364fd 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -83,7 +83,7 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
      * different transit times across the face of the PMT, it seems better to
      * convolve first which is what we do here. In addition, the problem is not
      * analytically tractable if you do things the other way around. */
-    return log(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma)) + log(f(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma));
+    return ln(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma)) + log(f(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma));
 }
 
 static double gsl_muon_time(double x, void *params)
@@ -414,7 +414,7 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
                 logp[j] = log(pq(ev->pmt_hits[i].qhs,j)) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect, &mu_direct[i], 1, &ts[i], tmean, 1.5);
                 if (j == 1 || logp[j] > max_logp) max_logp = logp[j];
 
-                if (logp[j] - max_logp < MIN_RATIO*log(10)) {
+                if (logp[j] - max_logp < MIN_RATIO*ln(10)) {
                     j++;
                     break;
                 }
diff --git a/src/misc.c b/src/misc.c
index 82d114b..823bc78 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -6,6 +6,113 @@
 static struct {
     int n;
     double f;
+} ln_table[LN_MAX + 1] = {
+    {0,-INFINITY},
+    {1,0},
+    {2,0.69314718055994529},
+    {3,1.0986122886681098},
+    {4,1.3862943611198906},
+    {5,1.6094379124341003},
+    {6,1.791759469228055},
+    {7,1.9459101490553132},
+    {8,2.0794415416798357},
+    {9,2.1972245773362196},
+    {10,2.3025850929940459},
+    {11,2.3978952727983707},
+    {12,2.4849066497880004},
+    {13,2.5649493574615367},
+    {14,2.6390573296152584},
+    {15,2.7080502011022101},
+    {16,2.7725887222397811},
+    {17,2.8332133440562162},
+    {18,2.8903717578961645},
+    {19,2.9444389791664403},
+    {20,2.9957322735539909},
+    {21,3.044522437723423},
+    {22,3.0910424533583161},
+    {23,3.1354942159291497},
+    {24,3.1780538303479458},
+    {25,3.2188758248682006},
+    {26,3.2580965380214821},
+    {27,3.2958368660043291},
+    {28,3.3322045101752038},
+    {29,3.3672958299864741},
+    {30,3.4011973816621555},
+    {31,3.4339872044851463},
+    {32,3.4657359027997265},
+    {33,3.4965075614664802},
+    {34,3.5263605246161616},
+    {35,3.5553480614894135},
+    {36,3.5835189384561099},
+    {37,3.6109179126442243},
+    {38,3.6375861597263857},
+    {39,3.6635616461296463},
+    {40,3.6888794541139363},
+    {41,3.713572066704308},
+    {42,3.7376696182833684},
+    {43,3.7612001156935624},
+    {44,3.784189633918261},
+    {45,3.8066624897703196},
+    {46,3.8286413964890951},
+    {47,3.8501476017100584},
+    {48,3.8712010109078911},
+    {49,3.8918202981106265},
+    {50,3.912023005428146},
+    {51,3.9318256327243257},
+    {52,3.9512437185814275},
+    {53,3.970291913552122},
+    {54,3.9889840465642745},
+    {55,4.0073331852324712},
+    {56,4.0253516907351496},
+    {57,4.0430512678345503},
+    {58,4.0604430105464191},
+    {59,4.0775374439057197},
+    {60,4.0943445622221004},
+    {61,4.1108738641733114},
+    {62,4.1271343850450917},
+    {63,4.1431347263915326},
+    {64,4.1588830833596715},
+    {65,4.1743872698956368},
+    {66,4.1896547420264252},
+    {67,4.2046926193909657},
+    {68,4.219507705176107},
+    {69,4.2341065045972597},
+    {70,4.2484952420493594},
+    {71,4.2626798770413155},
+    {72,4.2766661190160553},
+    {73,4.290459441148391},
+    {74,4.3040650932041702},
+    {75,4.3174881135363101},
+    {76,4.3307333402863311},
+    {77,4.3438054218536841},
+    {78,4.3567088266895917},
+    {79,4.3694478524670215},
+    {80,4.3820266346738812},
+    {81,4.3944491546724391},
+    {82,4.4067192472642533},
+    {83,4.4188406077965983},
+    {84,4.4308167988433134},
+    {85,4.4426512564903167},
+    {86,4.4543472962535073},
+    {87,4.4659081186545837},
+    {88,4.4773368144782069},
+    {89,4.4886363697321396},
+    {90,4.499809670330265},
+    {91,4.5108595065168497},
+    {92,4.5217885770490405},
+    {93,4.5325994931532563},
+    {94,4.5432947822700038},
+    {95,4.5538768916005408},
+    {96,4.5643481914678361},
+    {97,4.5747109785033828},
+    {98,4.5849674786705723},
+    {99,4.5951198501345898},
+    {100,4.6051701859880918},
+};
+
+static struct {
+    int n;
+    double f;
 } ln_fact_table[LNFACT_MAX + 1] = {
     {0,0},
     {1,0},
@@ -110,6 +217,18 @@ static struct {
     {100,363.73937555556347},
 };
 
+double ln(unsigned int n)
+{
+    /* Returns the logarithm of n.
+     *
+     * Uses a lookup table to return results for n < 100. */
+
+    if (n <= LN_MAX)
+        return ln_table[n].f;
+
+    return log(n);
+}
+
 double lnfact(unsigned int n)
 {
     /* Returns the logarithm of n!.
diff --git a/src/misc.h b/src/misc.h
index eb62e46..6e9095a 100644
--- a/src/misc.h
+++ b/src/misc.h
@@ -3,8 +3,10 @@
 
 #include <stdlib.h> /* for size_t */
 
+#define LN_MAX 100
 #define LNFACT_MAX 100
 
+double ln(unsigned int n);
 double lnfact(unsigned int n);
 double kahan_sum(double *x, size_t n);
 double interp1d(double x, double *xp, double *yp, size_t n);
diff --git a/src/test.c b/src/test.c
index fa4256f..9976f1f 100644
--- a/src/test.c
+++ b/src/test.c
@@ -489,6 +489,28 @@ err:
     return 1;
 }
 
+int test_ln(char *err)
+{
+    /* Tests the lnfact() function. */
+    size_t i;
+    double x, expected;
+
+    for (i = 1; i < 1000; i++) {
+        x = ln(i);
+        expected = log(i);
+
+        if (!isclose(x, expected, 1e-9, 1e-9)) {
+            sprintf(err, "ln(%zu) returned %.5g, but expected %.5g", i, x, expected);
+            goto err;
+        }
+    }
+
+    return 0;
+
+err:
+    return 1;
+}
+
 struct tests {
     testFunction *test;
     char *name;
@@ -505,6 +527,7 @@ struct tests {
     {test_interp1d, "test_interp1d"},
     {test_kahan_sum, "test_kahan_sum"},
     {test_lnfact, "test_lnfact"},
+    {test_ln, "test_ln"},
 };
 
 int main(int argc, char **argv)