#include <iostream>
 #include <string>
 #include <cstring>
 #include <cmath>
 #include <cstdlib>
 
 const double Pi = 3.14159265358979323846;
 const double ALMOST_ZERO = 9.0e-307;
 const double ALMOST_ONE  = 0.9999999999999997;
 
 
 std::string progname(const char* fullname)
 {
     const char* p;
     p = strrchr(fullname, '\\');
     if (p != 0)
     {
         p++;
     }
     else
     {
         p = strrchr(fullname, '/');
         if (p != 0)
         {
             p++;
         }
         else
         {
             p = fullname;
         }
     }
     std::string retval(p);
     return retval;
 }
 
 
 
 //--------------------------------------------------------------------
 //
 //  Name: FastNorm - Approximation for the area under the normal distribution
 //
 //  Syntax:         double FastNorm(double z);
 //  Description:    Area under the normal distribution to the left of z.
 //                  Adapted from Statistical computation,
 //                  J.H. Maindonald, 1984 Wiley & Sons, New York 1984.
 //                  (ISBN 0-471-86452-8)
 // Note:            This function should usually NOT be called by user
 //                  programs, because it is not very accurate for the central
 //                  part of the normal distribution (i.e. +/- ~4). Therefore,
 //                  user programs should normally use cumulative_normal_distribution()
 //                  unless speed is essential. However, in the tails of the distribution,
 //                  cumulative_normal_distribution() calls FastNorm().
 // ---------------------------------------------------------------------------------
 
 static double FastNorm(double Z)
 {   double X, M0, M1, M2, M3, M4, P;
 
     //  Tail calculation: Peizer & Pratt 1968
     //  Maindonald fig. 9.4 p. 292
 
     X = std::abs(Z);
     if (X < 1.9)
     {   M0 = 1.0/(1+0.33267*X);
         M1 = 0.4361836;
         M2 = -0.1201676;
         M3 = 0.937298;
         M4 = (X*X)/2.0;
         if (M4 <= ALMOST_ZERO) P = 0.5; // Avoid log(zero)
         else
         {   M4 = std::exp(M4) * std::sqrt(2.0*Pi);
             P  = M0*(M1+M2*M0+M3*std::sqrt(M0))/M4;
         }
     }
     else
     {   M1 = X*X;
         M2 = 1.0-1.0/(M1+3-1.0/(0.22*(M1+3.2)));
         M3 = std::log(M2) - std::log(X) - M1/2.0 - 0.9189386848;
         P  = std::exp(M3);
     }
     return (Z < 0.0) ? (P) : (1.0-P);
 }
 
 
 
 
 // ---------------------------------------------------------------------------------
 //
 //  Name:           cumulative_normal_distribution
 //  Purpose:        Approximation for the area under the normal distribution
 //  Syntax:         double cumulative_normal_distribution(double z);
 //  Description:    Calculates the area under the normal distribution to the left of z.
 //                  Adapted from Statistical computation, J.H. Maindonald,
 //                  1984 Wiley & Sons, New York 1984, ISBN 0-471-86452-8
 // ---------------------------------------------------------------------------------
 
 double cumulative_normal_distribution(double X)
 {   // Calculates area to the left of X, i.e. from (-infinity) to X.
     bool Negative_Argument;
     double P,S,C,I2;
     int I;
 
     Negative_Argument         = (X < 0.0);
     if (Negative_Argument)  X = -X; // calculate for positive argument
 
     if (fabs(X) < 4.1)
     {   // Morans approximation; fig 9.3 p 292, Maindonald
         S      = 0.0;
         C      = X * sqrt(2.0)/3.0;
         for (I = 0; I <= 12; I++)
         {   I2 = I+0.5;
             S  = S+sin(I2*C)*exp(-I2*I2/9)/I2;
         }
         P = 0.5 - S/Pi;
     }
     else P = FastNorm(-X);
     return (Negative_Argument) ? (P) : (1.0-P);
 }
 
 
 
 double d2(int n)
 {
     double x, dx;
     dx = 0.01;
     double sum = 0.0;
     for (x = -12; x <= 12; x = x + dx)
     {
         double alpha = cumulative_normal_distribution(x);
         sum += (1 - std::pow(1-alpha,n) - std::pow(alpha,n))*dx;
     }
     return sum;
 }
 
 
 double d3(int n)
 {
     double x_1, x_n, dx_1, dx_n;
     dx_1 = dx_n = 0.01;
     double sum = 0.0;
     for (x_n = -12; x_n <= 12; x_n = x_n + dx_n)
     {
         double alpha_n = cumulative_normal_distribution(x_n);
         for (x_1 = -12; x_1 <= x_n; x_1 = x_1 + dx_1)
         {
             double alpha_1 = cumulative_normal_distribution(x_1);
             // ERROR - mistake in original paper! sum += (1 - std::pow(alpha_1,n) - std::pow(1 - alpha_n,n) - std::pow(alpha_1 - alpha_n, n))*dx_1*dx_n;
             sum += (1 - std::pow(alpha_1,n) - std::pow(1 - alpha_n,n) + std::pow(alpha_1 - alpha_n, n))*dx_1*dx_n;
         }
     }
     double variance = 2*sum - std::pow(d2(n),2);
     return std::sqrt(variance);
 }
 
 
 
 main(int argc, char *argv[])
 {
     if (argc != 3)
     {
         std::cerr << "Syntax: " << progname(argv[0]) << " <option> " << "<integer>" << '\n';
         std::cerr << "where <option> is either -d2 or -d3, and <integer> is an integer greater than 2";
         exit(2);
     }
     std::string option(argv[1]);
     int integer = std::atoi(argv[2]);
 
     if (integer < 2)
     {
         std::cerr << "Error: the second argument should be an integer greater than 2." << '\n';
         exit(2);
     }
 
     double result(-1);
     if (option == "-d2")
     {
         result = d2(integer);
     }
     else if (option == "-d3")
     {
         result = d3(integer);
     }
     else
     {
         std::cerr << "Error: invalid option: " << option << '\n';
         exit(2);
     }
     std::cout << result;
     return 0;
 }