// @(#)root/mathmore:$Id$
// Authors: B. List 29.4.2010
 

 /**********************************************************************
  *                                                                    *
  * Copyright (c) 2004 ROOT Foundation,  CERN/PH-SFT                   *
  *                                                                    *
  * This library is free software; you can redistribute it and/or      *
  * modify it under the terms of the GNU General Public License        *
  * as published by the Free Software Foundation; either version 2     *
  * of the License, or (at your option) any later version.             *
  *                                                                    *
  * This library is distributed in the hope that it will be useful,    *
  * but WITHOUT ANY WARRANTY; without even the implied warranty of     *
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU   *
  * General Public License for more details.                           *
  *                                                                    *
  * You should have received a copy of the GNU General Public License  *
  * along with this library (see file COPYING); if not, write          *
  * to the Free Software Foundation, Inc., 59 Temple Place, Suite      *
  * 330, Boston, MA 02111-1307 USA, or contact the author.             *
  *                                                                    *
  **********************************************************************/

// Implementation file for class VavilovTest
// 
// Created by: blist  at Thu Apr 29 11:19:00 2010
// 
// Last update: Thu Apr 29 11:19:00 2010
// 


#include "VavilovTest.h"
#include "Math/Vavilov.h"
#include "Math/VavilovAccurate.h"
#include "Math/VavilovFast.h"
//#include "Math/SpecFuncMathCore.h"
//#include "Math/SpecFuncMathMore.h"

#include <cassert>
#include <iostream>
#include <cmath>
#include <iomanip>
#include <cmath>
#include <cstdlib>
#include <string>
#include <sstream>


namespace ROOT {
namespace Math {

static const double eu = 0.577215664901532860606;      // Euler's constant


//kappa = 0.01
static double vavilovPdfValues0[45][12]={
   {-3.00, 6.80E-4, 6.85E-4, 6.90E-4, 6.95E-4, 7.00E-4, 7.05E-4, 7.10E-4, 7.15E-4, 7.21E-4, 7.26E-4, 7.31E-4},               
   {-2.50, 9.73E-3, 9.80E-3, 9.87E-3, 9.93E-3, 1.00E-2, 1.01E-2, 1.01E-2, 1.02E-2, 1.03E-2, 1.03E-2, 1.04E-2},               
   {-2.00, 4.44E-2, 4.47E-2, 4.50E-2, 4.53E-2, 4.55E-2, 4.58E-2, 4.61E-2, 4.64E-2, 4.67E-2, 4.70E-2, 4.73E-2},               
   {-1.50, 1.02E-1, 1.02E-1, 1.03E-1, 1.03E-1, 1.04E-1, 1.04E-1, 1.05E-1, 1.06E-1, 1.06E-1, 1.07E-1, 1.07E-1},               
   {-1.00, 1.53E-1, 1.54E-1, 1.55E-1, 1.55E-1, 1.56E-1, 1.57E-1, 1.58E-1, 1.59E-1, 1.59E-1, 1.60E-1, 1.61E-1},               
   {-0.50, 1.79E-1, 1.80E-1, 1.81E-1, 1.82E-1, 1.83E-1, 1.83E-1, 1.84E-1, 1.85E-1, 1.86E-1, 1.87E-1, 1.88E-1},               
   { 0.00, 1.81E-1, 1.81E-1, 1.82E-1, 1.83E-1, 1.84E-1, 1.84E-1, 1.85E-1, 1.86E-1, 1.87E-1, 1.88E-1, 1.88E-1},               
   { 0.50, 1.67E-1, 1.68E-1, 1.68E-1, 1.69E-1, 1.69E-1, 1.70E-1, 1.71E-1, 1.71E-1, 1.72E-1, 1.73E-1, 1.73E-1},               
   { 1.00, 1.47E-1, 1.47E-1, 1.48E-1, 1.48E-1, 1.49E-1, 1.49E-1, 1.49E-1, 1.50E-1, 1.50E-1, 1.51E-1, 1.51E-1},               
   { 1.50, 1.25E-1, 1.26E-1, 1.26E-1, 1.26E-1, 1.27E-1, 1.27E-1, 1.27E-1, 1.28E-1, 1.28E-1, 1.29E-1, 1.29E-1},               
   { 2.00, 1.06E-1, 1.06E-1, 1.06E-1, 1.07E-1, 1.07E-1, 1.07E-1, 1.07E-1, 1.08E-1, 1.08E-1, 1.08E-1, 1.08E-1},               
   { 2.50, 8.91E-2, 8.93E-2, 8.94E-2, 8.96E-2, 8.97E-2, 8.99E-2, 9.00E-2, 9.02E-2, 9.03E-2, 9.04E-2, 9.06E-2},               
   { 3.00, 7.50E-2, 7.51E-2, 7.52E-2, 7.53E-2, 7.53E-2, 7.54E-2, 7.55E-2, 7.56E-2, 7.57E-2, 7.58E-2, 7.59E-2},               
   { 3.50, 6.34E-2, 6.34E-2, 6.34E-2, 6.35E-2, 6.35E-2, 6.36E-2, 6.36E-2, 6.36E-2, 6.37E-2, 6.37E-2, 6.37E-2},               
   { 4.00, 5.38E-2, 5.38E-2, 5.38E-2, 5.38E-2, 5.38E-2, 5.38E-2, 5.38E-2, 5.39E-2, 5.39E-2, 5.39E-2, 5.39E-2},               
   { 4.50, 4.60E-2, 4.60E-2, 4.60E-2, 4.59E-2, 4.59E-2, 4.59E-2, 4.59E-2, 4.59E-2, 4.58E-2, 4.58E-2, 4.58E-2},               
   { 5.00, 3.96E-2, 3.95E-2, 3.95E-2, 3.95E-2, 3.94E-2, 3.94E-2, 3.93E-2, 3.93E-2, 3.93E-2, 3.92E-2, 3.92E-2},               
   { 5.50, 3.43E-2, 3.42E-2, 3.42E-2, 3.41E-2, 3.41E-2, 3.40E-2, 3.40E-2, 3.39E-2, 3.39E-2, 3.38E-2, 3.38E-2},               
   { 6.00, 2.99E-2, 2.98E-2, 2.97E-2, 2.97E-2, 2.96E-2, 2.96E-2, 2.95E-2, 2.95E-2, 2.94E-2, 2.93E-2, 2.93E-2},               
   { 6.50, 2.62E-2, 2.61E-2, 2.61E-2, 2.60E-2, 2.59E-2, 2.59E-2, 2.58E-2, 2.57E-2, 2.57E-2, 2.56E-2, 2.55E-2},               
   { 7.00, 2.31E-2, 2.30E-2, 2.30E-2, 2.29E-2, 2.28E-2, 2.28E-2, 2.27E-2, 2.26E-2, 2.26E-2, 2.25E-2, 2.24E-2},               
   { 7.50, 2.05E-2, 2.04E-2, 2.04E-2, 2.03E-2, 2.02E-2, 2.01E-2, 2.01E-2, 2.00E-2, 1.99E-2, 1.99E-2, 1.98E-2},               
   { 8.00, 1.83E-2, 1.82E-2, 1.81E-2, 1.81E-2, 1.80E-2, 1.79E-2, 1.78E-2, 1.78E-2, 1.77E-2, 1.76E-2, 1.75E-2},               
   { 8.50, 1.64E-2, 1.63E-2, 1.62E-2, 1.62E-2, 1.61E-2, 1.60E-2, 1.59E-2, 1.59E-2, 1.58E-2, 1.57E-2, 1.56E-2},               
   { 9.00, 1.47E-2, 1.47E-2, 1.46E-2, 1.45E-2, 1.45E-2, 1.44E-2, 1.43E-2, 1.42E-2, 1.42E-2, 1.41E-2, 1.40E-2},               
   { 9.50, 1.33E-2, 1.33E-2, 1.32E-2, 1.31E-2, 1.30E-2, 1.30E-2, 1.29E-2, 1.28E-2, 1.27E-2, 1.27E-2, 1.26E-2},               
   {10.00, 1.21E-2, 1.20E-2, 1.20E-2, 1.19E-2, 1.18E-2, 1.17E-2, 1.17E-2, 1.16E-2, 1.15E-2, 1.14E-2, 1.14E-2},               
   {11.00, 1.01E-2, 1.00E-2, 9.94E-3, 9.87E-3, 9.80E-3, 9.73E-3, 9.66E-3, 9.59E-3, 9.51E-3, 9.44E-3, 9.37E-3},               
   {12.00, 8.51E-3, 8.44E-3, 8.37E-3, 8.31E-3, 8.24E-3, 8.17E-3, 8.10E-3, 8.03E-3, 7.96E-3, 7.89E-3, 7.82E-3},               
   {13.00, 7.26E-3, 7.20E-3, 7.14E-3, 7.07E-3, 7.00E-3, 6.94E-3, 6.87E-3, 6.81E-3, 6.74E-3, 6.67E-3, 6.61E-3},               
   {14.00, 6.27E-3, 6.20E-3, 6.14E-3, 6.08E-3, 6.02E-3, 5.95E-3, 5.89E-3, 5.83E-3, 5.76E-3, 5.70E-3, 5.64E-3},               
   {15.00, 5.46E-3, 5.40E-3, 5.34E-3, 5.28E-3, 5.22E-3, 5.16E-3, 5.10E-3, 5.04E-3, 4.97E-3, 4.91E-3, 4.85E-3},               
   {16.00, 4.79E-3, 4.73E-3, 4.67E-3, 4.62E-3, 4.56E-3, 4.50E-3, 4.44E-3, 4.39E-3, 4.33E-3, 4.27E-3, 4.21E-3},               
   {17.00, 4.23E-3, 4.18E-3, 4.12E-3, 4.07E-3, 4.01E-3, 3.96E-3, 3.90E-3, 3.85E-3, 3.79E-3, 3.73E-3, 3.68E-3},               
   {18.00, 3.77E-3, 3.72E-3, 3.66E-3, 3.61E-3, 3.56E-3, 3.50E-3, 3.45E-3, 3.40E-3, 3.34E-3, 3.29E-3, 3.24E-3},               
   {19.00, 3.37E-3, 3.32E-3, 3.27E-3, 3.22E-3, 3.17E-3, 3.12E-3, 3.07E-3, 3.02E-3, 2.97E-3, 2.91E-3, 2.86E-3},               
   {20.00, 3.04E-3, 2.99E-3, 2.94E-3, 2.89E-3, 2.84E-3, 2.79E-3, 2.74E-3, 2.69E-3, 2.64E-3, 2.60E-3, 2.55E-3},               
   {22.00, 2.49E-3, 2.45E-3, 2.40E-3, 2.36E-3, 2.31E-3, 2.27E-3, 2.22E-3, 2.18E-3, 2.13E-3, 2.09E-3, 2.04E-3},               
   {24.00, 2.08E-3, 2.04E-3, 2.00E-3, 1.96E-3, 1.92E-3, 1.88E-3, 1.83E-3, 1.79E-3, 1.75E-3, 1.71E-3, 1.66E-3},               
   {26.00, 1.77E-3, 1.73E-3, 1.69E-3, 1.65E-3, 1.61E-3, 1.57E-3, 1.53E-3, 1.49E-3, 1.45E-3, 1.41E-3, 1.37E-3},               
   {28.00, 1.51E-3, 1.48E-3, 1.44E-3, 1.41E-3, 1.37E-3, 1.33E-3, 1.30E-3, 1.26E-3, 1.22E-3, 1.18E-3, 1.15E-3},               
   {30.00, 1.31E-3, 1.28E-3, 1.24E-3, 1.21E-3, 1.18E-3, 1.14E-3, 1.11E-3, 1.07E-3, 1.04E-3, 1.00E-3, 9.68E-4},               
   {32.00, 1.15E-3, 1.12E-3, 1.08E-3, 1.05E-3, 1.02E-3, 9.87E-4, 9.54E-4, 9.22E-4, 8.89E-4, 8.56E-4, 8.24E-4},               
   {34.00, 1.01E-3, 9.81E-4, 9.51E-4, 9.21E-4, 8.90E-4, 8.60E-4, 8.29E-4, 7.98E-4, 7.68E-4, 7.37E-4, 7.06E-4},               
   {36.00, 8.98E-4, 8.70E-4, 8.41E-4, 8.12E-4, 7.83E-4, 7.54E-4, 7.25E-4, 6.96E-4, 6.67E-4, 6.38E-4, 6.09E-4}};         


//kappa = 0.04
static double vavilovPdfValues1[42][12]={
   {-3.00, 7.01E-4, 7.18E-4, 7.34E-4, 7.52E-4, 7.69E-4, 7.87E-4, 8.06E-4, 8.25E-4, 8.44E-4, 8.64E-4, 8.84E-4},
   {-2.50, 1.00E-2, 1.02E-2, 1.05E-2, 1.07E-2, 1.09E-2, 1.12E-2, 1.14E-2, 1.16E-2, 1.19E-2, 1.21E-2, 1.24E-2},
   {-2.00, 4.58E-2, 4.67E-2, 4.76E-2, 4.85E-2, 4.94E-2, 5.04E-2, 5.13E-2, 5.23E-2, 5.33E-2, 5.44E-2, 5.54E-2},
   {-1.50, 1.05E-1, 1.06E-1, 1.08E-1, 1.10E-1, 1.12E-1, 1.14E-1, 1.16E-1, 1.18E-1, 1.20E-1, 1.22E-1, 1.24E-1},
   {-1.00, 1.58E-1, 1.60E-1, 1.62E-1, 1.65E-1, 1.67E-1, 1.70E-1, 1.73E-1, 1.75E-1, 1.78E-1, 1.81E-1, 1.83E-1},
   {-0.50, 1.85E-1, 1.87E-1, 1.89E-1, 1.92E-1, 1.95E-1, 1.97E-1, 2.00E-1, 2.02E-1, 2.05E-1, 2.08E-1, 2.10E-1},
   { 0.00, 1.86E-1, 1.88E-1, 1.90E-1, 1.92E-1, 1.95E-1, 1.97E-1, 1.99E-1, 2.01E-1, 2.03E-1, 2.06E-1, 2.08E-1},
   { 0.50, 1.72E-1, 1.74E-1, 1.75E-1, 1.77E-1, 1.78E-1, 1.80E-1, 1.82E-1, 1.83E-1, 1.85E-1, 1.87E-1, 1.88E-1},
   { 1.00, 1.51E-1, 1.52E-1, 1.53E-1, 1.54E-1, 1.55E-1, 1.57E-1, 1.58E-1, 1.59E-1, 1.60E-1, 1.61E-1, 1.62E-1},
   { 1.50, 1.29E-1, 1.30E-1, 1.31E-1, 1.31E-1, 1.32E-1, 1.33E-1, 1.33E-1, 1.34E-1, 1.34E-1, 1.35E-1, 1.36E-1},
   { 2.00, 1.09E-1, 1.10E-1, 1.10E-1, 1.10E-1, 1.11E-1, 1.11E-1, 1.11E-1, 1.11E-1, 1.12E-1, 1.12E-1, 1.12E-1},
   { 2.50, 9.18E-2, 9.19E-2, 9.20E-2, 9.21E-2, 9.22E-2, 9.22E-2, 9.23E-2, 9.23E-2, 9.23E-2, 9.24E-2, 9.24E-2},
   { 3.00, 7.73E-2, 7.72E-2, 7.71E-2, 7.70E-2, 7.69E-2, 7.68E-2, 7.67E-2, 7.66E-2, 7.64E-2, 7.62E-2, 7.61E-2},
   { 3.50, 6.53E-2, 6.51E-2, 6.49E-2, 6.47E-2, 6.45E-2, 6.42E-2, 6.40E-2, 6.37E-2, 6.34E-2, 6.32E-2, 6.29E-2},
   { 4.00, 5.54E-2, 5.52E-2, 5.49E-2, 5.46E-2, 5.43E-2, 5.40E-2, 5.36E-2, 5.33E-2, 5.29E-2, 5.26E-2, 5.22E-2},
   { 4.50, 4.74E-2, 4.71E-2, 4.67E-2, 4.64E-2, 4.60E-2, 4.56E-2, 4.53E-2, 4.49E-2, 4.45E-2, 4.40E-2, 4.36E-2},
   { 5.00, 4.08E-2, 4.04E-2, 4.00E-2, 3.96E-2, 3.92E-2, 3.88E-2, 3.84E-2, 3.80E-2, 3.76E-2, 3.71E-2, 3.67E-2},
   { 5.50, 3.53E-2, 3.49E-2, 3.45E-2, 3.41E-2, 3.37E-2, 3.33E-2, 3.28E-2, 3.24E-2, 3.19E-2, 3.15E-2, 3.10E-2},
   { 6.00, 3.08E-2, 3.04E-2, 3.00E-2, 2.95E-2, 2.91E-2, 2.87E-2, 2.82E-2, 2.78E-2, 2.73E-2, 2.69E-2, 2.64E-2},
   { 6.50, 2.70E-2, 2.66E-2, 2.62E-2, 2.57E-2, 2.53E-2, 2.49E-2, 2.44E-2, 2.40E-2, 2.35E-2, 2.31E-2, 2.26E-2},
   { 7.00, 2.38E-2, 2.34E-2, 2.30E-2, 2.26E-2, 2.21E-2, 2.17E-2, 2.13E-2, 2.08E-2, 2.04E-2, 1.99E-2, 1.94E-2},
   { 7.50, 2.11E-2, 2.07E-2, 2.03E-2, 1.99E-2, 1.95E-2, 1.90E-2, 1.86E-2, 1.82E-2, 1.77E-2, 1.73E-2, 1.68E-2},
   { 8.00, 1.88E-2, 1.84E-2, 1.80E-2, 1.76E-2, 1.72E-2, 1.68E-2, 1.64E-2, 1.59E-2, 1.55E-2, 1.51E-2, 1.46E-2},
   { 8.50, 1.69E-2, 1.65E-2, 1.61E-2, 1.57E-2, 1.53E-2, 1.49E-2, 1.45E-2, 1.40E-2, 1.36E-2, 1.32E-2, 1.27E-2},
   { 9.00, 1.52E-2, 1.48E-2, 1.44E-2, 1.40E-2, 1.36E-2, 1.32E-2, 1.28E-2, 1.24E-2, 1.20E-2, 1.16E-2, 1.12E-2},
   { 9.50, 1.37E-2, 1.34E-2, 1.30E-2, 1.26E-2, 1.22E-2, 1.18E-2, 1.14E-2, 1.10E-2, 1.06E-2, 1.02E-2, 9.81E-3},
   {10.00, 1.25E-2, 1.21E-2, 1.17E-2, 1.14E-2, 1.10E-2, 1.06E-2, 1.02E-2, 9.84E-3, 9.45E-3, 9.05E-3, 8.65E-3},
   {11.00, 1.04E-2, 1.00E-2, 9.70E-3, 9.35E-3, 8.99E-3, 8.63E-3, 8.27E-3, 7.91E-3, 7.54E-3, 7.17E-3, 6.79E-3},
   {12.00, 8.77E-3, 8.44E-3, 8.11E-3, 7.78E-3, 7.45E-3, 7.11E-3, 6.77E-3, 6.43E-3, 6.08E-3, 5.73E-3, 5.38E-3},
   {13.00, 7.49E-3, 7.18E-3, 6.87E-3, 6.56E-3, 6.24E-3, 5.92E-3, 5.60E-3, 5.28E-3, 4.95E-3, 4.63E-3, 4.30E-3},
   {14.00, 6.46E-3, 6.17E-3, 5.87E-3, 5.58E-3, 5.28E-3, 4.98E-3, 4.68E-3, 4.38E-3, 4.07E-3, 3.76E-3, 3.45E-3},
   {15.00, 5.62E-3, 5.35E-3, 5.07E-3, 4.79E-3, 4.51E-3, 4.22E-3, 3.94E-3, 3.65E-3, 3.36E-3, 3.07E-3, 2.78E-3},
   {16.00, 4.93E-3, 4.67E-3, 4.41E-3, 4.14E-3, 3.88E-3, 3.61E-3, 3.34E-3, 3.07E-3, 2.80E-3, 2.52E-3, 2.24E-3},
   {17.00, 4.36E-3, 4.11E-3, 3.86E-3, 3.61E-3, 3.36E-3, 3.10E-3, 2.85E-3, 2.59E-3, 2.33E-3, 2.07E-3, 1.81E-3},
   {18.00, 3.88E-3, 3.65E-3, 3.41E-3, 3.17E-3, 2.93E-3, 2.69E-3, 2.44E-3, 2.20E-3, 1.95E-3, 1.71E-3, 1.46E-3},
   {19.00, 3.48E-3, 3.25E-3, 3.02E-3, 2.79E-3, 2.57E-3, 2.34E-3, 2.10E-3, 1.87E-3, 1.64E-3, 1.41E-3, 1.17E-3},
   {20.00, 3.13E-3, 2.91E-3, 2.70E-3, 2.48E-3, 2.26E-3, 2.04E-3, 1.82E-3, 1.60E-3, 1.38E-3, 1.16E-3, 9.32E-4},
   {22.00, 2.57E-3, 2.37E-3, 2.17E-3, 1.98E-3, 1.78E-3, 1.58E-3, 1.37E-3, 1.17E-3, 9.71E-4, 7.69E-4, 5.66E-4},
   {24.00, 1.97E-3, 1.80E-3, 1.63E-3, 1.47E-3, 1.30E-3, 1.13E-3, 9.69E-4, 8.04E-4, 6.39E-4, 4.75E-4, 3.11E-4},
   {26.00, 1.14E-3, 1.04E-3, 9.34E-4, 8.32E-4, 7.31E-4, 6.31E-4, 5.34E-4, 4.38E-4, 3.44E-4, 2.52E-4, 1.61E-4},
   {28.00, 6.50E-4, 5.85E-4, 5.22E-4, 4.61E-4, 4.02E-4, 3.45E-4, 2.89E-4, 2.35E-4, 1.84E-4, 1.34E-4, 8.60E-5},
   {30.00, 3.95E-4, 3.53E-4, 3.12E-4, 2.73E-4, 2.36E-4, 2.00E-4, 1.66E-4, 1.34E-4, 1.03E-4, 7.46E-5, 4.76E-5}};


//kappa = 0.07
static double vavilovPdfValues2[41][12]={
   {-3.50, 0,       0,       0,       0,       0,       0,       1.01E-5, 1.06E-5, 1.10E-5, 1.14E-5, 1.19E-5},
   {-3.00, 7.23E-4, 7.50E-4, 7.78E-4, 8.07E-4, 8.37E-4, 8.68E-4, 9.00E-4, 9.34E-4, 9.69E-4, 1.01E-3, 1.04E-3},
   {-2.50, 1.03E-2, 1.07E-2, 1.10E-2, 1.14E-2, 1.18E-2, 1.22E-2, 1.26E-2, 1.30E-2, 1.35E-2, 1.39E-2, 1.44E-2},
   {-2.00, 4.72E-2, 4.86E-2, 5.00E-2, 5.15E-2, 5.31E-2, 5.47E-2, 5.63E-2, 5.80E-2, 5.98E-2, 6.15E-2, 6.34E-2},
   {-1.50, 1.08E-1, 1.11E-1, 1.14E-1, 1.17E-1, 1.20E-1, 1.23E-1, 1.26E-1, 1.29E-1, 1.33E-1, 1.36E-1, 1.40E-1},
   {-1.00, 1.62E-1, 1.66E-1, 1.70E-1, 1.74E-1, 1.78E-1, 1.82E-1, 1.86E-1, 1.90E-1, 1.94E-1, 1.99E-1, 2.03E-1},
   {-0.50, 1.90E-1, 1.94E-1, 1.98E-1, 2.01E-1, 2.05E-1, 2.09E-1, 2.13E-1, 2.17E-1, 2.21E-1, 2.25E-1, 2.30E-1},
   { 0.00, 1.92E-1, 1.95E-1, 1.98E-1, 2.01E-1, 2.04E-1, 2.07E-1, 2.10E-1, 2.14E-1, 2.17E-1, 2.20E-1, 2.23E-1},
   { 0.50, 1.77E-1, 1.79E-1, 1.82E-1, 1.84E-1, 1.86E-1, 1.88E-1, 1.90E-1, 1.92E-1, 1.95E-1, 1.97E-1, 1.99E-1},
   { 1.00, 1.56E-1, 1.57E-1, 1.58E-1, 1.60E-1, 1.61E-1, 1.62E-1, 1.64E-1, 1.65E-1, 1.66E-1, 1.67E-1, 1.68E-1},
   { 1.50, 1.33E-1, 1.34E-1, 1.35E-1, 1.35E-1, 1.36E-1, 1.36E-1, 1.37E-1, 1.37E-1, 1.38E-1, 1.38E-1, 1.39E-1},
   { 2.00, 1.13E-1, 1.13E-1, 1.13E-1, 1.13E-1, 1.13E-1, 1.13E-1, 1.13E-1, 1.13E-1, 1.13E-1, 1.13E-1, 1.13E-1},
   { 2.50, 9.46E-2, 9.44E-2, 9.42E-2, 9.40E-2, 9.37E-2, 9.33E-2, 9.30E-2, 9.26E-2, 9.21E-2, 9.16E-2, 9.11E-2},
   { 3.00, 7.96E-2, 7.92E-2, 7.87E-2, 7.82E-2, 7.77E-2, 7.71E-2, 7.65E-2, 7.58E-2, 7.51E-2, 7.44E-2, 7.36E-2},
   { 3.50, 6.73E-2, 6.67E-2, 6.60E-2, 6.53E-2, 6.46E-2, 6.39E-2, 6.31E-2, 6.23E-2, 6.14E-2, 6.06E-2, 5.96E-2},
   { 4.00, 5.71E-2, 5.64E-2, 5.57E-2, 5.49E-2, 5.41E-2, 5.32E-2, 5.23E-2, 5.14E-2, 5.05E-2, 4.95E-2, 4.85E-2},
   { 4.50, 4.88E-2, 4.80E-2, 4.72E-2, 4.64E-2, 4.55E-2, 4.46E-2, 4.37E-2, 4.27E-2, 4.17E-2, 4.07E-2, 3.97E-2},
   { 5.00, 4.20E-2, 4.12E-2, 4.03E-2, 3.95E-2, 3.86E-2, 3.76E-2, 3.67E-2, 3.57E-2, 3.47E-2, 3.36E-2, 3.26E-2},
   { 5.50, 3.64E-2, 3.55E-2, 3.47E-2, 3.38E-2, 3.29E-2, 3.19E-2, 3.10E-2, 3.00E-2, 2.90E-2, 2.80E-2, 2.69E-2},
   { 6.00, 3.17E-2, 3.09E-2, 3.00E-2, 2.91E-2, 2.82E-2, 2.73E-2, 2.63E-2, 2.53E-2, 2.44E-2, 2.33E-2, 2.23E-2},
   { 6.50, 2.78E-2, 2.70E-2, 2.61E-2, 2.52E-2, 2.43E-2, 2.34E-2, 2.25E-2, 2.15E-2, 2.06E-2, 1.96E-2, 1.86E-2},
   { 7.00, 2.45E-2, 2.37E-2, 2.29E-2, 2.20E-2, 2.11E-2, 2.02E-2, 1.93E-2, 1.84E-2, 1.74E-2, 1.65E-2, 1.55E-2},
   { 7.50, 2.18E-2, 2.09E-2, 2.01E-2, 1.93E-2, 1.84E-2, 1.76E-2, 1.67E-2, 1.58E-2, 1.49E-2, 1.39E-2, 1.30E-2},
   { 8.00, 1.94E-2, 1.86E-2, 1.78E-2, 1.70E-2, 1.62E-2, 1.53E-2, 1.45E-2, 1.36E-2, 1.27E-2, 1.18E-2, 1.09E-2},
   { 8.50, 1.74E-2, 1.66E-2, 1.58E-2, 1.50E-2, 1.42E-2, 1.34E-2, 1.26E-2, 1.18E-2, 1.09E-2, 1.00E-2, 9.17E-3},
   { 9.00, 1.57E-2, 1.49E-2, 1.41E-2, 1.34E-2, 1.26E-2, 1.18E-2, 1.10E-2, 1.02E-2, 9.39E-3, 8.56E-3, 7.72E-3},
   { 9.50, 1.42E-2, 1.34E-2, 1.27E-2, 1.19E-2, 1.12E-2, 1.04E-2, 9.66E-3, 8.88E-3, 8.09E-3, 7.30E-3, 6.49E-3},
   {10.00, 1.28E-2, 1.21E-2, 1.14E-2, 1.07E-2, 9.99E-3, 9.25E-3, 8.51E-3, 7.75E-3, 6.99E-3, 6.23E-3, 5.45E-3},
   {11.00, 1.07E-2, 1.00E-2, 9.37E-3, 8.70E-3, 8.02E-3, 7.34E-3, 6.64E-3, 5.95E-3, 5.24E-3, 4.53E-3, 3.81E-3},
   {12.00, 9.01E-3, 8.39E-3, 7.77E-3, 7.14E-3, 6.50E-3, 5.87E-3, 5.22E-3, 4.58E-3, 3.93E-3, 3.27E-3, 2.61E-3},
   {13.00, 7.35E-3, 6.79E-3, 6.23E-3, 5.68E-3, 5.12E-3, 4.55E-3, 3.99E-3, 3.43E-3, 2.86E-3, 2.29E-3, 1.73E-3},
   {14.00, 5.54E-3, 5.08E-3, 4.63E-3, 4.18E-3, 3.73E-3, 3.28E-3, 2.84E-3, 2.40E-3, 1.97E-3, 1.54E-3, 1.11E-3},
   {15.00, 3.97E-3, 3.61E-3, 3.27E-3, 2.92E-3, 2.59E-3, 2.26E-3, 1.93E-3, 1.62E-3, 1.31E-3, 1.00E-3, 7.05E-4},
   {16.00, 2.81E-3, 2.54E-3, 2.28E-3, 2.02E-3, 1.78E-3, 1.54E-3, 1.30E-3, 1.08E-3, 8.59E-4, 6.49E-4, 4.48E-4},
   {17.00, 2.00E-3, 1.80E-3, 1.60E-3, 1.41E-3, 1.23E-3, 1.05E-3, 8.82E-4, 7.21E-4, 5.68E-4, 4.23E-4, 2.86E-4},
   {18.00, 1.44E-3, 1.29E-3, 1.14E-3, 9.93E-4, 8.56E-4, 7.26E-4, 6.03E-4, 4.87E-4, 3.79E-4, 2.77E-4, 1.83E-4},
   {19.00, 1.05E-3, 9.32E-4, 8.17E-4, 7.08E-4, 6.05E-4, 5.08E-4, 4.17E-4, 3.33E-4, 2.55E-4, 1.83E-4, 1.18E-4},
   {20.00, 7.77E-4, 6.83E-4, 5.94E-4, 5.10E-4, 4.32E-4, 3.59E-4, 2.91E-4, 2.29E-4, 1.72E-4, 1.21E-4, 7.58E-5},
   {22.00, 4.31E-4, 3.73E-4, 3.20E-4, 2.70E-4, 2.24E-4, 1.82E-4, 1.44E-4, 1.10E-4, 7.96E-5, 5.34E-5, 3.11E-5},
   {24.00, 2.40E-4, 2.05E-4, 1.73E-4, 1.43E-4, 1.16E-4, 9.22E-5, 7.08E-5, 5.32E-5, 3.63E-5, 2.30E-5, 1.23E-5},
   {26.00, 1.30E-4, 1.09E-4, 9.02E-5, 7.34E-5, 5.84E-5, 4.51E-5, 3.37E-5, 2.39E-5, 1.58E-5, 0,       0}};


//kappa = 0.10
static double vavilovPdfValues3[41][12]={
   {-3.50, 0,       0,       0,       0,       1.00E-5, 1.06E-5, 1.11E-5, 1.18E-5, 1.24E-5, 1.31E-5, 1.38E-5},
   {-3.00, 7.45E-4, 7.82E-4, 8.21E-4, 8.62E-4, 9.05E-4, 9.50E-4, 9.98E-4, 1.05E-3, 1.10E-3, 1.15E-3, 1.21E-3},
   {-2.50, 1.07E-2, 1.11E-2, 1.16E-2, 1.21E-2, 1.27E-2, 1.33E-2, 1.38E-2, 1.45E-2, 1.51E-2, 1.58E-2, 1.65E-2},
   {-2.00, 4.86E-2, 5.05E-2, 5.25E-2, 5.46E-2, 5.67E-2, 5.90E-2, 6.13E-2, 6.37E-2, 6.62E-2, 6.88E-2, 7.15E-2},
   {-1.50, 1.11E-1, 1.15E-1, 1.19E-1, 1.23E-1, 1.27E-1, 1.32E-1, 1.36E-1, 1.41E-1, 1.45E-1, 1.50E-1, 1.55E-1},
   {-1.00, 1.67E-1, 1.72E-1, 1.77E-1, 1.82E-1, 1.88E-1, 1.93E-1, 1.99E-1, 2.04E-1, 2.10E-1, 2.16E-1, 2.22E-1},
   {-0.50, 1.96E-1, 2.01E-1, 2.06E-1, 2.10E-1, 2.15E-1, 2.20E-1, 2.26E-1, 2.31E-1, 2.36E-1, 2.42E-1, 2.47E-1},
   { 0.00, 1.98E-1, 2.01E-1, 2.05E-1, 2.09E-1, 2.13E-1, 2.17E-1, 2.21E-1, 2.25E-1, 2.28E-1, 2.32E-1, 2.37E-1},
   { 0.50, 1.83E-1, 1.85E-1, 1.88E-1, 1.90E-1, 1.93E-1, 1.95E-1, 1.98E-1, 2.00E-1, 2.02E-1, 2.05E-1, 2.07E-1},
   { 1.00, 1.60E-1, 1.62E-1, 1.63E-1, 1.65E-1, 1.66E-1, 1.67E-1, 1.68E-1, 1.69E-1, 1.70E-1, 1.71E-1, 1.72E-1},
   { 1.50, 1.37E-1, 1.38E-1, 1.38E-1, 1.39E-1, 1.39E-1, 1.39E-1, 1.39E-1, 1.39E-1, 1.39E-1, 1.39E-1, 1.39E-1},
   { 2.00, 1.16E-1, 1.16E-1, 1.16E-1, 1.15E-1, 1.15E-1, 1.14E-1, 1.14E-1, 1.13E-1, 1.13E-1, 1.12E-1, 1.11E-1},
   { 2.50, 9.75E-2, 9.69E-2, 9.62E-2, 9.54E-2, 9.45E-2, 9.36E-2, 9.26E-2, 9.16E-2, 9.04E-2, 8.92E-2, 8.79E-2},
   { 3.00, 8.21E-2, 8.11E-2, 8.01E-2, 7.90E-2, 7.79E-2, 7.66E-2, 7.54E-2, 7.40E-2, 7.26E-2, 7.10E-2, 6.94E-2},
   { 3.50, 6.93E-2, 6.82E-2, 6.70E-2, 6.57E-2, 6.43E-2, 6.29E-2, 6.15E-2, 5.99E-2, 5.83E-2, 5.67E-2, 5.49E-2},
   { 4.00, 5.89E-2, 5.76E-2, 5.63E-2, 5.49E-2, 5.34E-2, 5.19E-2, 5.04E-2, 4.88E-2, 4.71E-2, 4.53E-2, 4.35E-2},
   { 4.50, 5.03E-2, 4.90E-2, 4.76E-2, 4.61E-2, 4.46E-2, 4.31E-2, 4.15E-2, 3.99E-2, 3.82E-2, 3.64E-2, 3.46E-2},
   { 5.00, 4.33E-2, 4.19E-2, 4.05E-2, 3.90E-2, 3.75E-2, 3.60E-2, 3.44E-2, 3.27E-2, 3.11E-2, 2.93E-2, 2.75E-2},
   { 5.50, 3.75E-2, 3.61E-2, 3.47E-2, 3.32E-2, 3.17E-2, 3.02E-2, 2.86E-2, 2.70E-2, 2.54E-2, 2.37E-2, 2.20E-2},
   { 6.00, 3.27E-2, 3.13E-2, 2.99E-2, 2.85E-2, 2.70E-2, 2.55E-2, 2.40E-2, 2.24E-2, 2.08E-2, 1.92E-2, 1.75E-2},
   { 6.50, 2.86E-2, 2.73E-2, 2.59E-2, 2.45E-2, 2.31E-2, 2.17E-2, 2.02E-2, 1.87E-2, 1.71E-2, 1.55E-2, 1.39E-2},
   { 7.00, 2.53E-2, 2.40E-2, 2.26E-2, 2.13E-2, 1.99E-2, 1.85E-2, 1.70E-2, 1.56E-2, 1.41E-2, 1.26E-2, 1.11E-2},
   { 7.50, 2.24E-2, 2.11E-2, 1.98E-2, 1.85E-2, 1.72E-2, 1.58E-2, 1.45E-2, 1.31E-2, 1.16E-2, 1.02E-2, 8.73E-3},
   { 8.00, 1.98E-2, 1.86E-2, 1.74E-2, 1.61E-2, 1.48E-2, 1.35E-2, 1.22E-2, 1.09E-2, 9.57E-3, 8.21E-3, 6.83E-3},
   { 8.50, 1.74E-2, 1.62E-2, 1.51E-2, 1.39E-2, 1.27E-2, 1.15E-2, 1.03E-2, 9.04E-3, 7.80E-3, 6.55E-3, 5.29E-3},
   { 9.00, 1.50E-2, 1.39E-2, 1.28E-2, 1.18E-2, 1.07E-2, 9.57E-3, 8.47E-3, 7.37E-3, 6.27E-3, 5.16E-3, 4.06E-3},
   { 9.50, 1.27E-2, 1.17E-2, 1.07E-2, 9.77E-3, 8.80E-3, 7.84E-3, 6.88E-3, 5.92E-3, 4.97E-3, 4.03E-3, 3.09E-3},
   {10.00, 1.05E-2, 9.69E-3, 8.84E-3, 7.99E-3, 7.16E-3, 6.33E-3, 5.51E-3, 4.70E-3, 3.90E-3, 3.11E-3, 2.33E-3},
   {11.00, 7.12E-3, 6.48E-3, 5.85E-3, 5.23E-3, 4.62E-3, 4.03E-3, 3.45E-3, 2.89E-3, 2.35E-3, 1.83E-3, 1.32E-3},
   {12.00, 4.77E-3, 4.30E-3, 3.84E-3, 3.39E-3, 2.96E-3, 2.54E-3, 2.15E-3, 1.77E-3, 1.41E-3, 1.07E-3, 7.45E-4},
   {13.00, 3.21E-3, 2.86E-3, 2.53E-3, 2.21E-3, 1.90E-3, 1.61E-3, 1.34E-3, 1.08E-3, 8.42E-4, 6.21E-4, 4.18E-4},
   {14.00, 2.18E-3, 1.92E-3, 1.68E-3, 1.45E-3, 1.23E-3, 1.03E-3, 8.38E-4, 6.65E-4, 5.06E-4, 3.62E-4, 2.34E-4},
   {15.00, 1.49E-3, 1.30E-3, 1.12E-3, 9.53E-4, 7.99E-4, 6.57E-4, 5.27E-4, 4.09E-4, 3.04E-4, 2.11E-4, 1.30E-4},
   {16.00, 1.01E-3, 8.76E-4, 7.47E-4, 6.28E-4, 5.19E-4, 4.20E-4, 3.31E-4, 2.51E-4, 1.81E-4, 1.21E-4, 7.12E-5},
   {17.00, 6.88E-4, 5.88E-4, 4.96E-4, 4.12E-4, 3.35E-4, 2.67E-4, 2.06E-4, 1.53E-4, 1.07E-4, 6.91E-5, 3.84E-5},
   {18.00, 4.60E-4, 3.89E-4, 3.24E-4, 2.66E-4, 2.13E-4, 1.67E-4, 1.26E-4, 9.16E-5, 6.25E-5, 3.87E-5, 2.03E-5},
   {19.00, 3.01E-4, 2.52E-4, 2.07E-4, 1.68E-4, 1.33E-4, 1.02E-4, 7.62E-5, 5.41E-5, 3.59E-5, 2.15E-5, 1.07E-5},
   {20.00, 1.94E-4, 1.61E-4, 1.31E-4, 1.05E-4, 8.18E-5, 6.21E-5, 4.55E-5, 3.17E-5, 2.06E-5, 1.20E-5, 0},
   {22.00, 7.97E-5, 6.47E-5, 5.16E-5, 4.03E-5, 3.07E-5, 2.27E-5, 1.61E-5, 1.08E-5, 0,       0,       0},
   {24.00, 3.26E-5, 2.59E-5, 2.01E-5, 1.53E-5, 1.12E-5, 0,       0,       0,       0,       0,       0},
   {26.00, 1.33E-5, 1.04E-5, 0,       0,       0,       0,       0,       0,       0,       0,       0}}; 


// kappa = 0.40
static double vavilovPdfValues4[28][12]={
   {-3.50, 1.07E-5, 1.26E-5, 1.47E-5, 1.73E-5, 2.03E-5, 2.37E-5, 2.78E-5, 3.26E-5, 3.82E-5, 4.48E-5, 5.24E-5},
   {-3.00, 1.00E-3, 1.16E-3, 1.33E-3, 1.53E-3, 1.75E-3, 2.02E-3, 2.32E-3, 2.66E-3, 3.05E-3, 3.50E-3, 4.02E-3},
   {-2.50, 1.44E-2, 1.62E-2, 1.83E-2, 2.06E-2, 2.32E-2, 2.61E-2, 2.93E-2, 3.30E-2, 3.71E-2, 4.17E-2, 4.68E-2},
   {-2.00, 6.56E-2, 7.25E-2, 8.00E-2, 8.83E-2, 9.74E-2, 1.07E-1, 1.18E-1, 1.30E-1, 1.43E-1, 1.57E-1, 1.73E-1},
   {-1.50, 1.50E-1, 1.62E-1, 1.76E-1, 1.90E-1, 2.05E-1, 2.21E-1, 2.38E-1, 2.57E-1, 2.76E-1, 2.97E-1, 3.19E-1},
   {-1.00, 2.26E-1, 2.40E-1, 2.54E-1, 2.69E-1, 2.84E-1, 3.00E-1, 3.16E-1, 3.32E-1, 3.49E-1, 3.66E-1, 3.83E-1},
   {-0.50, 2.65E-1, 2.75E-1, 2.85E-1, 2.96E-1, 3.05E-1, 3.15E-1, 3.24E-1, 3.32E-1, 3.40E-1, 3.47E-1, 3.52E-1},
   { 0.00, 2.66E-1, 2.71E-1, 2.76E-1, 2.79E-1, 2.82E-1, 2.84E-1, 2.84E-1, 2.84E-1, 2.82E-1, 2.78E-1, 2.73E-1},
   { 0.50, 2.44E-1, 2.43E-1, 2.42E-1, 2.39E-1, 2.36E-1, 2.31E-1, 2.25E-1, 2.18E-1, 2.09E-1, 1.99E-1, 1.88E-1},
   { 1.00, 2.07E-1, 2.03E-1, 1.97E-1, 1.91E-1, 1.83E-1, 1.75E-1, 1.65E-1, 1.55E-1, 1.44E-1, 1.31E-1, 1.18E-1},
   { 1.50, 1.66E-1, 1.59E-1, 1.51E-1, 1.43E-1, 1.34E-1, 1.24E-1, 1.14E-1, 1.03E-1, 9.21E-2, 8.05E-2, 6.86E-2},
   { 2.00, 1.26E-1, 1.18E-1, 1.10E-1, 1.01E-1, 9.25E-2, 8.35E-2, 7.43E-2, 6.51E-2, 5.58E-2, 4.66E-2, 3.75E-2},
   { 2.50, 9.11E-2, 8.37E-2, 7.62E-2, 6.87E-2, 6.11E-2, 5.36E-2, 4.63E-2, 3.91E-2, 3.22E-2, 2.56E-2, 1.94E-2},
   { 3.00, 6.35E-2, 5.71E-2, 5.09E-2, 4.48E-2, 3.88E-2, 3.31E-2, 2.77E-2, 2.26E-2, 1.78E-2, 1.35E-2, 9.60E-3},
   { 3.50, 4.28E-2, 3.77E-2, 3.29E-2, 2.82E-2, 2.39E-2, 1.98E-2, 1.60E-2, 1.26E-2, 9.52E-3, 6.84E-3, 4.55E-3},
   { 4.00, 2.79E-2, 2.41E-2, 2.06E-2, 1.72E-2, 1.42E-2, 1.14E-2, 8.95E-3, 6.78E-3, 4.91E-3, 3.35E-3, 2.08E-3},
   { 4.50, 1.77E-2, 1.50E-2, 1.25E-2, 1.02E-2, 8.21E-3, 6.42E-3, 4.87E-3, 3.55E-3, 2.46E-3, 1.59E-3, 9.22E-4},
   { 5.00, 1.10E-2, 9.10E-3, 7.42E-3, 5.93E-3, 4.63E-3, 3.51E-3, 2.58E-3, 1.81E-3, 1.20E-3, 7.34E-4, 3.96E-4},
   { 5.50, 6.65E-3, 5.39E-3, 4.30E-3, 3.35E-3, 2.55E-3, 1.88E-3, 1.33E-3, 9.02E-4, 5.72E-4, 3.30E-4, 1.66E-4},
   { 6.00, 3.93E-3, 3.13E-3, 2.44E-3, 1.85E-3, 1.37E-3, 9.83E-4, 6.75E-4, 4.39E-4, 2.66E-4, 1.45E-4, 6.73E-5},
   { 6.50, 2.28E-3, 1.78E-3, 1.35E-3, 1.01E-3, 7.25E-4, 5.04E-4, 3.34E-4, 2.09E-4, 1.21E-4, 6.24E-5, 2.68E-5},
   { 7.00, 1.30E-3, 9.91E-4, 7.38E-4, 5.35E-4, 3.76E-4, 2.53E-4, 1.63E-4, 9.79E-5, 5.41E-5, 2.64E-5, 1.05E-5},
   { 7.50, 7.27E-4, 5.43E-4, 3.96E-4, 2.80E-4, 1.91E-4, 1.25E-4, 7.76E-5, 4.49E-5, 2.36E-5, 1.08E-5, 0},
   { 8.00, 4.00E-4, 2.92E-4, 2.08E-4, 1.44E-4, 9.57E-5, 6.08E-5, 3.64E-5, 2.02E-5, 1.02E-5, 0,       0},
   { 8.50, 2.17E-4, 1.55E-4, 1.08E-4, 7.29E-5, 4.72E-5, 2.91E-5, 1.69E-5, 0,       0,       0,       0},
   { 9.00, 1.16E-4, 8.12E-5, 5.53E-5, 3.63E-5, 2.29E-5, 1.37E-5, 0,       0,       0,       0,       0},
   { 9.50, 6.08E-5, 4.19E-5, 2.79E-5, 1.79E-5, 1.09E-5, 0,       0,       0,       0,       0,       0},
   {10.00, 3.17E-5, 2.13E-5, 1.39E-5, 0,       0,       0,       0,       0,       0,       0,       0}};


//kappa = 0.70
static double vavilovPdfValues5[22][12]={
   {-3.50, 1.44E-5, 1.83E-5, 2.32E-5, 2.95E-5, 3.75E-5, 4.76E-5, 6.04E-5, 7.69E-5, 9.72E-5, 1.23E-4, 1.56E-4},
   {-3.00, 1.36E-3, 1.67E-3, 2.04E-3, 2.50E-3, 3.07E-3, 3.76E-3, 4.60E-3, 5.62E-3, 6.87E-3, 8.39E-3, 1.02E-2},
   {-2.50, 1.94E-2, 2.30E-2, 2.72E-2, 3.22E-2, 3.81E-2, 4.49E-2, 5.29E-2, 6.23E-2, 7.33E-2, 8.61E-2, 1.01E-1},
   {-2.00, 8.86E-2, 1.01E-1, 1.16E-1, 1.32E-1, 1.50E-1, 1.71E-1, 1.94E-1, 2.19E-1, 2.47E-1, 2.79E-1, 3.13E-1},
   {-1.50, 2.02E-1, 2.23E-1, 2.46E-1, 2.70E-1, 2.96E-1, 3.23E-1, 3.52E-1, 3.82E-1, 4.13E-1, 4.44E-1, 4.76E-1},
   {-1.00, 3.03E-1, 3.23E-1, 3.42E-1, 3.62E-1, 3.81E-1, 3.99E-1, 4.15E-1, 4.30E-1, 4.42E-1, 4.52E-1, 4.57E-1},
   {-0.50, 3.44E-1, 3.54E-1, 3.62E-1, 3.68E-1, 3.71E-1, 3.72E-1, 3.70E-1, 3.64E-1, 3.54E-1, 3.40E-1, 3.22E-1},
   { 0.00, 3.23E-1, 3.20E-1, 3.15E-1, 3.09E-1, 2.97E-1, 2.85E-1, 2.69E-1, 2.51E-1, 2.30E-1, 2.07E-1, 1.81E-1},
   { 0.50, 2.60E-1, 2.49E-1, 2.36E-1, 2.21E-1, 2.05E-1, 1.87E-1, 1.68E-1, 1.48E-1, 1.27E-1, 1.06E-1, 8.54E-2},
   { 1.00, 1.86E-1, 1.72E-1, 1.57E-1, 1.41E-1, 1.25E-1, 1.09E-1, 9.28E-2, 7.71E-2, 6.20E-2, 4.79E-2, 3.50E-2},
   { 1.50, 1.21E-1, 1.08E-1, 9.44E-2, 8.15E-2, 6.91E-2, 5.73E-2, 4.63E-2, 3.62E-2, 2.72E-2, 1.93E-2, 1.28E-2},
   { 2.00, 7.20E-2, 6.19E-2, 5.22E-2, 4.33E-2, 3.51E-2, 2.79E-2, 2.11E-2, 1.55E-2, 1.09E-2, 7.11E-3, 4.22E-3},
   { 2.50, 3.99E-2, 3.31E-2, 2.69E-2, 2.13E-2, 1.65E-2, 1.24E-2, 8.97E-3, 6.19E-3, 4.02E-3, 2.41E-3, 1.28E-3},
   { 3.00, 2.07E-2, 1.66E-2, 1.30E-2, 9.87E-3, 7.30E-3, 5.21E-3, 3.56E-3, 2.31E-3, 1.39E-3, 7.61E-4, 3.60E-4},
   { 3.50, 1.02E-2, 7.84E-3, 5.89E-3, 4.31E-3, 3.04E-3, 2.06E-3, 1.33E-3, 8.09E-4, 4.52E-4, 2.26E-4, 9.47E-5},
   { 4.00, 4.74E-3, 3.52E-3, 2.55E-3, 1.78E-3, 1.20E-3, 7.76E-4, 4.73E-4, 2.69E-4, 1.39E-4, 6.32E-5, 2.34E-5},
   { 4.50, 2.10E-3, 1.51E-3, 1.05E-3, 7.05E-4, 4.54E-4, 2.78E-4, 1.60E-4, 8.52E-5, 4.08E-5, 1.68E-5, 0},       
   { 5.00, 8.98E-4, 6.21E-4, 4.15E-4, 2.67E-4, 1.64E-4, 9.55E-5, 5.19E-5, 2.58E-5, 1.14E-5, 0,       0},       
   { 5.50, 3.68E-4, 2.45E-4, 1.58E-4, 9.71E-5, 5.70E-5, 3.15E-5, 1.61E-5, 0,       0,       0,       0},       
   { 6.00, 1.45E-4, 9.32E-5, 5.76E-5, 3.41E-5, 1.91E-5, 1.00E-5, 0,       0,       0,       0,       0},       
   { 6.50, 5.53E-5, 3.43E-5, 2.04E-5, 1.15E-5, 0,       0,       0,       0,       0,       0,       0},       
   { 7.00, 2.04E-5, 1.22E-5, 0,       0,       0,       0,       0,       0,       0,       0,       0}};


// kappa = 1.00
static double vavilovPdfValues6[19][12]={
   {-3.50, 1.94E-5, 2.64E-5, 3.59E-5, 4.87E-5, 6.61E-5, 8.96E-5, 1.21E-4, 1.64E-4, 2.22E-4, 3.00E-4, 4.05E-4},       
   {-3.00, 1.83E-3, 2.37E-3, 3.06E-3, 3.94E-3, 5.08E-3, 6.53E-3, 8.39E-3, 1.08E-2, 1.38E-2, 1.76E-2, 2.25E-2},       
   {-2.50, 2.62E-2, 3.22E-2, 3.95E-2, 4.84E-2, 5.91E-2, 7.20E-2, 8.76E-2, 1.06E-1, 1.28E-1, 1.55E-1, 1.86E-1},       
   {-2.00, 1.19E-1, 1.39E-1, 1.63E-1, 1.89E-1, 2.18E-1, 2.51E-1, 2.88E-1, 3.29E-1, 3.74E-1, 4.23E-1, 4.75E-1},       
   {-1.50, 2.69E-1, 2.99E-1, 3.31E-1, 3.64E-1, 3.97E-1, 4.31E-1, 4.65E-1, 4.97E-1, 5.26E-1, 5.52E-1, 5.73E-1},       
   {-1.00, 3.85E-1, 4.07E-1, 4.26E-1, 4.43E-1, 4.56E-1, 4.66E-1, 4.69E-1, 4.67E-1, 4.58E-1, 4.42E-1, 4.17E-1},       
   {-0.50, 4.01E-1, 4.02E-1, 4.00E-1, 3.92E-1, 3.80E-1, 3.63E-1, 3.42E-1, 3.15E-1, 2.84E-1, 2.49E-1, 2.11E-1},       
   { 0.00, 3.29E-1, 3.13E-1, 2.94E-1, 2.73E-1, 2.49E-1, 2.22E-1, 1.94E-1, 1.65E-1, 1.36E-1, 1.08E-1, 8.10E-2},       
   { 0.50, 2.23E-1, 2.02E-1, 1.80E-1, 1.57E-1, 1.34E-1, 1.12E-1, 9.10E-2, 7.13E-2, 5.35E-2, 3.80E-2, 2.50E-2},       
   { 1.00, 1.29E-1, 1.11E-1, 9.38E-2, 7.73E-2, 6.21E-2, 4.84E-2, 3.64E-2, 2.61E-2, 1.78E-2, 1.12E-2, 6.42E-3},       
   { 1.50, 6.60E-2, 5.39E-2, 4.30E-2, 3.34E-2, 2.51E-2, 1.83E-2, 1.27E-2, 8.37E-3, 5.15E-3, 2.88E-3, 1.42E-3},       
   { 2.00, 3.01E-2, 2.33E-2, 1.76E-2, 1.29E-2, 9.10E-3, 6.15E-3, 3.95E-3, 2.38E-3, 1.32E-3, 6.54E-4, 2.74E-4},       
   { 2.50, 1.24E-2, 9.15E-3, 6.53E-3, 4.50E-3, 2.98E-3, 1.88E-3, 1.11E-3, 6.12E-4, 3.05E-4, 1.33E-4, 4.73E-5},       
   { 3.00, 4.71E-3, 3.29E-3, 2.22E-3, 1.44E-3, 8.94E-4, 5.24E-4, 2.87E-4, 1.44E-4, 6.44E-5, 2.46E-5, 0},             
   { 3.50, 1.65E-3, 1.09E-3, 6.99E-4, 4.28E-4, 2.48E-4, 1.35E-4, 6.81E-5, 3.11E-5, 1.55E-5, 0,       0},             
   { 4.00, 5.38E-4, 3.39E-4, 2.05E-4, 1.18E-4, 6.41E-5, 3.24E-5, 1.51E-5, 0,       0,       0,       0},             
   { 4.50, 1.65E-4, 9.86E-5, 5.64E-5, 3.05E-5, 1.55E-5, 0,       0,       0,       0,       0,       0},             
   { 5.00, 4.75E-5, 2.70E-5, 1.46E-5, 0,       0,       0,       0,       0,       0,       0,       0},             
   { 5.50, 1.30E-5, 0,       0,       0,       0,       0,       0,       0,       0,       0,       0}};     


// kappa = 4.00
static double vavilovPdfValues7[20][12]={
   {-4.00, 0,       0,       0,       0,       0,       1.32E-5, 3.04E-5, 6.90E-5, 1.55E-4, 3.44E-4, 7.55E-4},       
   {-3.75, 1.38E-5, 2.98E-5, 6.38E-5, 1.35E-4, 2.83E-4, 5.85E-4, 1.20E-3, 2.41E-3, 4.78E-3, 9.35E-3, 1.80E-2},       
   {-3.50, 3.81E-4, 7.44E-4, 1.44E-3, 2.73E-3, 5.12E-3, 9.45E-3, 1.71E-2, 3.05E-2, 5.33E-2, 9.11E-2, 1.52E-1},       
   {-3.25, 4.78E-3, 8.44E-3, 1.47E-2, 2.50E-2, 4.19E-2, 6.88E-2, 1.10E-1, 1.73E-1, 2.64E-1, 3.90E-1, 5.57E-1},       
   {-3.00, 3.19E-2, 5.08E-2, 7.95E-2, 1.22E-1, 1.82E-1, 2.64E-1, 3.73E-1, 5.11E-1, 6.75E-1, 8.56E-1, 1.03E+0},       
   {-2.75, 1.27E-1, 1.83E-1, 2.56E-1, 3.51E-1, 4.66E-1, 6.00E-1, 7.44E-1, 8.86E-1, 1.01E+0, 1.08E+0, 1.09E+0},       
   {-2.50, 3.28E-1, 4.26E-1, 5.38E-1, 6.58E-1, 7.77E-1, 8.81E-1, 9.56E-1, 9.85E-1, 9.56E-1, 8.63E-1, 7.14E-1},       
   {-2.25, 5.89E-1, 6.91E-1, 7.84E-1, 8.57E-1, 8.97E-1, 8.95E-1, 8.46E-1, 7.51E-1, 6.18E-1, 4.65E-1, 3.12E-1},       
   {-2.00, 7.75E-1, 8.22E-1, 8.37E-1, 8.15E-1, 7.56E-1, 6.63E-1, 5.44E-1, 4.14E-1, 2.88E-1, 1.78E-1, 9.59E-2},       
   {-1.75, 7.79E-1, 7.45E-1, 6.81E-1, 5.91E-1, 4.85E-1, 3.73E-1, 2.65E-1, 1.72E-1, 1.01E-1, 5.10E-2, 2.17E-2},       
   {-1.50, 6.16E-1, 5.32E-1, 4.36E-1, 3.38E-1, 2.45E-1, 1.64E-1, 1.01E-1, 5.60E-2, 2.73E-2, 1.13E-2, 3.74E-3},       
   {-1.25, 3.95E-1, 3.07E-1, 2.26E-1, 1.56E-1, 9.96E-2, 5.85E-2, 3.11E-2, 1.46E-2, 5.90E-3, 1.97E-3, 5.05E-4},       
   {-1.00, 2.09E-1, 1.47E-1, 9.68E-2, 5.94E-2, 3.35E-2, 1.72E-2, 7.85E-3, 3.12E-3, 1.05E-3, 2.80E-4, 5.49E-5},       
   {-0.75, 9.33E-2, 5.91E-2, 3.49E-2, 1.91E-2, 9.47E-3, 4.23E-3, 1.66E-3, 5.59E-4, 1.54E-4, 3.29E-5, 0},       
   {-0.50, 3.56E-2, 2.04E-2, 1.08E-2, 5.22E-3, 2.29E-3, 8.90E-4, 3.00E-4, 8.49E-5, 1.93E-5, 0,       0},       
   {-0.25, 1.18E-2, 6.07E-3, 2.88E-3, 1.24E-3, 4.78E-4, 1.62E-4, 4.67E-5, 1.12E-5, 0,       0,       0},       
   { 0.00, 3.41E-3, 1.59E-3, 6.74E-4, 2.58E-4, 8.75E-5, 2.57E-5, 0,       0,       0,       0,       0},       
   { 0.25, 8.74E-4, 3.67E-4, 1.40E-4, 4.74E-5, 1.42E-5, 0,       0,       0,       0,       0,       0},       
   { 0.50, 2.00E-4, 7.57E-5, 2.58E-5, 0,       0,       0,       0,       0,       0,       0,       0},       
   { 0.75, 4.11E-5, 1.41E-5, 0,       0,       0,       0,       0,       0,       0,       0,       0}};


// kappa = 7.00
static double vavilovPdfValues8[21][12]={
   {-4.40, 0,       0,       0,       0,       0,       0,       0,       0,       1.11E-5, 3.90E-5, 1.33E-4},       
   {-4.20, 0,       0,       0,       0,       0,       2.17E-5, 6.93E-5, 2.15E-4, 6.46E-4, 1.88E-3, 5.28E-3},       
   {-4.00, 0,       1.07E-5, 3.24E-5, 9.54E-5, 2.73E-4, 7.60E-4, 2.04E-3, 5.31E-3, 1.33E-2, 3.18E-2, 7.28E-2},       
   {-3.80, 1.11E-4, 2.99E-4, 7.80E-4, 1.97E-3, 4.82E-3, 1.14E-2, 2.57E-2, 5.57E-2, 1.15E-1, 2.24E-1, 4.12E-1},       
   {-3.60, 1.77E-3, 4.11E-3, 9.26E-3, 2.01E-2, 4.18E-2, 8.32E-2, 1.58E-1, 2.83E-1, 4.78E-1, 7.51E-1, 1.09E+0},       
   {-3.40, 1.54E-2, 3.10E-2, 6.01E-2, 1.12E-1, 1.97E-1, 3.31E-1, 5.23E-1, 7.74E-1, 1.06E+0, 1.33E+0, 1.50E+0},       
   {-3.20, 7.96E-2, 1.39E-1, 2.33E-1, 3.69E-1, 5.54E-1, 7.81E-1, 1.02E+0, 1.24E+0, 1.37E+0, 1.35E+0, 1.17E+0},       
   {-3.00, 2.63E-1, 3.99E-1, 5.74E-1, 7.78E-1, 9.89E-1, 1.17E+0, 1.27E+0, 1.25E+0, 1.10E+0, 8.46E-1, 5.51E-1},       
   {-2.80, 5.86E-1, 7.71E-1, 9.54E-1, 1.10E+0, 1.18E+0, 1.17E+0, 1.04E+0, 8.35E-1, 5.83E-1, 3.46E-1, 1.68E-1},       
   {-2.60, 9.21E-1, 1.05E+0, 1.12E+0, 1.10E+0, 9.97E-1, 8.19E-1, 6.02E-1, 3.88E-1, 2.14E-1, 9.71E-2, 3.44E-2},       
   {-2.40, 1.06E+0, 1.04E+0, 9.55E-1, 8.02E-1, 6.12E-1, 4.18E-1, 2.52E-1, 1.30E-1, 5.63E-2, 1.94E-2, 4.95E-3},       
   {-2.20, 9.18E-1, 7.85E-1, 6.16E-1, 4.40E-1, 2.83E-1, 1.61E-1, 7.89E-2, 3.27E-2, 1.10E-2, 2.84E-3, 5.18E-4},       
   {-2.00, 6.17E-1, 4.57E-1, 3.08E-1, 1.87E-1, 1.01E-1, 4.75E-2, 1.90E-2, 6.29E-3, 1.64E-3, 3.16E-4, 4.05E-5},       
   {-1.80, 3.28E-1, 2.10E-1, 1.22E-1, 6.30E-2, 2.85E-2, 1.11E-2, 3.62E-3, 9.50E-4, 1.91E-4, 2.78E-5, 0},       
   {-1.60, 1.41E-1, 7.84E-2, 3.89E-2, 1.71E-2, 6.49E-3, 2.09E-3, 5.52E-4, 1.15E-4, 1.77E-5, 0,       0},       
   {-1.40, 4.99E-2, 2.40E-2, 1.02E-2, 3.80E-3, 1.21E-3, 3.22E-4, 6.88E-5, 1.13E-5, 0,       0,       0},       
   {-1.20, 1.47E-2, 6.10E-3, 2.23E-3, 7.05E-4, 1.88E-4, 4.12E-5, 0,       0,       0,       0,       0},       
   {-1.00, 3.64E-3, 1.31E-3, 4.11E-4, 1.10E-4, 2.46E-5, 0,       0,       0,       0,       0,       0},       
   {-0.80, 7.71E-4, 2.40E-4, 6.46E-5, 1.47E-5, 0,       0,       0,       0,       0,       0,       0},       
   {-0.60, 1.41E-4, 3.80E-5, 0,       0,       0,       0,       0,       0,       0,       0,       0},       
   {-0.40, 2.24E-5, 0,       0,       0,       0,       0,       0,       0,       0,       0,       0}};


// kappa = 10.00
static double vavilovPdfValues9[21][12]={
   {-4.60, 0,       0,       0,       0,       0,       0,       0,       0,       0,       3.85E-5, 1.79E-4},       
   {-4.40, 0,       0,       0,       0,       0,       1.18E-5, 5.14E-5, 2.12E-4, 8.32E-4, 3.08E-3, 1.07E-2},       
   {-4.20, 0,       0,       1.41E-5, 5.57E-5, 2.10E-4, 7.50E-4, 2.54E-3, 8.11E-3, 2.42E-2, 6.69E-2, 1.70E-1},       
   {-4.00, 5.27E-5, 1.84E-4, 6.15E-4, 1.95E-3, 5.83E-3, 1.64E-2, 4.32E-2, 1.06E-1, 2.37E-1, 4.82E-1, 8.79E-1},       
   {-3.80, 1.43E-3, 4.07E-3, 1.10E-2, 2.78E-2, 6.61E-2, 1.46E-1, 2.96E-1, 5.49E-1, 9.16E-1, 1.35E+0, 1.73E+0},       
   {-3.60, 1.79E-2, 4.17E-2, 9.09E-2, 1.85E-1, 3.46E-1, 5.96E-1, 9.30E-1, 1.30E+0, 1.59E+0, 1.68E+0, 1.47E+0},       
   {-3.40, 1.16E-1, 2.20E-1, 3.88E-1, 6.29E-1, 9.31E-1, 1.24E+0, 1.48E+0, 1.54E+0, 1.38E+0, 1.02E+0, 5.99E-1},       
   {-3.20, 4.21E-1, 6.50E-1, 9.23E-1, 1.19E+0, 1.39E+0, 1.44E+0, 1.30E+0, 1.01E+0, 6.46E-1, 3.31E-1, 1.28E-1},       
   {-3.00, 9.12E-1, 1.15E+0, 1.31E+0, 1.35E+0, 1.24E+0, 9.87E-1, 6.74E-1, 3.84E-1, 1.76E-1, 6.16E-2, 1.53E-2},       
   {-2.80, 1.25E+0, 1.28E+0, 1.18E+0, 9.66E-1, 6.92E-1, 4.25E-1, 2.18E-1, 9.11E-2, 2.95E-2, 6.96E-3, 1.09E-3},       
   {-2.60, 1.13E+0, 9.45E-1, 7.01E-1, 4.56E-1, 2.55E-1, 1.20E-1, 4.64E-2, 1.41E-2, 3.19E-3, 5.02E-4, 4.87E-5},       
   {-2.40, 7.06E-1, 4.80E-1, 2.87E-1, 1.48E-1, 6.46E-2, 2.33E-2, 6.71E-3, 1.47E-3, 2.33E-4, 2.41E-5, 0},       
   {-2.20, 3.13E-1, 1.73E-1, 8.33E-2, 3.41E-2, 1.16E-2, 3.19E-3, 6.84E-4, 1.08E-4, 1.18E-5, 0,       0},       
   {-2.00, 1.02E-1, 4.59E-2, 1.77E-2, 5.74E-3, 1.52E-3, 3.19E-4, 5.06E-5, 0,       0,       0,       0},       
   {-1.80, 2.48E-2, 9.10E-3, 2.82E-3, 7.24E-4, 1.49E-4, 2.37E-5, 0,       0,       0,       0,       0},       
   {-1.60, 4.64E-3, 1.38E-3, 3.45E-4, 6.99E-5, 1.11E-5, 0,       0,       0,       0,       0,       0},       
   {-1.40, 6.77E-4, 1.64E-4, 3.29E-5, 0,       0,       0,       0,       0,       0,       0,       0},       
   {-1.20, 7.84E-5, 1.55E-5, 0,       0,       0,       0,       0,       0,       0,       0,       0}};

static double (*vavilovPdfValues[10])[12]={vavilovPdfValues0, vavilovPdfValues1, vavilovPdfValues2, vavilovPdfValues3, vavilovPdfValues4, vavilovPdfValues5, vavilovPdfValues6, vavilovPdfValues7, vavilovPdfValues8, vavilovPdfValues9};

static double vavilovKappaValues[10] = {.01, .04, .07, .1, .4, .7, 1, 4, 7, 10};

static int vavilovNLambda[10] = {45, 42, 41, 41, 28, 22, 19, 20, 21, 21};
   
int VavilovTest::GetSBNKappa () {
  return 10;
}

double VavilovTest::GetSBKappa (int ikappa) {
  if (ikappa < 0 || ikappa >= VavilovTest::GetSBNKappa()) return -1;
  return vavilovKappaValues[ikappa];
}

int VavilovTest::GetSBNBeta2 () {
  return 11;
}

double VavilovTest::GetSBBeta2 (int ibeta2) {
  if (ibeta2 < 0 || ibeta2 >= VavilovTest::GetSBNBeta2()) return -1;
  return 0.1*ibeta2;
}

int VavilovTest::GetSBNLambda (int ikappa) {
  if (ikappa < 0 || ikappa >= VavilovTest::GetSBNKappa()) return -1;
  return vavilovNLambda[ikappa];
}

double VavilovTest::GetSBLambda (int ikappa, int ilambda) {
  if (ikappa < 0 || ikappa >= VavilovTest::GetSBNKappa()) return 0;
  if (ilambda < 0 || ilambda >= VavilovTest::GetSBNLambda(ikappa)) return 0;
  return vavilovPdfValues[ikappa][ilambda][0];
}

double VavilovTest::GetSBVavilovPdf (int ikappa, int ibeta2, int ilambda) {
  if (ikappa < 0 || ikappa >= VavilovTest::GetSBNKappa()) return 0;
  if (ibeta2 < 0 || ibeta2 >= VavilovTest::GetSBNBeta2()) return 0;
  if (ilambda < 0 || ilambda >= VavilovTest::GetSBNLambda(ikappa)) return 0;
  return vavilovPdfValues[ikappa][ilambda][ibeta2+1];
}

static double myRound (double x, double y, double& xmantissa, int digits) {
  int exponent;
  if (y)      exponent = std::floor(std::log10(fabs(y)));
  else if (x) exponent = std::floor(std::log10(fabs(x)));
  else exponent = 0;
  double power    = std::pow (10.0, exponent);
  double mantissa = y/power;
  double dpower   = std::pow (10.0, digits-1);
  mantissa = roundl (mantissa*dpower)/dpower;
  if (mantissa >= 10) {
    mantissa *= 0.1;
    exponent += 1;
    power    = std::pow (10.0, exponent);
  }
  xmantissa = x/power;
  mantissa = roundl (xmantissa*dpower)/dpower;
  return mantissa*power;
}
double myRound (double x, double y, int digits) {
  double xmantissa;
  return myRound (x, y, xmantissa, digits);

}

static std::string format (double x, double y, int digits, int width) {
  int exponent;
  if (y)      exponent = std::floor(std::log10(fabs(y)));
  else if (x) exponent = std::floor(std::log10(fabs(x)));
  else exponent = 0;
  double power    = std::pow (10.0, exponent);
  double mantissa = y/power;
  double dpower   = std::pow (10.0, digits-1);
  mantissa = roundl (mantissa*dpower)/dpower;
  if (mantissa >= 10) {
    mantissa *= 0.1;
    exponent += 1;
  }
  mantissa = roundl (x/power*dpower)/dpower;
  
  std::stringstream out;
  out << std::setw(width-4) << std::fixed << std::setprecision(digits-1) << mantissa;
  out << "e" << std::showpos << std::setw(3)<< std::internal << std::setfill('0') << exponent;
  return out.str();

}

int VavilovTest::PdfTest (ROOT::Math::Vavilov& v, std::ostream& os) {
   double maxabsdiff, maxdiffmantissa, agreefraction, agreediffmantissa;
   GetPdfTestParams (v, maxabsdiff, maxdiffmantissa, agreefraction, agreediffmantissa);
   return VavilovTest::PdfTest (v, os, maxabsdiff, maxdiffmantissa, agreefraction, agreediffmantissa);
}   

int VavilovTest::PdfTest (ROOT::Math::Vavilov& v, std::ostream& os, 
                          double maxabsdiff, double maxdiffmantissa, 
                          double agreefraction, double agreediffmantissa) {
                          
   std::ios::fmtflags defaultflags = os.flags();
   int defaultwidth = os.width();
   int defaultprecision = os.precision();
   
   int nfail = 0;
  
   os << "\n\nTesting Pdf\n\n";
   
   for (int ikappa = 0; ikappa < GetSBNKappa(); ++ikappa) {
      double kappa = GetSBKappa (ikappa);
      os << "\n\nkappa = " << kappa << "\n\n";
      
      double maxreldiff = 0;
      double maxmandiff = 0;
      bool pass = true;
      int agree = 0;
      int disagree = 0;
         
      for (int ilambda = 0; ilambda < GetSBNLambda (ikappa); ++ilambda) {
         double lambda = GetSBLambda (ikappa, ilambda);
         os << std::setw(5) << std::fixed << std::setprecision(2) << lambda;
         double x = lambda;
         for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
            double beta2 = GetSBBeta2 (ibeta);
            v.SetKappaBeta2 (kappa, beta2);
            double pdf = v.Pdf(x);
            double val = GetSBVavilovPdf (ikappa, ibeta, ilambda);
            double diffmantissa;
            myRound (pdf-val, val, diffmantissa, 3);
            
            if (val > 0 && pdf > 0) {
              double absdiff = fabs(pdf-val);
              if (absdiff > maxabsdiff && fabs(diffmantissa) > maxdiffmantissa) {
                //os << "FAIL";
                pass = false;
              }
              if (fabs(diffmantissa) > 0.006) {
                if (absdiff > maxabsdiff)            maxabsdiff = absdiff;
                if (absdiff/val > maxreldiff)        maxreldiff = absdiff/val;
                if (fabs(diffmantissa) > maxmandiff) maxmandiff = fabs(diffmantissa);
              } 
              if (fabs(diffmantissa) > agreediffmantissa) ++disagree; else ++agree;
            }
            if (val == 0)
               os << "          ";
            else if (pdf == 0)   
               os << "  ---     ";
            else if (fabs(diffmantissa) > 0.006) 
               os << format (pdf-val, val, 3, 10);
            else 
               os << "  0       ";
         }
         os << "\n";
      }
      os.flags (defaultflags);
      if (agree < disagree*agreefraction) {
        pass = false;
        //os << "agreefraction test failed.\n";
      }  
      if (!pass) ++nfail;
      os << "Max abs diff: " << maxabsdiff << ", max rel diff: " << maxreldiff
           << ", max diff mantissa: " << std::fixed << std::setprecision(2) << maxmandiff
           << ", agree/disagree=" << agree << "/" << disagree
           << ", pass=" << pass << std::endl;
  }
  os << "\n\nNumber of failed tests: " << nfail << std::endl;
  
  os.flags (defaultflags);
  os.width(defaultwidth);
  os.precision(defaultprecision);
  
  return nfail;
}

static void moments (ROOT::Math::Vavilov& v, double& integral,
                     double& mean, double& variance,
                     double& skewness, double& kurtosis) {
  int nsteps = 10000;
  double t0 = v.GetLambdaMin();
  double t1 = v.GetLambdaMax();
  double t = t1 - t0;
  double dt = t/nsteps;
  
  double sum = 0;
  double sumx = 0;
  for (int i = 0; i < nsteps; ++i) {
    double x = (i+0.5)*dt + t0;
    double pdf = v.Pdf(x);
    sum += pdf;
    sumx += x*pdf;
  }              
  integral = sum*dt;   
  mean = sumx/sum;    
  double sumx2 = 0;
  double sumx3 = 0;
  double sumx4 = 0;
  for (int i = 0; i < nsteps; ++i) {
    double x = (i+0.5)*dt + t0;
    double pdf = v.Pdf(x);
    sumx2 += pow(x-mean, 2)*pdf;
    sumx3 += pow(x-mean, 3)*pdf;
    sumx4 += pow(x-mean, 4)*pdf;
  }              
  variance =  sumx2/sum;   
  skewness =  sumx3/sum*pow (variance, -1.5);
  kurtosis =  sumx4/sum/(variance*variance)-3;
}              

void VavilovTest::PrintPdfTable (ROOT::Math::Vavilov& v, std::ostream& os, int digits) {
   std::ios::fmtflags defaultflags = os.flags();
   int defaultwidth = os.width();
   int defaultprecision = os.precision();
     
   for (int ikappa = 0; ikappa < GetSBNKappa(); ++ikappa) {
      double kappa = GetSBKappa (ikappa);
      os << "\n\nkappa = " << kappa << "\n\n";
               
      for (int ilambda = 0; ilambda < GetSBNLambda (ikappa); ++ilambda) {
         double lambda = GetSBLambda (ikappa, ilambda);
         os << std::setw(5) << std::fixed << std::setprecision(2) << lambda;
         double x = lambda;
         for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
            double beta2 = GetSBBeta2 (ibeta);
            v.SetKappaBeta2 (kappa, beta2);
            double pdf = v.Pdf(x);
            if (pdf > 0) 
               os << std::setw(digits+7) << std::scientific << std::setprecision(digits-1) << pdf;
            else
               os << std::setw(digits+7) << " ";
         }
         os << "\n";
      }
      os << "Xmin:";
      for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
         double beta2 = GetSBBeta2 (ibeta);
         v.SetKappaBeta2 (kappa, beta2);
         os << std::setw(digits+7) << std::fixed << std::setprecision(digits-1) << v.GetLambdaMin();
      }
      os << "\n";
      os << "Xmax:";
      for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
         double beta2 = GetSBBeta2 (ibeta);
         v.SetKappaBeta2 (kappa, beta2);
         os << std::setw(digits+7) << std::fixed << std::setprecision(digits-1) << v.GetLambdaMax();
      }
      os << "\n";
      os << "int: ";
      double integral[11], calcmean[11], calcvariance[11],  calcskewness[11], calckurtosis[11];        
      for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
         double beta2 = GetSBBeta2 (ibeta);
         v.SetKappaBeta2 (kappa, beta2);
         moments (v, integral[ibeta], calcmean[ibeta], calcvariance[ibeta],  calcskewness[ibeta], calckurtosis[ibeta]);
         os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << integral[ibeta];
      }
      os << "\nmean:";
      for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
         os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << calcmean[ibeta];
      }
      os << "\n     ";
      for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
         double beta2 = GetSBBeta2 (ibeta);
         os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << v.Mean (kappa, beta2);
      }
      os << "\nvar: ";
      for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
         os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << calcvariance[ibeta];
      }
      os << "\n     ";
      for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
         double beta2 = GetSBBeta2 (ibeta);
         os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << v.Variance (kappa, beta2);
      }
      os << "\nskew:";
      for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
         os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << calcskewness[ibeta];
      }
      os << "\n     ";
      for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
         double beta2 = GetSBBeta2 (ibeta);
         os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << v.Skewness (kappa, beta2);
      }
      os << "\nkurt:";
      for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
         os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << calckurtosis[ibeta];
      }
      os << "\n     ";
      for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
         double beta2 = GetSBBeta2 (ibeta);
         os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << v.Kurtosis (kappa, beta2);
      }
      os << "\n";
  }
  
  os.flags (defaultflags);
  os.width(defaultwidth);
  os.precision(defaultprecision);
}

int VavilovTest::CdfTest (ROOT::Math::Vavilov& v, std::ostream& os) {
   double maxabsdiff, maxcdfdiff;
   GetCdfTestParams (v, maxabsdiff, maxcdfdiff);
   return VavilovTest::CdfTest (v, os, maxabsdiff, maxcdfdiff);
}
   
int VavilovTest::CdfTest (ROOT::Math::Vavilov& v, std::ostream& os, double maxabsdiff, double maxcdfdiff) {   
   std::ios::fmtflags defaultflags = os.flags();
   int defaultwidth = os.width();
   int defaultprecision = os.precision();
   
   int nfail = 0;
   
   os << "\n\nTesting Cdf and Cdf_c\n\n";
   
   for (int ikappa = 0; ikappa < GetSBNKappa(); ++ikappa) {
      double kappa = GetSBKappa (ikappa);
      os << "\n\nkappa = " << kappa << "\n\n";
      
      double absdiffmax = 0;
      double reldiffmax = 0;
      double cdfdiffmax = 0;
      bool pass = true;
         
      double cdf_calc[11];
      for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
         cdf_calc[ibeta] = 0;
      }  
      
      for (int ilambda = 0; ilambda < GetSBNLambda (ikappa); ++ilambda) {
         double lambda1 = GetSBLambda (ikappa, ilambda);
         os << std::setw(5) << std::fixed << std::setprecision(2) << lambda1;
         for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
            double beta2 = GetSBBeta2 (ibeta);
            v.SetKappaBeta2 (kappa, beta2);
            double lambda0 = (ilambda > 0) ?  GetSBLambda (ikappa, ilambda-1) : v.GetLambdaMin();
            if (lambda1 > lambda0) {
              int n = 100;
              double dlambda = (lambda1 - lambda0)/n;
              for (int i = 0; i < n; ++i) {
                double lambda = lambda0 + (i+0.5)*dlambda;
                cdf_calc[ibeta] += v.Pdf(lambda)*dlambda;
              }
            }
            
            double cdf = v.Cdf(lambda1);
            double cdf_c = v.Cdf_c(lambda1);
            double val = cdf_calc[ibeta];
            
            if (fabs(cdf-val) > absdiffmax) absdiffmax = fabs(cdf-val);
            if (fabs(cdf+cdf_c-1) > cdfdiffmax) cdfdiffmax = cdf+cdf_c-1;
            if (val > 0 && fabs(cdf/val-1) > reldiffmax) reldiffmax = fabs(cdf/val-1);
            
            if (val == 0)
               os << "          ";
            else 
               os << std::scientific << std::setw(10) << std::setprecision(2) << cdf-val;
         }
         os << "\n";
      }
      os.flags (defaultflags);
      if (absdiffmax > maxabsdiff) pass = false;
      if (cdfdiffmax > maxcdfdiff) pass = false;
      if (!pass) ++nfail;
      os << "Max abs diff: " << absdiffmax << ", max rel diff: " << reldiffmax
         << ", max diff cdf+cdf_c-1: " << cdfdiffmax
           << ", pass=" << pass << std::endl;
  }
  os << "\n\nNumber of failed tests: " << nfail << std::endl;
  
  os.flags (defaultflags);
  os.width(defaultwidth);
  os.precision(defaultprecision);
  
  return nfail;
}

int VavilovTest::QuantileTest (ROOT::Math::Vavilov& v, std::ostream& os) {
   double maxabsdiff;
   GetQuantileTestParams (v, maxabsdiff);
   return VavilovTest::QuantileTest (v, os, maxabsdiff);
}

int VavilovTest::QuantileTest (ROOT::Math::Vavilov& v, std::ostream& os, double maxabsdiff) {
   double qmin = 0;
   
   std::ios::fmtflags defaultflags = os.flags();
   int defaultwidth = os.width();
   int defaultprecision = os.precision();
   
   int nfail = 0;
      
   static const double qvalues[45] = {0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009,
                                      0.01,  0.02,  0.03,  0.04,  0.05,  0.06,  0.07,  0.08,  0.09,
                                      0.1,   0.2,   0.3,   0.4,   0.5,   0.6,   0.7,   0.8,   0.9,
                                      0.91,  0.92,  0.93,  0.94,  0.95,  0.96,  0.97,  0.98,  0.99,
                                      0.991, 0.992, 0.993, 0.994, 0.995, 0.996, 0.997, 0.998, 0.999};
   
   os << "\n\nTesting Quantile\n\n";
   
   for (int ikappa = 0; ikappa < GetSBNKappa(); ++ikappa) {
      double kappa = GetSBKappa (ikappa);
      os << "\n\nkappa = " << kappa << "\n\n";
      
      double absdiffmax = 0;
      bool pass = true;
      
      for (int iq = 0; iq < 45; ++iq) {
         double qval = qvalues[iq];
         os << std::setw(5) << std::fixed << std::setprecision(3) << qval;
         for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
            double beta2 = GetSBBeta2 (ibeta);
            v.SetKappaBeta2 (kappa, beta2);
            double lambda =  v.Quantile (qval);
            double cdfval = v.Cdf(lambda);
            if (qval > qmin && (1-qval) > qmin) {            
            if (fabs(cdfval-qval) > absdiffmax) absdiffmax = fabs(cdfval-qval);
              os << "  " << std::fixed << std::setw(7) << std::setprecision(4) << cdfval-qval << " ";
            }
            else {
               os << " (" << std::fixed << std::setw(7) << std::setprecision(4) << cdfval-qval << ")";
            }
         }
         os << "\n";
      }
      os.flags (defaultflags);
      if (absdiffmax > maxabsdiff) pass = false;
      if (!pass) ++nfail;
      os << "Max abs diff: " << absdiffmax 
           << ", pass=" << pass << std::endl;
   }
      
   os << "\n\nTesting Quantile_c\n\n";
   
   for (int ikappa = 0; ikappa < GetSBNKappa(); ++ikappa) {
      double kappa = GetSBKappa (ikappa);
      os << "\n\nkappa = " << kappa << "\n\n";
      
      double absdiffmax = 0;
      bool pass = true;
    
      for (int iq = 0; iq < 45; ++iq) {
         double qval = qvalues[iq];
         os << std::setw(5) << std::fixed << std::setprecision(3) << qval;
         for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
            double beta2 = GetSBBeta2 (ibeta);
            v.SetKappaBeta2 (kappa, beta2);
            double lambda_c =  v.Quantile_c (qval);
            double cdf_c_val = v.Cdf_c(lambda_c);
            if (qval > qmin && (1-qval) > qmin) {            
            if (fabs(cdf_c_val-qval) > absdiffmax) absdiffmax = fabs(cdf_c_val-qval);
            
            os << "  " << std::fixed << std::setw(7) << std::setprecision(4) << cdf_c_val-qval << " ";
            }
            else {
               os << " (" << std::fixed << std::setw(7) << std::setprecision(4) << cdf_c_val-qval << ")";
            }
         }
         os << "\n";
      }
      os.flags (defaultflags);
      if (absdiffmax > maxabsdiff) pass = false;
      if (!pass) ++nfail;
      os << "Max abs diff: " << absdiffmax 
           << ", pass=" << pass << std::endl;
   }
   os << "\n\nNumber of failed tests: " << nfail << std::endl;
  
   os.flags (defaultflags);
   os.width(defaultwidth);
   os.precision(defaultprecision);
  
   return nfail;
}

void VavilovTest::GetPdfTestParams (const Vavilov& v, double& maxabsdiff, double& maxdiffmantissa, double& agreefraction, double& agreediffmantissa) {
   if (dynamic_cast <const VavilovFast *>(&v)) {
      maxabsdiff = 0.08;
      maxdiffmantissa = 0.1;
      agreefraction = 1;
      agreediffmantissa = 0.9;
   }
   else {
      maxabsdiff = 2E-3;
      maxdiffmantissa = 0.03;
      agreefraction = 5;
      agreediffmantissa = 0.015;
   }
}

void VavilovTest::GetCdfTestParams (const Vavilov& v, double& maxabsdiff, double& maxcdfdiff) {
   if (dynamic_cast <const VavilovFast *>(&v)) {
      maxabsdiff = 0.018;
      maxcdfdiff = 1E-16;
   }
   else {
      maxabsdiff = 1E-5;
      maxcdfdiff = 5E-15;
   }
}

void VavilovTest::GetQuantileTestParams (const Vavilov& v, double& maxabsdiff) {
   if (dynamic_cast <const VavilovFast *>(&v)) {
      maxabsdiff = 0.03;
   }
   else {
      maxabsdiff = 2E-10;
   }
}


} // namespace Math
} // namespace ROOT