From 1c14db584d9a90589ae04889975944e4e7ca3a77 Mon Sep 17 00:00:00 2001 From: John Darrington Date: Wed, 18 Sep 2013 20:37:59 +0200 Subject: [PATCH] Crosstabs: Use logarithms to calculate the product/quotients of factorials Avoid overflow of large numbers resulting from factorials, by using the logarithms as intermediate results --- src/language/stats/crosstabs.q | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/language/stats/crosstabs.q b/src/language/stats/crosstabs.q index 93c28ca7e7..2246d3af17 100644 --- a/src/language/stats/crosstabs.q +++ b/src/language/stats/crosstabs.q @@ -2090,16 +2090,17 @@ display_directional (struct crosstabs_proc *proc, struct pivot_table *pt, /* Statistical calculations. */ -/* Returns the value of the gamma (factorial) function for an integer +/* Returns the value of the logarithm of gamma (factorial) function for an integer argument PT. */ static double -gamma_int (double pt) +log_gamma_int (double pt) { - double r = 1; + double r = 0; int i; for (i = 2; i < pt; i++) - r *= i; + r += log(i); + return r; } @@ -2108,11 +2109,11 @@ gamma_int (double pt) static inline double Pr (int a, int b, int c, int d) { - return (gamma_int (a + b + 1.) / gamma_int (a + 1.) - * gamma_int (c + d + 1.) / gamma_int (b + 1.) - * gamma_int (a + c + 1.) / gamma_int (c + 1.) - * gamma_int (b + d + 1.) / gamma_int (d + 1.) - / gamma_int (a + b + c + d + 1.)); + return exp (log_gamma_int (a + b + 1.) - log_gamma_int (a + 1.) + + log_gamma_int (c + d + 1.) - log_gamma_int (b + 1.) + + log_gamma_int (a + c + 1.) - log_gamma_int (c + 1.) + + log_gamma_int (b + d + 1.) - log_gamma_int (d + 1.) + - log_gamma_int (a + b + c + d + 1.)); } /* Swap the contents of A and B. */ -- 2.30.2