Advanced Search

Difference of two numbers in log form

/ Published in: C
Save to your folder(s)
Given log(x) and log(y) compute log(x - y)

Expand | Embed | Plain Text

Copy this code and paste it in your HTML 
  1. /* Copyright (C) 2009 by Andrew Miller ([email protected]).
  2.  * You may use this code under the GNU GPL version 3 (or at your option, any later version)
  3.  *   See http://www.gnu.org/licenses/gpl-3.0.txt
  4.  * You may alternatively elect to use this code under the GNU LGPL version 3  (or at your 
  5.  *   option, any later version) at http://www.gnu.org/licenses/lgpl-3.0.txt .
  6.  * You may alternatively elect to use this code under the Mozilla MPL version 1.1 or later
  7.  *  (http://www.mozilla.org/MPL/MPL-1.1.html)
  8.  *
  9.  */
  10. /* Given log(x) and log(y) compute log(x - y). */
  11. double
  12. log_form_subtract(double log_x, double log_y)
  13. {
  14.   if (log_x <= log_y)
  15.     return 0.0/0.0;
  16.  
  17.   double diff = log_x - log_y;
  18.  
  19.   // If log(x) dominates, return it...
  20.   if (diff > 708.0)
  21.     return log_x;
  22.  
  23.   /* We use the following trick:
  24.    * log(x-y) = log(a(x/a - y/a)) = log(a) + log(x/a - y/a)
  25.    *          = log(a) + log(exp(log(x)-log(a)) - exp(log(y)-log(a)))
  26.    * We pick log(a) = (log(x) + log(y)) / 2. So
  27.    * log(x-y) = (log(x) + log(y)) / 2 +
  28.    *          = log(exp((log(x) - log(y)) / 2) - exp((log(y) - log(x)) / 2))
  29.    */
  30.  
  31.   diff /= 2.0;
  32.  
  33.   return (log_x + log_y) / 2.0 + log(exp(diff) - exp(-diff));
  34. }

Report this snippet


Comments

RSS Icon Subscribe to comments

You need to login to post a comment.