Compiler Patterns

/*
name:    expression_publics.c
purpose: set of functions normally contained in other files, but 
         collected here to help a "public" user. Prevents them 
         from having to have access to my whole source tree.
author:  David Qualls  (davidlqualls@yahoo.com)
date:    Nov 2011
license: CodeProject Open License (except lngamma_lanczos which 
         is GNU GPL.  Don't use that function unless you are 
         willing to propogate the license.  I've included the 
         lngamma_lanczos function here just in case you are 
         willing to propogate their license; saves you from 
         having to write your own).
*/

#include "expression_publics.h"


char *ltrim( const char* sp )
{
   while( *sp and isspace(*sp) )
   {
      sp++;
   }
   return (char*)sp;
}


int alphanumcmp( const char *s1, const char *s2 )
{
   int ret = 0;
   
   while( *s1 == *s2 )
   {
      if( not isalnum(*s1) )    // both strings hit the same delimiter,
         goto END;              // the strings are equal.
      s1++;
      s2++;
   }
// at this point, *s1 != *s2.  
   if( not isalnum(*s1) )       // s1 is terminated
   {  
      if( not isalnum(*s2) )    // and so is s2
         goto END;              // the strings are equal.
   // else                      // s2 is continuing
      ret = -1;                 // s2 is longer than s1.
   }
   else                         // s1 is continuing
   {
      if( not isalnum(*s2) )    // s2 is terminated
         ret = +1;              // s1 is longer than s2
      else                      // s2 is also continuing
      {
         if( *s1 < *s2 )
            ret = -1;
         else
            ret = +1;           // already know they're not equal.
      }
   }
END:
   return ret;
}


double mypow( double d1, double d2 )
{//if d2 is an integer, uses multiplication instead of transcendental 
 //calls.  This preserves the sign of d1. Plus, for small exponents, may be 
 //substantially faster than exp(log(x)).
 //
 //I wrote this before I realized that the standard pow() function may have pretty 
 //much the same functionality.  I had thought pow wouldn't work if d1 was < 0.
 //I'm leaving it in in-case some pow()s don't have this behavior.
 //
 //This func *may* be slightly non-standard as mypow(0,0) returns unity instead 
 //of raising an exception.  I think this is slightly preferred for some 
 //formulas that can be expressed more simply by assuming 0^0 == 1.
 
   double ret;
   long val = (long)d2;  //truncate the exponent.
   
   if( val == d2 ) //if exponent is an integer,
   {
      int oneover;
      unsigned long expon;
      double oppow; //oppow is d1 raised to an integer power of 2: d1^(2^i)
      
      if( val < 0 )
      {  oneover = 1;
         expon = (unsigned long) -val;
      }
      else
      {  oneover = 0;
         expon = (unsigned long)val;
      }
         
      oppow = d1; // d1^(2^0)
      ret = 1;
      while( expon )
      {
         if( expon bitand 1 ) //getting ready to shift a '1' off the right end of expon.
            ret *= oppow;     //then include it in the answer.
         oppow *= oppow;      // oppow = oppow^2
         expon >>= 1;
      }
      if( oneover )  //exponent was negative,
         ret = 1.0/ret;
   }
   else if( 0.0 == d1 ) //can't take log(0)
      ret = 0.0;
   else
      ret =  exp(log(d1)*d2);
   return ret;
}

#define pi       3.1415926535897932384626433832795
#define sqrt_pi  1.7724538509055160272981674833411
#define piover2  1.5707963267948966192313216916398
#define piover4  0.78539816339744830961566084581988
#define LogRootTwoPi_ 0.91893853320467274178032973640562
#define e        2.71828182845904523536028747135266
#define euler    0.57721566490153286060651209008240
#define sqrt_2   1.4142135623730950488016887242097
#define sqrt_3   1.7320508075688772935274463415059

const static double gp[] = //for use with gamma func.
{0.99999999999980993, 676.5203681218851, -1259.1392167224028,
 771.32342877765313, -176.61502916214059, 12.507343278686905,
 -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 
};
double gamma( double z )
{/*returns the gamma of REAL z.
   Gamma(z) = integrate( t^(z-1) * exp(-t), 'dt', 0, infinity )
   Adapted from en.wikipedia.org/wiki/Lanczos_approximation.
   I added the code to deal with near-zero conditions.
 */
    
   double    x, t;
   int       ndx;
   const int g = lengthof(gp)-2;
   
   if( z < 0.0 )
      return pi / sin(pi*z) * gamma(1-z) ;
   if( z < 0.5 ) //the approximation seems to struggle near zero.
      return gamma(z+3)/(z*(z+1)*(z+2));
   z -= 1.0;
   x = gp[0];
   for( ndx=1; ndx<g+2; ndx++ )
   {
      x += gp[ndx] / (z+ndx) ;
   } 
   t = z + g + 0.5;
//   return sqrt_2*sqrt_pi * pow(t, z+0.5) * exp(-t) * x;
   return sqrt_2*sqrt_pi * exp((z+0.5)*log(t) -t) * x;
}



static double lanczos_7_c[9] = {
  0.99999999999980993227684700473478,
  676.520368121885098567009190444019,
 -1259.13921672240287047156078755283,
  771.3234287776530788486528258894,
 -176.61502916214059906584551354,
  12.507343278686904814458936853,
 -0.13857109526572011689554707,
  9.984369578019570859563e-6,
  1.50563273514931155834e-7
};

static double lngamma_lanczos(double x)
{//copied (adapted) from GNU Scientific Library (GSL).  
 //CAUTION! You may not have the rights to use this!  Be 
 //sure to check the GNU General Public License 
 //at http://www.gnu.org/licenses/
  int k;
  double Ag;
  double term1, term2;
  
  x -= 1.0; /* Lanczos writes z! instead of Gamma(z) */

  Ag = lanczos_7_c[0];
  for(k=1; k<=8; k++) { Ag += lanczos_7_c[k]/(x+k); }

  /* (x+0.5)*log(x+7.5) - (x+7.5) + LogRootTwoPi_ + log(Ag(x)) */
  term1 = (x+0.5)*log((x+7.5)/M_E);
  term2 = LogRootTwoPi_ + log(Ag);
  return term1 + (term2 - 7.0);
}


double lngamma( double z )
{
   if( z <= 0.0 )
   {
      RaiseD( ArgumentOutOfRange, z );
      return 0.0;
   }

   if( z < 2.1 ) //need to get argument away from 1.0 or 2.0
      return lngamma_lanczos( z+3 ) - log( z*(z+1)*(z+2) );
   
   return lngamma_lanczos( z );
}


const static double fp[] = //for use with the factorial func
{// 0    1    2    3    4      5      6      7        8
   1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0
 //   9         10          11          12            13
 , 362880.0, 3628800.0, 39916800.0, 479001600.0, 6227020800.0
 //    14               15              16               17
 , 87178291200.0, 1307674368.0e3, 20922789888.0e3, 355687428096.0e3
 //       18                 19
 , 6402373705728.0e3, 121645100408832.0e3
};
double factorial( unsigned j )
{
   const unsigned g = lengthof(fp);
   double ret;
   
   if( j < g )
      ret = fp[j];
   else if( j < (g*7)>>2 ) //1.75*g
   {
      ret = j;
      while( --j >= g )
         ret *= j;
      ret *= fp[j];
   }
   else
      ret = gamma( j + 1.0 );
   
   return ret;
}



#ifndef __DMC__
//The Digital Mars Compiler library already includes error function and complementary.
//Perhaps because it's C99 compliant.  C99 contains these funcs.  If included within 
//the compiler, probably better to use it as the following approximation may not have 
//enough accuracy.

double erfc( double x )
{/*returns the complementary error function defined by
   erfc(x) = 2/sqrt_pi * integrate( exp(-t^2), 'dt', x, infinity )
   This is the area under the standard normal curve, to the right of x.
   Effectively, 1 - CDF(x) where CDF is the cumulative distribution func
   for the standard normal.  See random.c for putting this back to non-
   standard (i.e. not 0 == mean and 1 == sigma).
*/
   double t, z, ret;
   
   z = fabs(x);
   t = 1.0 / (1.0 + z/2) ;
   ret = t*exp( -z*z   //uses an embedded 10'th order poly and horner's rule.
                 -1.26551223 + t*(
                  1.00002368 + t*(
                  0.37409196 + t*(
                  0.09678418 + t*(
                 -0.18628806 + t*(
                  0.27886807 + t*(
                 -1.13520398 + t*(
                  1.48851587 + t*(
                 -0.82215223 + t*(
                  0.17087277)))))))))) ;
   if( x >= 0.0 )
      return ret;
   return 2.0 - ret;
}

double erf( double x )
{/*returns the error function defined by
   erfc(x) = 2/sqrt_pi * integrate( exp(-t^2), 'dt', 0.0, x )
   This is used in computing the CDF for a standard normal curve.
*/
   return 1.0 - erfc( x );
}

#endif
//end of #ifndef __DMC__