The SQLite Database Engine

Listing 3: Creating a custom function.

void normal_curve_area(sqlite_func* context, int argc, const char **argv)
{
   char** endptr;
   char result[65];
   double x1, x2, mu, sigma;
   if(argc != 4)
   {
      return;
   }
   x1 = strtod((char*)argv[0], endptr);
   if((x1==0) && (argv[0]==*endptr))
      return;
   x2 = strtod((char*)argv[1], endptr);
   if((x2==0) && (argv[1]==*endptr))

      return;
   mu = strtod((char*)argv[2], endptr);
   if((x1==0) && (argv[2]==*endptr))
      return;
   sigma = strtod((char*)argv[3], endptr);
   if((x1==0) && (argv[3]==*endptr))
      return;
   sprintf(result, "%f", GetNormalCurveArea(x1,x2,mu,sigma));
   sqlite_set_result_string(context, result, -1);
}
double GetNormalCurveArea(double x1, double x2, double mu, double sigma)
{
   /* Maclaurin Series Expansion for exp(-x2/2)
   Michael Owens
   Description: This function takes two random variables, a lower 
   limit (x1) and an upper limit (x2), on a Gaussian distribution and 
   computes the total area between them.
   User Input Parameters: 
   x2: upper limit
   x1: lower limit
   mu: population mean
   sigma: variance
   Nomenclature:
   sz: dummy variable for series expansion
   z = (x-mu)/sig
   cum: the cumulative value of z, or integral
   cum1 is the area from -r1 to 0 while
   cum2 is the area from 0 to r2.
   Limitations:
   The Limiting Values of z: A value greater than z=5 will give exactly 50% of
   the normal curve to four decimal places, and larger values will only
   encumber series convergence, therefore any values greater than 4 will be
   reset to 4.
   */
   double j = 10; // Initialized for priming the while() block.
   double bound = 4.2;
   double z1 = (x1 - mu) / sigma;
   double z2 = (x2 - mu) / sigma;
   if (z1 < -bound)
      z1 = (double)-bound;
   if (z1 > bound)
      z1 = (double)bound;
   if (z2 < -bound)
      z2 = (double)-bound;
   if (z2 > bound)
      z2 = (double)bound;
   double cum1 = fabs(z1);
   double cum2 = fabs(z2);
   // Use absolute values for computing terms
   x1 = fabs(z1);
   x2 = fabs(z2);
   // Computations
   // Maclaurin Series: term by term addition
   // Area of lower limit
   if(cum1)
      SeriesExpansion(x1,cum1);
   else
      cum1 = 0;
   // Area of upper limit
   if(cum2)
      SeriesExpansion(x2,cum2);
   else
      cum2 = 0;
   // Determine the total area:
   double Area;
   if ((z2 + z2) < (fabs(z2 + z2))) // if z2 is negative
      Area = cum1 - cum2; // then z1 must be negative too.
   else
      if ((z1 + z1) < (fabs(z1 + z1))) // z2 is positve and if z1 negative
         Area = cum1 + cum2;
      else
         Area = fabs(cum2 - cum1); // if z1 is positive
   // Limiting area from origin to +infinity
   double CA;
   CA = pow(2*3.1415926535, 0.5);
   // Normalized area
   Area = Area/CA; // Area from origin to lower limit.
   return Area;
}
short SeriesExpansion(double &x, double &cum)
{
   double SeriesTerm;
   double j = 10;
   for (int i = 1; j > 0.0001; i++)
   {
      int f = i;
      double factorial = f;
      if(f-1)
      {
      while(f-1)
         factorial *= --f;
      }
      if(!factorial)
      return 0;
      SeriesTerm = (pow(-1,i));
      SeriesTerm *= (double)1/((2*i)+1);
      SeriesTerm *= (double)pow(x,(2*i+1));
      SeriesTerm *= (double)1/((pow(2,i)*factorial));
      cum += SeriesTerm;
      j = fabs(SeriesTerm);
   }
   return 1;
}