/********************************************************************/
/* Copyright (c) 2017 System fugen G.K. and Yuzi Mizuno          */
/* All rights reserved.                                             */
/********************************************************************/
#ifndef _MGGausp_HH_
#define _MGGausp_HH_
/** @addtogroup ALGORITHM
 *  @{
 */

/** 
 *	@brief Legendre-Gauss quadratuer formula over (a,b) .
 * The DE formula (double exponential formula) is applied.
 * *Note* Original Fortran program codes are from DDGAUSP of the book 
 * "Fortran77による数値計算ソフトウェア" by M.Sugihara, M.Mori, published by Maruzen K.K.
 *
 *	@retval result of integration.
 */
template<class func> double mgGausp(
	func& f,		///<Function object to evaluate the function f(t).
				///<f(t) must return double value.
	double a,	///<lower bound of integration.
	double b,	///<upper bound of integration.
	int n=10	///<number of points of the formula, must not be larger than 16.
){
    /* Initialized data */
    const static struct {
	double e_1;
	double fill_2[17];
	double e_3;
	double fill_4[15];
	double e_5[2];
	double fill_6[16];
	double e_7[2];
	double fill_8[14];
	double e_9[3];
	double fill_10[15];
	double e_11[3];
	double fill_12[13];
	double e_13[4];
	double fill_14[14];
	double e_15[4];
	double fill_16[12];
	double e_17[5];
	double fill_18[13];
	double e_19[5];
	double fill_20[11];
	double e_21[6];
	double fill_22[12];
	double e_23[6];
	double fill_24[10];
	double e_25[7];
	double fill_26[11];
	double e_27[7];
	double fill_28[9];
	double e_29[8];
	double fill_30[10];
	double e_31[8];
	double fill_32[8];
	} equiv_7 = { 0., {0}, .5773502691896257, {0}, 0., .7745966692414832, 
		{0}, .3399810435848563, .8611363115940525, {0}, 0., 
		.5384693101056831, .9061798459386639, {0}, .2386191860831969, 
		.6612093864662644, .932469514203152, {0}, 0., 
		.4058451513773972, .7415311855993944, .9491079123427584, {0}, 
		.1834346424956498, .5255324099163289, .7966664774136267, 
		.9602898564975361, {0}, 0., .3242534234038089, 
		.6133714327005903, .8360311073266357, .968160239507626, {0}, 
		.1488743389816312, .4333953941292473, .6794095682990244, 
		.8650633666889845, .9739065285171717, {0}, 0., 
		.269543155952345, .5190961292068118, .7301520055740493, 
		.8870625997680953, .9782286581460569, {0}, .1252334085114689, 
		.3678314989981802, .5873179542866174, .7699026741943046, 
		.9041172563704747, .9815606342467192, {0}, 0., 
		.2304583159551348, .4484927510364469, .6423493394403402, 
		.8015780907333098, .9175983992229779, .9841830547185881, {0}, 
		.1080549487073436, .3191123689278898, .515248636358154, 
		.6872929048116854, .827201315069765, .9284348836635735, 
		.9862838086968123, {0}, 0., .2011940939974345, 
		.3941513470775634, .5709721726085388, .72441773136017, 
		.8482065834104272, .9372733924007058, .9879925180204854, {0}, 
		.09501250983763744, .2816035507792589, .4580167776572274, 
		.6178762444026437, .7554044083550029, .8656312023878317, 
		.9445750230732326, .9894009349916499 };

#define c ((double *)&equiv_7)

    const static struct {
	double e_1;
	double fill_2[17];
	double e_3;
	double fill_4[15];
	double e_5[2];
	double fill_6[16];
	double e_7[2];
	double fill_8[14];
	double e_9[3];
	double fill_10[15];
	double e_11[3];
	double fill_12[13];
	double e_13[4];
	double fill_14[14];
	double e_15[4];
	double fill_16[12];
	double e_17[5];
	double fill_18[13];
	double e_19[5];
	double fill_20[11];
	double e_21[6];
	double fill_22[12];
	double e_23[6];
	double fill_24[10];
	double e_25[7];
	double fill_26[11];
	double e_27[7];
	double fill_28[9];
	double e_29[8];
	double fill_30[10];
	double e_31[8];
	double fill_32[8];
	} equiv_8 = { 2., {0}, .9999999999999998, {0}, .8888888888888888, 
		.5555555555555558, {0}, .6521451548625459, .3478548451374539, 
		{0}, .5688888888888888, .4786286704993666, .2369268850561892, 
		{0}, .4679139345726909, .3607615730481386, .1713244923791703, 
		{0}, .4179591836734694, .3818300505051189, .2797053914892767, 
		.1294849661688699, {0}, .3626837833783622, .3137066458778873, 
		.2223810344533746, .1012285362903763, {0}, .3302393550012598, 
		.3123470770400028, .2606106964029356, .1806481606948573, 
		.08127438836157435, {0}, .2955242247147527, .2692667193099963,
		 .2190863625159819, .1494513491505806, .06667134430868799, {0}
		, .2729250867779006, .2628045445102467, .2331937645919905, 
		.1862902109277342, .1255803694649046, .05566856711617373, {0},
		 .2491470458134029, .2334925365383548, .2031674267230659, 
		.1600783285433463, .1069393259953185, .04717533638651187, {0},
		 .2325515532308739, .2262831802628972, .2078160475368886, 
		.1781459807619456, .1388735102197873, .09212149983772848, 
		.04048400476531587, {0}, .2152638534631578, .2051984637212956,
		 .1855383974779379, .1572031671581936, .1215185706879031, 
		.08015808715976016, .03511946033175195, {0}, 
		.2025782419255613, .1984314853271116, .1861610000155623, 
		.166269205816994, .1395706779261542, .1071592204671719, 
		.07036604748810814, .0307532419961171, {0}, .1894506104550686,
		 .1826034150449236, .1691565193950025, .1495959888165767, 
		.1246289712555339, .09515851168249285, .06225352393864789, 
		.0271524594117541 };

#define w ((double *)&equiv_8)

    if (n < 1 || n > 16) return 0.;

    double c1 = (b + a) / 2, c2 = (b - a) / 2;

   int nh;
   double v;
   if(n % 2 == 0){
		nh = n / 2;
		v = 0.;
    }else{
		nh = (n - 1) / 2;
		v = w[n * 17 - 17] * f(c1);
    }

    for (int i = 1; i <= nh; ++i) {
		double t1 = c1 + c2 * c[i + n * 17 - 17];
		double t2 = c1 - c2 * c[i + n * 17 - 17];
		v += w[i + n * 17 - 17] * (f(t1) + f(t2));
    }

    v = c2 * v;
    return v;
}

#undef w
#undef c

/** @} */ // end of ALGORITHM group
#endif