Analysis.cpp

#pragma once
#include "stdafx.h"

double epsilon;

/*=================
DIFFERENTIATION
==================*/

double nthDerivative1D(double D, double x, Function f)
{   //where D is the number of derivatives to take
	if (D == 0) { return f.evaluate(x); }
	if (D == 1) { epsilon = pow(10, -1 * (10 - getLength(x) - D)); return (0.5 / epsilon)*(f.evaluate(x + epsilon) - f.evaluate(x - epsilon)); }//use secant method

	double answer = 0;
	epsilon = pow(10, -1 * (7 - getLength(x) - D));
	double lastbit = (1 / pow(epsilon, D));

	if (D != floor(D)) {	//fractional derivative
		epsilon = 0.01;
		double derivs = (x / epsilon);
		lastbit = (1 / pow(epsilon, D));
		for (double i = 0; i <= derivs; ++i) {
			double temp = pow(-1, i) * (tgamma(D + 1) / (factorial(i)*tgamma(D + 1 - i)))  * (f.evaluate((x + ((D / 2) - i)*epsilon)));   //limit def derivative
			answer += temp;
		}
		answer *= lastbit;
		return answer;
	}

	for (int i = 0; i <= D; ++i) {
		double temp = pow(-1, i)*combination(D, i) * (f.evaluate((x + ((ceil(D / 2) - i)*epsilon))));   //limit def derivative
		answer += temp;
	}
	answer *= lastbit;
	return answer;
}

ComplexNumber nthDerivative1D(double D, ComplexNumber x, ComplexFunction f)//for ComplexFunctions that are univariate 
{   //where D is the number of derivatives to take
	if (D == 0) { return f.evaluateComplex(x); }
	if (D == 1) { //use secant method
		epsilon = pow(10, -1 * (10 - getLength(x.Re) - D)); 
		ComplexNumber c1(x + ComplexNumber(epsilon));
		ComplexNumber c2(x - ComplexNumber(epsilon));
		return ComplexNumber(0.5 / epsilon)*(f.evaluateComplex(c1) - f.evaluateComplex(c2)); 
	}

	ComplexNumber answer(0);
	epsilon = pow(10, -1 * (7 - getLength(x.Re) - D));
	double lastbit = (1 / pow(epsilon, D));

	if (D != floor(D)) {	//fractional derivative
		epsilon = 0.01;
		double derivs = (x.Re / epsilon);
		lastbit = (1 / pow(epsilon, D));
		for (double i = 0; i <= derivs; ++i) {
			ComplexNumber temp = ComplexNumber(pow(-1, i) * (tgamma(D + 1) / (factorial(i)*tgamma(D + 1 - i))))  * (f.evaluateComplex((x + ComplexNumber((D / 2) - i)*ComplexNumber(epsilon))));   //limit def derivative
			answer = answer + temp;
		}
		answer = answer * lastbit;
		return answer;
	}

	for (int i = 0; i <= D; ++i) {
		ComplexNumber temp = ComplexNumber(pow(-1, i)*combination(D, i)) * (f.evaluateComplex((x + ComplexNumber((ceil(D / 2) - i)*epsilon))));   //limit def derivative
		answer = answer + temp;
	}
	answer = answer * lastbit;
	return answer;
}

double nthDerivative1DEpsilon(double D, double x, double e, Function f) {  //where D is the number of derivatives to take, e is the epsilon value
	if (D == 0) { return f.evaluate(x); }
	if (D == 1) { return (0.5 / epsilon)*(f.evaluate(x + e) - f.evaluate(x - e)); }//use secant method

	double answer = 0;
	double lastbit = (1 / pow(e, D));

	if (D != floor(D)) {	//fractional derivative
		double derivs = (x / e);
		lastbit = (1 / pow(e, D));
		for (double i = 0; i <= derivs; ++i) {
			double temp = pow(-1, i) * (tgamma(D + 1) / (factorial(i)*tgamma(D + 1 - i)))  * (f.evaluate((x + ((D / 2) - i)*e)));   //limit def derivative
			answer += temp;
		}
		answer *= lastbit;
		return answer;
	}

	for (int i = 0; i <= D; ++i) {
		double temp = pow(-1, i)*combination(D, i) * (f.evaluate((x + ((ceil(D / 2) - i)*e))));   //limit def derivative
		answer += temp;
	}
	answer *= lastbit;
	return answer;
}

double nthPartialDerivative(double D, int wrt, std::vector<double> x, Function f) {   //where D is the number of derivatives to take, wrt is the index of the differentiating variable
	if (D == 0) { return f.evaluate(x); }
	if (D == 1) {
		epsilon = pow(10, -1 * (9 - D));
		std::vector<double> x2 = x;
		x[wrt] += epsilon;
		x2[wrt] -= epsilon;
		return (0.5 / epsilon)*(f.evaluate(x) - f.evaluate(x2));
	}//use secant method

	double answer = 0;
	epsilon = pow(10, -1 * (9 - D));
	double lastbit = (1 / pow(epsilon, D));

	if (D != floor(D)) {	//fractional derivative
		double derivs = (x[wrt] / epsilon);
		while (derivs > 100) { derivs /= 10; }
		for (double i = 0; i <= derivs; ++i) {
			std::vector<double> x2 = x;
			x2[wrt] += ((ceil(D / 2) - i)*epsilon);
			double temp = pow(-1, i) * (tgamma(D + 1) / (factorial(i)*tgamma(D = 1 - i)))  * (f.evaluate(x2));   //limit def derivative
			answer += temp;
		}
		answer *= lastbit;
		return answer;
	}

	for (int i = 0; i <= D; ++i) {
		std::vector<double> x2 = x;
		x2[wrt] += ((ceil(D / 2) - i)*epsilon);
		double temp = pow(-1, i)*combination(D, i) * (f.evaluate(x2));   //limit def derivative
		answer += temp;
	}
	answer *= lastbit;
	return answer;
}

/*================
	INTEGRATION
=================*/

std::vector<double> GaussLegendreTable() {
	double x[100] = {
		//  GIANT TABLE OF WEIGHTS AND ABSCISSA POINTS
		//  FOR GAUSS QUADRATURE INTEGRAL OF FUNCTION
		//  NOTE:  Also, include the negative value for
		//  x_i for a total of 100 nodes.
		//  ============================================
		//			x_i								w_i
		0.0156289844215430828722167,	0.0312554234538633569476425,
		0.0468716824215916316149239,	0.0312248842548493577323765,
		0.0780685828134366366948174,	0.0311638356962099067838183,
		0.1091892035800611150034260,	0.0310723374275665165878102,
		0.1402031372361139732075146,	0.0309504788504909882340635,
		0.1710800805386032748875324,	0.0307983790311525904277139,
		0.2017898640957359972360489,	0.0306161865839804484964594,
		0.2323024818449739696495100,	0.0304040795264548200165079,
		0.2625881203715034791689293,	0.0301622651051691449190687,
		0.2926171880384719647375559,	0.0298909795933328309168368,
		0.3223603439005291517224766,	0.0295904880599126425117545,
		0.3517885263724217209723438,	0.0292610841106382766201190,
		0.3808729816246299567633625,	0.0289030896011252031348762,
		0.4095852916783015425288684,	0.0285168543223950979909368,
		0.4378974021720315131089780,	0.0281027556591011733176483,
		0.4657816497733580422492166,	0.0276611982207923882942042,
		0.4932107892081909335693088,	0.0271926134465768801364916,
		0.5201580198817630566468157,	0.0266974591835709626603847,
		0.5465970120650941674679943,	0.0261762192395456763423087,
		0.5725019326213811913168704,	0.0256294029102081160756420,
		0.5978474702471787212648065,	0.0250575444815795897037642,
		0.6226088602037077716041908,	0.0244612027079570527199750,
		0.6467619085141292798326303,	0.0238409602659682059625604,
		0.6702830156031410158025870,	0.0231974231852541216224889,
		0.6931491993558019659486479,	0.0225312202563362727017970,
		0.7153381175730564464599671,	0.0218430024162473863139537,
		0.7368280898020207055124277,	0.0211334421125276415426723,
		0.7575981185197071760356680,	0.0204032326462094327668389,
		0.7776279096494954756275514,	0.0196530874944353058653815,
		0.7968978923903144763895729,	0.0188837396133749045529412,
		0.8153892383391762543939888,	0.0180959407221281166643908,
		0.8330838798884008235429158,	0.0172904605683235824393442,
		0.8499645278795912842933626,	0.0164680861761452126431050,
		0.8660146884971646234107400,	0.0156296210775460027239369,
		0.8812186793850184155733168,	0.0147758845274413017688800,
		0.8955616449707269866985210,	0.0139077107037187726879541,
		0.9090295709825296904671263,	0.0130259478929715422855586,
		0.9216092981453339526669513,	0.0121314576629794974077448,
		0.9332885350430795459243337,	0.0112251140231859771172216,
		0.9440558701362559779627747,	0.0103078025748689695857821,
		0.9539007829254917428493369,	0.0093804196536944579514182,
		0.9628136542558155272936593,	0.0084438714696689714026208,
		0.9707857757637063319308979,	0.0074990732554647115788287,
		0.9778093584869182885537811,	0.0065469484508453227641521,
		0.9838775407060570154961002,	0.0055884280038655151572119,
		0.9889843952429917480044187,	0.0046244500634221193510958,
		0.9931249370374434596520099,	0.0036559612013263751823425,
		0.9962951347331251491861317,	0.0026839253715534824194396,
		0.9984919506395958184001634,	0.0017093926535181052395294,
		0.9997137267734412336782285,	0.0007346344905056717304063
	};
	std::vector<double> x2 = toSTLVector(x, 100);
	return x2;
}

std::vector<double> GaussHermiteTable() {	//table of nodes and weights (x_i and w_i) for Gauss-Hermite quadratre, use for integrals of bound (-inf, +inf).
	double x[150] = {
		//  GIANT TABLE OF WEIGHTS AND ABSCISSA POINTS FOR GAUSS QUADRATURE INTEGRAL OF FUNCTION
		//  NOTE:  These exclusively positive values are saved, but their negative versions can
		//  be used in integration.  This method uses 100 nodes by this method, not 50.
		//  ====================================================================================
		//			x_i					w_i					w_i*e^-x_i^2
		13.4064873381449101385,	5.9080678650312068153E-79,	0.67435355242090814391,
		12.82379974948780890634,	1.9728605748794525545E-72,	0.5185068072702001877,
		12.3429642228596742951,	3.083028990003274812E-67,	0.4499931710544165908,
		11.91506194311416580198,	9.01922230369355618E-63,	0.40868865844190154715,
		11.52141540078703024169,	8.518883081761633787E-59,	0.38013743880057595596,
		11.1524043855851252649,	3.459477936475550445E-55,	0.35881840907111208345,
		10.80226075368471459482,	7.1915294634633710296E-52,	0.34208926231424714546,
		10.46718542134281214178,	8.59756395482527161E-49,	0.32849948464545747856,
		10.14450994129284546989,	6.420725205348472483E-46,	0.31717511097185385938,
		9.832269807777969094355,	3.1852178778359179307E-43,	0.307552728501673729425,
		9.528965823390114804697,	1.10047068271422366946E-40,	0.299249958039079093578,
		9.233420890219161550478,	2.748784884357112492E-38,	0.2919965639602760822598,
		8.944689217325474478799,	5.1162326043852221804E-36,	0.285595214531367884937,
		8.661996168134517714376,	7.27457259688776757451E-34,	0.27989787254508465858,
		8.384696940416265075086,	8.0674342787093771737E-32,	0.274790938326829146069,
		8.112247311162791917212,	7.10181222638493422966E-30,	0.2701855435339907877,
		7.844182384460821168792,	5.0377911662131877842E-28,	0.266011005285956696413,
		7.580100807857488884286,	2.9173500726293324178E-26,	0.262210289443758959518,
		7.319652822304535316332,	1.39484152606876708048E-24,	0.258736790904380865771,
		7.062531060248865437466,	5.5610269616591673171E-23,	0.25555200056201943528,
		6.808463352858796414479,	1.86499767513025225817E-21,	0.2526237833914175732842,
		6.55720703192153931598,	5.3023161831318486853E-20,	0.249925086601289282222,
		6.308544361112135121639,	1.28683292112115327575E-18,	0.247432956126835293416,
		6.062278832614302638665,	2.68249216476037608E-17,	0.2451277779135678348954,
		5.818232135203517047362,	4.8298353217030333477E-16,	0.242992685579107730478,
		5.576241649329924103303,	7.54889687791524329227E-15,	0.24101309292238722418,
		5.336158360138360497277,	1.028874937350992546798E-13,	0.239176321299559035959,
		5.097845105089136247,	1.22787851441012497003E-12,	0.2374712999204616475105,
		4.861175091791210210046,	1.28790382573155823281E-11,	0.235888322794699254658,
		4.626030635787155773074,	1.19130063492907294982E-10,	0.23441885012175822832,
		4.392302078682684016745,	9.747921253871621245759E-10,	0.233055344869710310074,
		4.159886855131030540068,	7.0758572838895729074E-9,	0.231791137453703842101,
		3.92868868342767097201,	4.5681275084849395136E-8,	0.230620313034512848758,
		3.698616859318491939797,	2.62909748375372507934E-7,	0.229537617164869894481,
		3.469585636418589169768,	1.3517971591103672866E-6,	0.22853837642624680354,
		3.241513679631012950359,	6.2215248177778633174E-6,	0.22761843139840206564,
		3.014323580331155516715,	2.5676159384549063056E-5,	0.22677407984368572366,
		2.78794142398198931319,	9.5171627785509664702E-5,	0.226002028407800085009,
		2.562296402372608025056,	3.172919710433003055501E-4,	0.2252993514677081961876,
		2.337320463906878505013,	9.52692188548619117479E-4,	0.2246634560172451970005,
		2.11294799637118795203,	0.002579273260059090173396,	0.224092051687791109734,
		1.889115537427008371494,	0.00630300028560805254939,	0.2235831251671776091963,
		1.665761508741509469867,	0.01391566522023180641788,	0.2231349184139684965647,
		1.442825970215932787703,	0.0277791273859335142701,	0.222745910173357235135,
		1.22025039121895305882,	0.0501758126774286956979,	0.2224148003905443508997,
		0.9979774360981052439241,	0.0820518273912244646804,	0.2221404971917734901661,
		0.7759507615401457819751,	0.1215379868441041819837,	0.2219221061654981101362,
		0.554114823591616988233,	0.1631300305027829414193,	0.221758921729023785805,
		0.3324146923422318070459,	0.1984628502541864777089,	0.22165042041151554711,
		0.1107958724224394828876,	0.2188926295874391250627,	0.2215962559241832702211,
	};
	std::vector<double> x2 = toSTLVector(x, 150);
	return x2;
}

std::vector<long double> GaussLaguerreTable() {//table of nodes and weights (x_i and w_i) for Gauss-Laguerre quadratre, use for integrals of bound (-inf, +inf).
	long double x[256] = {
		//  GIANT TABLE OF WEIGHTS AND ABSCISSA POINTS FOR GAUSS QUADRATURE INTEGRAL OF FUNCTION
		//  NOTE:  This method uses only positive values, for a total of 100 nodes. ALPHA = 0
		//  ====================================================================================
		//x_i												w_i													  w_i/w(x)	
		1.125138826367596296085184031621e-02, 2.855184445323972862907317736119e-02,
		5.928474126902645428792200896141e-02, 6.335021178450511877979781272594e-02,
		1.457079665943124651418540591021e-01, 9.130838136613431442316163259027e-02,
		2.705531787586650661907608971002e-01, 1.099139004109117461010139158327e-01,
		4.338414075538368030560965807545e-01, 1.182741710341736987898096888739e-01,
		6.355976657816219383408676779691e-01, 1.170457390004067215664584392067e-01,
		8.758523845465207791553460132612e-01, 1.080899875455684154364737831251e-01,
		1.154641701974397950081533557076e+00, 9.394288863892850178780884363559e-02,
		1.472007563166735475544466330382e+00, 7.725366879789805320778002523588e-02,
		1.827997778312350285359845287177e+00, 6.032705626566157053893030860025e-02,
		2.222666071561902448178964529136e+00, 4.484734824719521406824246579976e-02,
		2.656072129883481195228853293088e+00, 3.179694793687684617396324848206e-02,
		3.128281655027916954983103697375e+00, 2.153014945379444392611072854376e-02,
		3.639366419852403212210745221687e+00, 1.393695173384634832775768859749e-02,
		4.189404329594047974784930798646e+00, 8.631585380202247148844730964892e-03,
		4.778479488434876091657242392131e+00, 5.117777013669228528739367228449e-03,
		5.406682271600499185278938201046e+00, 2.906347436485955858179800772191e-03,
		6.074109403196536843091558445058e+00, 1.581432943316679397234160134889e-03,
		6.780864039975625411048049299433e+00, 8.247389850988121504354385932530e-04,
		7.527055861229375885855122798418e+00, 4.123260885396947309702908308044e-04,
		8.312801165007770603378843811915e+00, 1.976494264425914986205298897833e-04,
		9.138222970880392396002626419694e+00, 9.085157887824515080228263061528e-05,
		1.000345112946822208926827614353e+01, 4.004849278358052989778876604415e-05,
		1.090862243899088254784886130100e+01, 1.693076239808178557551028884746e-05,
		1.185388076909185680383326445378e+01, 6.864525291110682089386362784119e-06,
		1.283937719222324968745519356727e+01, 2.669216598142102660158722285839e-06,
		1.386527012289208030295367999714e+01, 9.953640102863844771774833321963e-07,
		1.493172546509192747374735531331e+01, 3.559435753003065439880201665633e-07,
		1.603891676826707935092137834277e+01, 1.220532551948811948312056157344e-07,
		1.718702539218126510275850445906e+01, 4.012791920935635068901677660238e-08,
		1.837624068108969493338275233695e+01, 1.264811414747597864456501109076e-08,
		1.960676014764164672790546909894e+01, 3.821489729426572290234113720030e-09,
		2.087878966697137291589320144029e+01, 1.106641059227341699944800440243e-09,
		2.219254368146783690667639231820e+01, 3.071009237097423195822900346389e-10,
		2.354824541674892056092497300970e+01, 8.165499384154489560264378850036e-11,
		2.494612710940348862795103966404e+01, 2.079853632781377842341896121163e-11,
		2.638643024710529089762693051319e+01, 5.073953770839870404329698640195e-12,
		2.786940582174639022956968185640e+01, 1.185314377179611230509373313098e-12,
		2.939531459628491372156562883813e+01, 2.650927523728873586005654881953e-13,
		3.096442738605275400232208613174e+01, 5.674632215757658766810656061614e-14,
		3.257702535532375017814194864562e+01, 1.162373814707515892215294349006e-14,
		3.423340033000224266047941087526e+01, 2.277766292702386379197331044506e-15,
		3.593385512735615381077229249628e+01, 4.268831970297649277391721041258e-16,
		3.767870390378837443006555820165e+01, 7.649288799363275105259484578031e-17,
		3.946827252171576412714890330036e+01, 1.310131391983824641880828868213e-17,
		4.130289893670708963800804176367e+01, 2.144145234124663634370678869197e-18,
		4.318293360612038324387836352250e+01, 3.351944287208847808014707290438e-19,
		4.510873992057723174415061485072e+01, 5.003733086459473768231793651206e-20,
		4.708069465971721685607253511276e+01, 7.130030641958562120497024646264e-21,
		4.909918847379102680215358528598e+01, 9.694540740397266403532090582853e-22,
		5.116462639277665944464043359157e+01, 1.257284755639784598440599274317e-22,
		5.327742836484077391610853679441e+01, 1.554661095563063448220273119882e-23,
		5.543802982611789186830896382911e+01, 1.832109793253421778719084254378e-24,
		5.764688230394522881442498112197e+01, 2.056797978136734920722781371545e-25,
		5.990445405587205569652926350617e+01, 2.198666052623291192576574486891e-26,
		6.221123074696145824567915529620e+01, 2.236916007324289367294062222812e-27,
		6.456771616812121542904105154672e+01, 2.164960644633905440025630940452e-28,
		6.697443299844156105480271561954e+01, 1.992227680618793787387725075369e-29,
		6.943192361478342995576210977421e+01, 1.742153886325439585907653178392e-30,
		7.194075095215437515730184810623e+01, 1.446949786106284637699605377666e-31,
		7.450149941873402779302798318554e+01, 1.140751706123082283418897029015e-32,
		7.711477586977057052831989243538e+01, 8.531805097810209072211639296988e-34,
		7.978121064496855285445821249912e+01, 6.049701177938858435049974561125e-35,
		8.250145867443145291403917688451e+01, 4.064343243264800301779467322272e-36,
		8.527620065871535873779640425817e+01, 2.585349374987909630702838798761e-37,
		8.810614432909950369403173932576e+01, 1.556028762522623447585274555958e-38,
		9.099202579479261315603030302452e+01, 8.854625849663330011025592370634e-40,
		9.393461098447969449552446429249e+01, 4.760457517364580680322165132837e-41,
		9.693469719038194045161995512395e+01, 2.416078510661232205003307206936e-42,
		9.999311472386427152162130002668e+01, 1.156647050338737493210994431814e-43,
		1.031107286925939873923197491579e+02, 5.218510619492375995216995904088e-45,
		1.062884409103454426684261296594e+02, 2.216974335336180330548414691337e-46,
		1.095271919517775508066180569177e+02, 8.860102756613696060255660465478e-48,
		1.128279633659041938773334872636e+02, 3.327811159201095553289565856322e-49,
		1.161917800635567809402358717078e+02, 1.173490043078302543621673297784e-50,
		1.196197128959320104623488874197e+02, 3.880967726420921430711844828839e-52,
		1.231128814433601900609118145086e+02, 1.202426327933061417922628432257e-53,
		1.266724570357601836623386949566e+02, 3.486023044105546382152946824426e-55,
		1.302996660289134625872174928638e+02, 9.445545221595566814883993731545e-57,
		1.339957933637479643435821208364e+02, 2.388884274559683949225854100867e-58,
		1.377621864393393809643026665775e+02, 5.631884750754630520125877615504e-60,
		1.416002593343930403057896427217e+02, 1.235928611912160188213664051764e-61,
		1.455114974166595923936405970079e+02, 2.521004202377267433009913724667e-63,
		1.494974623851777070881731754510e+02, 4.772224621999805220517054020724e-65,
		1.535597977965664401179827482608e+02, 8.370019891999578341619585259759e-67,
		1.577002351339781050590953365465e+02, 1.357824341120209853637619284541e-68,
		1.619206004859756341637536290312e+02, 2.033688727153154158744017519558e-70,
		1.662228219127688750928757391598e+02, 2.806838480695353838788461214029e-72,
		1.706089375892422346465506633095e+02, 3.562567607062095973328239094080e-74,
		1.750811048284146048806174055016e+02, 4.149452749293770593460266457155e-76,
		1.796416101058669946026349646393e+02, 4.425007965766321930949020891633e-78,
		1.842928802258468053418345050204e+02, 4.310084261289849665223321554111e-80,
		1.890374947939541090012929983448e+02, 3.824661016761739822392023682989e-82,
		1.938782001904729675408028759397e+02, 3.083547842598792754625941584836e-84,
		1.988179252737206028047459449936e+02, 2.252139822170620842056423253276e-86,
		2.038597990857698445716648248975e+02 ,1.485514740645043123171789185206e-88,
		2.090071708855108678533875111813e+02, 8.819635476372656372128494267554e-91,
		2.142636328987880212805277584919e+02, 4.696417822125985072051048332170e-93,
		2.196330462555781740383874010243e+02, 2.234393825454772742691079434952e-95,
		2.251195706842090277566597965661e+02, 9.458787038220740320973770981125e-98,
		2.307276986582036196816588686801e+02, 3.546960831240672613931702385830e-100,
		2.364622948501776659660189041578e+02, 1.172532130034887230041858420158e-102,
		2.423286419497026982678645198658e+02, 3.399090555639915547941644793682e-105,
		2.483324941623571788920166017798e+02, 8.591907200623898045079541036460e-108,
		2.544801400448691318935438033580e+02, 1.881891397353535964729395644556e-110,
		2.607784767735797365385600645383e+02, 3.547358632306256523712374269015e-113,
		2.672350985289538367639924720288e+02, 5.711482228283600474459214964334e-116,
		2.738584024626936097934146026484e+02, 7.789473788044460956107072732965e-119,
		2.806577167763234923975041009771e+02, 8.915898699491269351479711891865e-122,
		2.876434568992193306384736778999e+02, 8.476856358868403207417548329016e-125,
		2.948273177876477391798066721040e+02, 6.617326935494900345407717728242e-128,
		3.022225132464494653805359817114e+02, 4.186216357415709519007664178836e-131,
		3.098440773266126634477723636429e+02, 2.114385168981142071200931970133e-134,
		3.177092489549062894956780523400e+02, 8.382163501367869536416750179622e-138,
		3.258379701211949496504012779308e+02, 2.557202302197677884687797768369e-141,
		3.342535420676541353751844501744e+02, 5.866768642191204372046123603405e-145,
		3.429835062738253164085089133285e+02, 9.849861030064843801988568904805e-149,
		3.520608535465261850830434269835e+02, 1.171383943342068942456857274084e-152,
		3.615257263923250475990668518389e+02, 9.483963265567383663821702301497e-157,
		3.714278892143275230785179848674e+02, 4.977096381123802811665397634330e-161,
		3.818304441190610801962076168823e+02, 1.590898527750997654819806386954e-165,
		3.928156712408080988093778198985e+02, 2.856303829119002923206075680445e-170,
		4.044947247505150743896660716596e+02, 2.582250719691489992650314591223e-175,
		4.170249029779890158201972775939e+02, 1.007350250050797409529831872550e-180,
		4.306434441665973815583235516679e+02, 1.344252500443816318217729833630e-186,
		4.457430969739279896521717207263e+02, 4.182962214036834733897266272207e-193,
		4.630800341094462582080137934057e+02, 1.457165307726186315944816631883e-200,
		4.846155439864439760440631311098e+02, 8.640591690468708676928914223540e-210
	};
	std::vector<long double> x2 = toSTLVector(x, 256);
	return x2;
}

std::vector<double> GaussLaguerreTableWithAlpha() {//table of nodes and weights (x_i and w_i) for Gauss-Laguerre quadratre, use for integrals of bound (-inf, +inf).
	double x[300] = {
		//  GIANT TABLE OF WEIGHTS AND ABSCISSA POINTS FOR GAUSS QUADRATURE INTEGRAL OF FUNCTION
		//  NOTE:  This method uses only positive values, for a total of 128 nodes. ALPHA = 3
		//  ====================================================================================
		//x_i						w_i						  w_i/w(x)	
		0.099782081069278731671,	9.9101125967843946E-5,	0.11021881154909085,
		0.23357552071597127401,		0.00159063146538549611,	0.157662406767508028,
		0.41533917986659868788,		0.00973898600796234939,	0.205913373050166629,
		0.64546122544105426427,		0.0358693267435527446,	0.254350736289247745,
		0.92407392100409047388,		0.0948602631482722415,	0.3028896185371095304,
		1.251269685395905049491,	0.197051212319999157,	0.351516618536502434,
		1.627137547496871919,		0.33879027754366553552,	0.40023493709983758699,
		2.051772781608406001,		0.4984387300877394884,	0.4490529825496380824,
		2.525280221996693004469,	0.641839441421756335,	0.497981185619600222,
		3.047775659335826366678,	0.7350700312152139973,	0.5470309252709534374,
		3.619386553175121876902,	0.7575783483828034958,	0.5962141188319223784,
		4.240252474260967089977,	0.7088953209706128649,	0.6455430540543899113,
		4.91052543507886991805,		0.606420041118840738,	0.695030319778331347,
		5.630370175871896587881,	0.4768117614680361653,	0.7446887806901526023,
		6.39996443732240497393,		0.3460784865157818125,	0.7945315734810861967,
		7.21949923557392243746,		0.2326854793201042904,	0.8445721142510163469,
		8.089179148146018484695,	0.1453333124608997458,	0.894824112345974061,
		9.00922261585072538749,		0.08452334471101406683,	0.9453015882492854973,
		9.979862264092981464119,	0.045860850857009976982,0.9960188943182878783,
		11.00134524605888180235,	0.02325180584014413378,	1.0469907377531165426,
		12.0739336098562341774,		0.011030565657988878249,1.0982322055007874619,
		13.19790469147190957022,	0.004901675107654013331,1.1497587909750855,
		14.37355153534770431968,	0.002042175993162277856,1.201586422577618271,
		15.60118334439720879368,	7.983058651357365944E-4,1.2537314940715985633,
		16.88112596136162867306,	2.929798015489975228E-4,1.3062108969051761372,
		18.21372238351696351034,	1.00997333520519311799E-4,	1.35904205461506932827,
		19.59933331289004514927,	3.2715307576848360707E-5,1.4122429594692722772,
		21.03833774431250005645,	9.96057795891044773519E-6,1.46583221153292348618,
		22.53113359383806907496,	2.85099779493774203816E-6,1.51982906036598840035,
		24.07813837026984443496,	7.6725679369126769334E-7,1.5742534495865350869,
		25.67978989279091566825,	1.94150289949960210198E-7,	1.62912606456000219821,
		27.33654705796652795469,	4.6193651467770271145E-8,1.6844683835037115549,
		29.04889065969070960669,	1.03333357049606930633E-8,	1.7403027323276111678,
		30.81732426598861046831,	2.17297613606434979747E-9,	1.79665234356746843535,
		32.64237515696132593745	,4.2948118613088091042E-10,	1.8535414198061044213,
		34.52459532857723716918,	7.9763489622468595281E-11,	1.91099520202245477279,
		36.46456256747807434478,	1.39158348004958596225E-11,	1.96904004335802407145,
		38.46288160248501641946,	2.27988354720336844794E-12,	2.0277034888465337316,
		40.52018533906711481998	,3.506281729265240310588E-13,	2.087014361716243190853,
		42.63713618367917789197,	5.0596879558506090623E-14,	2.1470028569467017796,
		44.81442746559822916248,	6.8475495460660578626E-15,	2.20770064284390636419,
		47.05278496469746038567,	8.68665573040780409288E-16,	2.26914097149156862938,
		49.35296855450664884186,	1.03235020766015589065E-16,	2.33135879904328161833,
		51.71577397093271242603,	1.148653074595729443869E-17,	2.394390916942994091744,
		54.1420347181702194518	,1.195769467200037419086E-18,	2.45827609530192711275,
		56.63262412463884244766,	1.16383298753891711585E-19,	2.52305523982195313083,
		59.18845756326586997183,	1.05824027216311573473E-20,	2.58877156384213669703,
		61.81049485211386430307	,8.98196482919520136378E-22,	2.65547077730094413945,
		64.49974285326802270469,	7.1100276413329832905E-23,	2.723201294656759424493,
		67.25725829008211020258,	5.2441970164689876315E-24,	2.79201446410003269245,
		70.08415080538018316268,	3.60050794697161506664E-25,	2.8619648207291552813,
		72.98158628607619524409,	2.298624844394128631E-26,	2.9331103667580948803,
		75.9507904829674832324,	1.36303273374174782539E-27,	3.00551288228797701767,
		78.99305295825578537761,	7.4983029706389601006E-29,	3.0792382707206170536,
		82.10973139774046383519,	3.82200696176891277302E-30,	3.15435694353594644269,
		85.30225632972086386553,	1.802651618975398610936E-31,	3.2309442499175503204,
		88.57213629856860438059,	7.8561308189462436152E-33,	3.30908095761607047248,
		91.92096354784420675989,	3.15885213032958450021E-34,	3.388853792519913073655,
		95.35042027592849095449,	1.16998554455408382049E-35,	3.47035604569497134628,
		98.86228553665331383359,	3.98498040028886572921E-37,	3.553688258207995298074,
		102.4584428686381874369,	1.2459024896164527259E-38,	3.6389589959222214615,
		106.1408887503268089859,	3.5687984696327369925E-40,	3.72628572872523092172,
		109.9117419935144530256,	9.3466150944843386956E-42,	3.81579583141467631308,
		113.7732542070172064594,	2.2332494352118337472E-43,	3.9076277268513597796,
		117.7278214847515653362,	4.8569520028689815528E-45,	4.0019321961503397373,
		121.7779974997437966299,	9.590885576114712386E-47,	4.0988738858251912425,
		125.9265082185840308578,	1.71503014678406870744E-48,	4.19863304819633330399,
		130.1762684909993652459,	2.7693061397288907284E-50,	4.30140755937298969448,
		134.53040081836745625,	4.0256364953738563535E-52,	4.4074152691842049466,
		138.9922566654885255002,	5.250998384289604018E-54,	4.51689675018555394866,
		143.5654407548590235339,	6.1243940074211687239E-56,	4.6301185291345209443,
		148.2538388760925493334,	6.3627873826741540341E-58,	4.7473769052325555617,
		153.0616498603790054758,	5.8642063229286694737E-60,	4.86900248651115931651,
		157.9934225181377921915,	4.7732339872364888973E-62,	4.995365611112620248186,
		163.0540985269941946428	,3.4147059587775692538E-64,	5.12688286683982830098,
		168.2490625001046571486,	2.1356814308688541949E-66,	5.26402498439479981429,
		173.5842007798958737597,	1.16105456095395335552E-68,	5.40732646314702154152,
		179.0659709148782858129,	5.4518961454699877677E-71,	5.55739740168315224668,
		184.7014843232417118325,	2.1957720432869313415E-73,	5.71493816141996101412,
		190.4986053778277342983,	7.52688988395321245484E-76,	5.88075770901668023266,
		196.4660711374518854682,	2.1772254305737710361E-78,	6.05579679065364558716,
		202.6136373095998902803,	5.26364416132325476957E-81,	6.2411575322689533595,
		208.9522579247043192547	,1.05217781064922846145E-83,	6.43814170335563372516,
		215.4943088857363937353,	1.7180225026741450365E-86,	6.64830083821729420811,
		222.2538694243907660249,	2.25994066684589121565E-89,	6.87350285865557578295,
		229.247081179342998443,	2.35725401191294123402E-92,	7.1160220906282304276,
		236.4926131506474225181,	1.91422737211540851185E-95,	7.37866314227616216647,
		244.0122739288232081931,	1.18456639003624991023E-98,	7.664934955701435739,
		251.831833400410306004,	5.4465453981346662057E-102,	7.9793012094569962296,
		259.9821501127195400405,	1.8052388934916472248E-105,	8.3275505107102251628,
		268.5007580425110914794,	4.1574208560099965176E-109,	8.71736130981819037962,
		277.4341682408380130416,	6.3563656585386461405E-113,	9.15919679548899079766,
		286.8413298329903379918,	6.0903962100755566096E-117,	9.6677874331715106096,
		296.7990677823919950087,	3.392638180689950864E-121,	10.2647250206185286954,
		307.4111069834275428207,	9.929902186555398412392E-126,	10.983322725519597443,
		318.824138863590525703,	1.32389387700966592627E-130,	11.8785599413584008405,
		331.2592480196265960103,	6.4871050676259301377E-136,	13.0499963855324313345,
		345.0821517127607298274,	8.2023622573924094412E-142,	14.7043471092340390983,
		360.9963463997363328048,	1.36086554781135209814E-148,	17.37912987104347053464,
		380.824562950739007213,	5.2407105091012422329E-157,	23.2927934794332467763
};
	std::vector<double> x2 = toSTLVector(x, 300);
	return x2;
}

std::vector<double> GaussKronodTable() {//table of nodes and weights (x_i and w_i) for Gauss-Kronod quadrature
	double x[198] = {
		//  GIANT TABLE OF WEIGHTS AND ABSCISSA POINTS FOR GAUSS QUADRATURE INTEGRAL OF FUNCTION
		//  NOTE: there is total of 98 Gauss nodes, 198 Kronod nodes
		//  ====================================================================================
		//Kronrod x_i>				Kronrod w_i					Corresponding Gauss w_i	
		-0.999804199685638727995,	5.274769683783323143E-4,
		-0.9988201506066353793618,	0.0014779639281743620209,	//0.003027278988922905077
		-0.996819814299826469401,	0.0025218765731496496845,
		-0.9937886619441677907601,	0.0035329557014832599892,	//0.0070350995900864514735
		-0.9897642384140644710632,	0.0045141331625998620836,
		-0.984757895914213004359,	0.005501420399338078204,	//0.0110205510315935804975
		-0.9787573030312066845299,	0.0064998725321664845168,
		-0.97176220090155538014,	0.0074869102764140445198,	//0.0149621449356246510296
		-0.9637900181363554282611,	0.0084551856397750456701,
		-0.954853658674137233555,	0.0094175722296862066762,	//0.0188435958530894584445
		-0.944955059221329473156,	0.010378732924116607707,
		-0.934100294755810149059,	0.0113278438578780228795,	//0.0226492015874466764988
		-0.922305475453936168755,	0.0122591748799473589077,
		-0.909585655828073285213,	0.0131792068212079366057,	//0.0263636189270660169609
		-0.8959496245835934207352,	0.0140911110472705440377,
		-0.881408445573008910037,	0.0149880989346802956593,	//0.0299718846205838253507
		-0.8659800043151457264443,	0.01586572369887289313,
		-0.8496821198441657010349,	0.016727898553777318682,	//0.0334594667916221743425
		-0.832528501487460168517,	0.017576872641244826238,
		-0.8145344273598554315395,	0.0184077538519825820281,	//0.03681232096300068981947
		-0.7957203207625248381361,	0.0192169329826556442117,
		-0.7761068943454466350181,	0.0200070643859274929265,	//5.E-47
		-0.7557118903695143082457,	0.0207798535198561693337,
		-0.734554254237402696214,	0.0215314818077816867538,	//5.E-40
		-0.7126570667050088308474,	0.022258913930129946636,
		-0.6900438244251321135048,	0.0229641216311680166298,	//9.E-33
		-0.666735700841571667866,	0.0236484943738631350276,
		-0.642754832419237664057,	0.02430890292981194970869,	//3.E-25
		-0.6181268178008792991612,	0.0249427312116380137806,
		-0.5928776941089007124559,	0.0255515676298983967962,	//1.E-20
		-0.5670315494953917328631,	0.0261366295727170561997,
		-0.5406132469917260665582,	0.0266952731071580377427,	//2.E-14
		-0.5136506284017343823111,	0.027225205808796711996,
		-0.486171941452492042177,	0.0277278075508688728381,	//3.E-9
		-0.4582036915390298548496,	0.0282042205700292005214,
		-0.4297729933415765246586,	0.0286521667109019998496,	//7.E-5
		-0.4009095739292798809373,	0.0290696140121504756931,
		-0.3716435012622848888637,	0.0294578454056420734076,	//0.05
		-0.3420031959918559601078,	0.0298179973237107521896,
		-0.3120175321197487622079,	0.030148081066799933145	,	//0.06
		-0.2817177082410294178775,	0.0304462797396858877386,
		-0.2511351786125772735072,	0.030713855857739083392,	//0.06
		-0.2202997655098053353243,	0.030951992615566392528,
		-0.1892415924618135864853,	0.03115893749659741587,		//0.0623
		-0.157992877664368358666,	0.0313330525912310597589,
		-0.126585997269672051068,	0.03147563582699890351852,	//0.062953
		-0.09505164998612337566,	0.031587959458685973161,
		-0.0634206849826867860288,	0.03166846605661171584462,	//0.06333551
		-0.03172586345907315363082,	0.031715664503973344041,
		0,							0.0317309313985218944091,	//0.0634632814047905977183
		0.03172586345907315363082,	0.031715664503973344041,
		0.06342068498268678602884,	0.0316684660566117158446,	//0.0633355093
		0.09505164998612337566,		0.031587959458685973161,
		0.126585997269672051068,	0.031475635826998903519,	//0.0629527074652
		0.157992877664368358666,	0.0313330525912310597589,
		0.189241592461813586485,	0.03115893749659741587,		//0.062316417320057
		0.2202997655098053353243,	0.03095199261556639252752,
		0.2511351786125772735072,	0.030713855857739083392,	//0.0614292009791929
		0.2817177082410294178775,	0.0304462797396858877386,
		0.3120175321197487622079,	0.030148081066799933145,	//0.06029463095315202
		0.3420031959918559601078,	0.0298179973237107521896,
		0.3716435012622848888637,	0.029457845405642073408,	//0.058917275760027266
		0.4009095739292798809373,	0.0290696140121504756931,
		0.4297729933415765246586,	0.0286521667109019998496,	//0.0573026815301874755
		0.4582036915390298548496,	0.02820422057002920052136,
		0.486171941452492042177,	0.027727807550868872838,	//0.05545734967480358869
		0.5136506284017343823111,	0.027225205808796711996,
		0.540613246991726066558,	0.026695273107158037743,	//0.0533887107082589685279
		0.567031549495391732863,	0.0261366295727170561997,
		0.5928776941089007124559,	0.025551567629898396796,	//0.05110509433014459067462
		0.618126817800879299161,	0.0249427312116380137806,
		0.6427548324192376640569,	0.0243089029298119497087,	//0.04861569588782824027765
		0.666735700841571667866,	0.0236484943738631350276,
		0.6900438244251321135048,	0.02296412163116801663,		//0.0459305393555958535425
		0.712657066705008830847,	0.02225891393012994663599,
		0.7345542542374026962137,	0.0215314818077816867538,	//0.04306043698125959798835
		0.7557118903695143082457,	0.020779853519856169334,
		0.7761068943454466350181,	0.0200070643859274929265,	//0.04001694576637302136861
		0.7957203207625248381361,	0.0192169329826556442117,
		0.81453442735985543154,		0.0184077538519825820281,	//0.0368123209630006898195
		0.8325285014874601685173,	0.017576872641244826238,
		0.8496821198441657010349,	0.016727898553777318682,	//0.0334594667916221743425
		0.8659800043151457264443,	0.01586572369887289313,
		0.881408445573008910037,	0.014988098934680295659,	//0.0299718846205838253507
		0.8959496245835934207352,	0.0140911110472705440377,
		0.909585655828073285213,	0.0131792068212079366057,	//0.0263636189270660169609
		0.9223054754539361687554,	0.012259174879947358908,
		0.934100294755810149059,	0.0113278438578780228795,	//0.0226492015874466764988
		0.944955059221329473156,	0.010378732924116607707,
		0.9548536586741372335552,	0.0094175722296862066762,	//0.0188435958530894584445
		0.9637900181363554282611,	0.0084551856397750456701,
		0.97176220090155538014,		0.0074869102764140445198,	//0.0149621449356246510296
		0.97875730303120668453,		0.0064998725321664845168,
		0.9847578959142130043593,	0.005501420399338078204,	//0.0110205510315935804975
		0.9897642384140644710632,	0.0045141331625998620836,
		0.9937886619441677907601,	0.0035329557014832599892,	//0.007035099590086451473
		0.9968198142998264694013,	0.002521876573149649685,
		0.9988201506066353793618,	0.0014779639281743620209,	//0.003027278988922905077
		0.999804199685638727995,	5.2747696837833231426E-4
	};
	std::vector<double> x2 = toSTLVector(x, 198);
	return x2;
}

std::vector<double> GaussKronodTableCorrespondingGaussTable() {//table of corresponding Gauss weights (x_i and w_i) for Gauss-Kronod quadratre
	double x[98] = {
		-0.9988201506066353793618312727041106551236809978645,	0.00302727898892290507748069817582714097283004200383,
		-0.99378866194416779076011385925795361884130739493865,	0.00703509959008645147345067831134841086286385703903,
		-0.98475789591421300435929899536857471035330746214525,	0.011020551031593580497508288166833332752193560544958,
		-0.97176220090155538013997242150477413349786254718091,	0.01496214493562465102958431912020150403628101637277,
		-0.95485365867413723355524367966818441905469445190027,	0.01884359585308945844445065339110734926774990201925,
		-0.93410029475581014905898246063654490246042295281019,	0.02264920158744667649877096421604288732679749902326,
		-0.90958565582807328521301964035580007104277123111672,	0.02636361892706601696094574523974204769446097176574,
		-0.88140844557300891003703153581700492615924003726239,	0.02997188462058382535069055807414721006207203595088,
		-0.84968211984416570103488187224297933332660433246472,	0.03345946679162217434248715089168920184947986786023,
		-0.81453442735985543153950078776477704449467213582275,	0.036812320963000689819467236698882295001361552406,
		-0.77610689434544663501814281224845607126445046182918,	0.040016945766373021368605036743103321427864607465,
		-0.7345542542374026962136742128585248846528354533237,	0.04306043698125959798834548400103849580671358123982,
		-0.69004382442513211350475180727494192763350512618904,	0.04593053935559585354249961997212767100882432235462,
		-0.6427548324192376640568569485839001095403354722533,	0.04861569588782824027765119728786853088136892861859,
		-0.59287769410890071245586433679262853463300101531885,	0.0511050943301445906746228060493185613896124493297,
		-0.54061324699172606655822546722943591692410134623837,	0.0533887107082589685279429337301762061655860694482,
		-0.48617194145249204217697609605857506207511519433494,	0.05545734967480358869043158148269806720425152407379,
		-0.42977299334157652465858414180742303306341789787807,	0.05730268153018747548516459235081653635058043880312,
		-0.3716435012622848888637345409469605276147447867339,	0.05891727576002726602452765074249251626249611385948,
		-0.31201753211974876220786067650641626509683340799439,	0.06029463095315201730310611687554159049777650891973,
		-0.25113517861257727350715594282461127458298417457366,	0.0614292009791929362968266461403010749431900091921,
		-0.18924159246181358648531017385998651117967887767882,	0.0623164173200572674010768252159497042241335362213,
		-0.12658599726967205106798528857553619032448712106687,	0.06295270746519569947439955578289100636372542776957,
		-0.06342068498268678602883482087848497315278687277888,	0.0633355092964917485908369274000905201479714165648,
		0,														0.06346328140479059771825346787952963499962952230825,
		0.06342068498268678602883482087848497315278687277888,	0.06333550929649174859083692740009052014797141656483,
		0.12658599726967205106798528857553619032448712106687,	0.0629527074651956994743995557828910063637254277696,
		0.18924159246181358648531017385998651117967887767882,	0.06231641732005726740107682521594970422413353622132,
		0.2511351786125772735071559428246112745829841745737,	0.0614292009791929362968266461403010749431900091921,
		0.31201753211974876220786067650641626509683340799439,	0.06029463095315201730310611687554159049777650891973,
		0.37164350126228488886373454094696052761474478673393,	0.0589172757600272660245276507424925162624961138595,
		0.42977299334157652465858414180742303306341789787807,	0.0573026815301874754851645923508165363505804388031,
		0.48617194145249204217697609605857506207511519433494,	0.0554573496748035886904315814826980672042515240738,
		0.54061324699172606655822546722943591692410134623837,	0.05338871070825896852794293373017620616558606944823,
		0.59287769410890071245586433679262853463300101531885,	0.0511050943301445906746228060493185613896124493297,
		0.64275483241923766405685694858390010954033547225328,	0.048615695887828240277651197287868530881368928618593,
		0.69004382442513211350475180727494192763350512618904,	0.04593053935559585354249961997212767100882432235462,
		0.7345542542374026962136742128585248846528354533237,	0.04306043698125959798834548400103849580671358123982,
		0.77610689434544663501814281224845607126445046182918,	0.040016945766373021368605036743103321427864607465,
		0.81453442735985543153950078776477704449467213582275,	0.03681232096300068981946723669888229500136155240598,
		0.84968211984416570103488187224297933332660433246472,	0.033459466791622174342487150891689201849479867860229,
		0.88140844557300891003703153581700492615924003726239,	0.0299718846205838253506905580741472100620720359509,
		0.90958565582807328521301964035580007104277123111672,	0.0263636189270660169609457452397420476944609717657,
		0.93410029475581014905898246063654490246042295281019,	0.0226492015874466764987709642160428873267974990233,
		0.95485365867413723355524367966818441905469445190027,	0.01884359585308945844445065339110734926774990201925,
		0.97176220090155538013997242150477413349786254718091,	0.01496214493562465102958431912020150403628101637277,
		0.98475789591421300435929899536857471035330746214525,	0.011020551031593580497508288166833332752193560545,
		0.99378866194416779076011385925795361884130739493865,	0.00703509959008645147345067831134841086286385703903,
		0.9988201506066353793618312727041106551236809978645,	0.003027278988922905077480698175827140972830042003832,
	};
	std::vector<double> x2 = toSTLVector(x, 98);
	return x2;
}

double integrateSimpsonsRule(double a, double b, Function f) {//fast, imprecise method of integration
	double h = std::abs(b - a) / 2000.0;	//always use 2000 abscissa points -- the system limit for quick integration
	double counter = a;
	int index = 0;
	double answer = 0;
	bool skip = false;
	while (counter < b) {
		if (counter == a || counter == (b - h)) { answer += (f.evaluate(counter)); skip = true; }
		if (counter != a && index % 2 == 0 && skip == false) { answer += (2 * f.evaluate(counter)); skip = true; }
		if (skip == false) { answer += (4 * f.evaluate(counter)); }
		counter += h;
		++index;
		skip = false;
	}
	answer *= (h / 3);
	return answer;
}

double integrateGaussLegendreQuadrature(double a, double b, Function f) {//uses Chebyshev-Gauss quadrature rules -- slower but more precise
	const int n = 100;	//use n abscissa points, the roots of Legendre Polynomial P_n(x)

	const double lastbit = 0.5*(b - a);
	const double secondbit = 0.5*(a + b);

	std::vector<double> x = GaussLegendreTable();

	double answer = 0;
	for (int i = 0; i < n; i += 2) {
		double temp = x[i + 1] * f.evaluate((lastbit * (-1 * x[i])) + secondbit);
		if (temp == temp) { answer += temp; }//filter out singularities
		temp = x[i + 1] * f.evaluate((lastbit*x[i]) + secondbit);
		if (temp == temp) { answer += temp; }//filter out singularities
	}
	answer *= lastbit;
	x.clear();
	return answer;
}

double integrateGaussLegendreQuadrature(double a, double b, Function f, std::vector<double> vals, int pos) {//uses Chebyshev-Gauss quadrature rules -- slower but more precise
	//evaluates w.r.t. variable at position 'pos' in array
	const int n = 100;	//use n abscissa points, the roots of Legendre Polynomial P_n(x)

	const double lastbit = 0.5*(b - a);
	const double secondbit = 0.5*(a + b);

	std::vector<double> x = GaussLegendreTable();

	double answer = 0;
	for (int i = 0; i < n; i += 2) {
		double bit1 = (lastbit * (-1 * x[i])) + secondbit;
		vals.insert(vals.begin() + pos, bit1);//add value to vector of values to evaluate at
		double temp = x[i + 1] * f.evaluate(vals);
		vals.erase(vals.begin() + pos);//remove value
		if (temp == temp) { answer += temp; }//filter out singularities
		double bit2 = (lastbit*x[i]) + secondbit;
		vals.insert(vals.begin() + pos, bit2);//add value to vector of values to evaluate at
		temp = x[i + 1] * f.evaluate(vals);
		vals.erase(vals.begin() + pos);//remove value
		if (temp == temp) { answer += temp; }//filter out singularities
	}
	answer *= lastbit;
	x.clear();
	return answer;
}

double integrateGaussHermiteQuadrature(Function f) {//used to integrate bounds (-inf, +inf)
	const int n = 150;	//use n abscissa points, the roots of Hermite Polynomial H_n(x)
	std::vector<double> x = GaussHermiteTable();

	double answer = 0;
	for (int i = 0; i < n; i += 3) {
		double temp = x[i + 2] * f.evaluate(-x[i]);//sum w_i*f(-x_i)
		if (temp == temp) { answer += temp; }//filter out singularities
		temp = x[i + 2] * f.evaluate(x[i]);//sum w_i*f(x_i)
		if (temp == temp) { answer += temp; }//filter out singularities
	}
	x.clear();
	return (double)answer;
}

double integrateGaussLaguerreQuadrature(Function f) {//used to integrate functions where f(x) = g(x)*exp(-x) from bounds [0, +inf)
	const int n = 125 * 2;	//use n abscissa points, the roots of Laguerre Polynomial L_n(x)
	std::vector<long double> x = GaussLaguerreTable();
	double answer = 0;
	for (int i = 0; i < n; i += 2) {
		double temp = x[i + 1] * exp(x[i]) * f.evaluate(x[i]);//sum w_i*f(x_i), where w_i is a weight from the table
		if (temp == temp) { answer += temp; }//filter out singularities
	}
	x.clear();
	return answer;
}

double integrateGaussLaguerreQuadrature(Function f, std::vector<double> vals, int pos) {//used to integrate bounds [0, +inf) w.r.t. one variable
	const int n = 125 * 2;	//use n abscissa points, the roots of Polynomial L_n(x)
	std::vector<long double> x = GaussLaguerreTable();
	double answer = 0;
	for (int i = 0; i < n; i += 2) {
		vals.insert(vals.begin()+pos, x[i]);//add value to vector of values to evaluate at
		double temp = x[i + 1] * exp(x[i])* f.evaluate(vals);//sum w_i*f(x_i) where w_i is a weight from the table
		if (temp == temp) { answer += temp; }//filter out singularities
		vals.erase(vals.begin() + pos);//remove value
	}
	x.clear();
	return answer;
}

double integrateGaussKronodQuadrature(double a, double b, Function f) {//uses Gauss-Kronod quadrature rules to estimate error
	const int n = 98;	//use n abscissa points, the roots of Stieltjes Polynomial E_n(x)

	const double lastbit = 0.5*(b - a);
	const double secondbit = 0.5*(a + b);

	std::vector<double> x = GaussKronodTableCorrespondingGaussTable();
	std::vector<double> x_k = GaussKronodTable();
	
	double answer = 0;
	double KronodTerm = 0;
	for (int i = 0; i < n; i += 2) {
		//calculate Gauss terms
		double temp = x[i + 1] * f.evaluate((lastbit*x[i]) + secondbit);
		if (temp == temp) { answer += temp; }//filter out singularities

		//calculate Kronod terms
		temp = x_k[2*i + 1] * f.evaluate((lastbit*x_k[2*i]) + secondbit);
		if (temp == temp) { KronodTerm += temp; }//filter out singularities
		temp = x_k[2 * i + 3] * f.evaluate((lastbit*x_k[2 * i+2]) + secondbit);
		if (temp == temp) { KronodTerm += temp; }//filter out singularities

	}
	KronodTerm *= lastbit;
	answer *= lastbit;
	x.clear();
	x_k.clear();
	return answer + abs(answer - KronodTerm);
}

double integrateGaussKronodQuadrature(double a, double b, Function f, std::vector<double> vals, int pos) {
	//evaluates w.r.t. variable at position 'pos' in array
	const int n = 98;	//use n abscissa points, the roots of Legendre Polynomial P_n(x)

	const double lastbit = 0.5*(b - a);
	const double secondbit = 0.5*(a + b);

	std::vector<double> x = GaussKronodTableCorrespondingGaussTable();
	std::vector<double> x_k = GaussKronodTable();

	double answer = 0;
	double KronodTerm = 0;
	for (int i = 0; i < n; i += 2) {
		double bit2 = (lastbit*x[i]) + secondbit;
		vals.insert(vals.begin() + pos, bit2);//add value to vector of values to evaluate at
		double temp = x[i + 1] * f.evaluate(vals);
		vals.erase(vals.begin() + pos);//remove value
		if (temp == temp) { answer += temp; }//filter out singularities

		bit2 = (lastbit*x_k[2*i]) + secondbit;
		vals.insert(vals.begin() + pos, bit2);//add value to vector of values to evaluate at
		temp = x[2*i + 1] * f.evaluate(vals);
		vals.erase(vals.begin() + pos);//remove value
		if (temp == temp) { KronodTerm += temp; }//filter out singularities

		bit2 = (lastbit*x_k[2 * i+2]) + secondbit;
		vals.insert(vals.begin() + pos, bit2);//add value to vector of values to evaluate at
		temp = x[2 * i + 3] * f.evaluate(vals);
		vals.erase(vals.begin() + pos);//remove value
		if (temp == temp) { KronodTerm += temp; }//filter out singularities

	}
	answer *= lastbit;
	KronodTerm *= lastbit;
	x.clear();
	return answer+abs(answer-KronodTerm);
}

/*====================================
	VECTOR-VALUED FUNCTION ANALYSIS
====================================*/

std::vector<double> gradient(std::vector<double> X, Function f) {
	std::vector<double> answer;
	for (int i = 0; i < X.size(); ++i) {
		answer.push_back(nthPartialDerivative(1, i, X, f));
	}
	return answer;
}

std::vector<double> gradient(std::vector<double> X, VectorValuedFunction vf) {
	std::vector<double> answer;
	for (int i = 0; i < X.size(); ++i) {
		Function f = vf.f[i];
		answer.push_back(nthPartialDerivative(1, i, X, f));
	}
	return answer;
}

double divergence(std::vector<double> X, Function f) {
	std::vector<double> v = gradient(X, f);
	return arrsum(v);
}


double divergence(std::vector<double> X, VectorValuedFunction vf) {
	std::vector<double> v = gradient(X, vf);
	return arrsum(v);
}

std::vector<double> directionalDerivative(std::vector<double> X, Function f, std::vector<double> v) {
	std::vector<double> answer = gradient(X, f);
	if (arrNorm(v) != 1) { v = normedArray(v); }//make v a unit vector
	for (int i = 0; i < answer.size(); ++i) {
		answer[i] *= v[i];
	}
	return answer;
}

std::vector<double> directionalDerivative(std::vector<double> X, VectorValuedFunction vf, std::vector<double> v) {
	std::vector<double> answer = gradient(X, vf);
	if (arrNorm(v) != 1) { v = normedArray(v); }//make v a unit vector
	for (int i = 0; i < answer.size(); ++i) {
		answer[i] *= v[i];
	}
	return answer;
}

/*

double compose(double(*f)(double), double(*g)(double), double x) {
return f(g(x));         //composite function of 1D functions
}

double* composeNthDVectorValuedFunction(double*(*g)(double[]), double*(*f)(double[]), double X[]) {
return f(g(X));     //don't forget to declare f and g as functions that pass static arrays or else this will not work
}

double* composeNthDParametricVectorValuedFunction(double*(*g)(double[], double), double*(*f)(double[], double), double X[], double t) {
return f(g(X, t), t); //don't forget to declare f and g as functions that pass static arrays or else this will not work
}

double arcLength(double a, double b, double(*funct)(double))
{
double h = 0.00001;
double counter = a;
double answer = 0;
while (counter < b)
{
double dydx = (funct(counter + h) - funct(counter)) / h;
double lineLength = sqrt(1 + (dydx*dydx));
answer += lineLength;
counter += h;
}
answer *= h;
return answer;
}

double arcLengthMonteCarlo(double a, double b, double(*funct)(double))
{   //this method of arc length integral is suitable for wildly varying functions
double h = 0.000003;
double counter = a;
double answer = 0;
while (counter < b)
{
double randVal = boundedRand(counter, counter + h);
double dydx = (funct(randVal + h) - funct(randVal)) / h;
double lineLength = sqrt(1 + (dydx*dydx));
answer += lineLength;
counter += h;
}
answer *= h;
return answer;
}

double ODEFirstOrder(double t, double IVP, double(*funct)(double, double))
//formula: y_n+1 = y_n + h*f(t_n,y_n)
{   //we assume the IVP is at t=0 here
double partitionwidth = 0.0001; //determines how accurate the calculation is
double t_step = 0;                  //t0 = 0
double lastvalue = IVP;             //y0 = IVP
for (int p = partitionwidth; p < 9999; p++)
{
if (t_step == t)
{
return lastvalue;
}
double y_next = lastvalue + ((partitionwidth)*funct(lastvalue, t_step));    //this is our value for y_n+1
lastvalue = y_next;
t_step += partitionwidth;
}
return lastvalue;
}

double ODEFirstOrderImprovedEuler(double t, double IVP, double(*funct)(double, double))
//formula: y_n+1 = y_n + h*f(t_n,y_n)
{   //we assume the IVP is at t=0 here
double partitionwidth = 0.000001;   //determines how accurate the integral is
double t_step = 0;                  //t0 = 0
double lastvalue = IVP;             //y0 = IVP
while (t_step < t)
{
double y_next = lastvalue + ((partitionwidth)*funct(lastvalue, t_step));    //this is our value for y_n+1
y_next = lastvalue + (partitionwidth / 2)*((funct(lastvalue, t_step)) + ((funct(y_next, t_step + partitionwidth))));
lastvalue = y_next;
t_step += partitionwidth;
}
return lastvalue;
}


double ODEFirstOrderRK(double t, double IVP, double(*funct)(double, double))
//formula: y_n+1 = y_n + h/6(k1 + 2k2 + 2k3 + k4)
{   //we assume the IVP is at t=0 here
double h = 0.0001;  //determines how accurate the calculation is
double t_step = 0;                  //t0 = 0
double y_n = IVP;               //y0 = IVP  (aka y_n)
for (int p = h; p < 9999; p++)
{
if (t_step == t)
{
return y_n;
}
double k_one = funct(y_n, t_step);
double k_two = funct((y_n + ((h / 2)*k_one)), (t_step + (h / 2)));
double k_three = funct((y_n + ((h / 2)*k_two)), (t_step + (h / 2)));
double k_four = funct((y_n + (h*k_three)), (t_step + h));
double y_next = y_n + ((h / 6)*(k_one + (2 * k_two) + (2 * k_three) + k_four)); //this is y_n+1
y_n = y_next;
t_step += h;
}
return y_n;
}



double ODESecondOrder(double t, double IVP_one, double IVP_two, double(*funct)(double, double, double))
//formula:
{   //we assume the IVP is at t=0 here
double partitionwidth = 0.0000001;  //determines how accurate the calculation is
double t_step = 0;                          //t0  = 0
double lastvalue_y = IVP_one;               //y0  = IVP1
double lastvalue_v = IVP_two;           //y'0 or v0 = IVP2
while (t_step < t)
{
if (t_step == t)
{
return lastvalue_y;
}
double y_next = lastvalue_y + ((partitionwidth)*(lastvalue_v)); //v = y', so y_n+1 = y_n + h*v_n
double y_prime_next = lastvalue_v + (partitionwidth*funct(lastvalue_y, lastvalue_v, t_step));   //v': v_n+1=v_n + h*f(y_n,v_n,t_n)
lastvalue_v = y_prime_next;
lastvalue_y = y_next;
t_step += partitionwidth;
}
return lastvalue_y;         //this value can be replaced with lastvalue_v if you would like to know the derivative y' at t instead.
}

double ODESecondOrderRK(double t, double IVP_one, double IVP_two, double(*funct)(double, double, double))
//formula:
{   //we assume the IVP is at t=0 here
double h = 0.000001;    //determines how accurate the calculation is
double t_step = 0;                          //t0  = 0
double y_n = IVP_one;               //y0  = IVP1
double v_n = IVP_two;           //y'0 or v0 = IVP2
while (t_step < t)
{
if (t_step == t)
{
return y_n;
}
//solving y_n+1
double k_one = v_n;
double k_two = v_n + ((h / 2)*k_one);
double k_three = v_n + ((h / 2)*k_two);
double k_four = v_n + (h*k_three);
double y_next = y_n + ((h / 6)*(k_one + (2 * k_two) + (2 * k_three) + k_four));

//solving v_n+1
k_one = funct(y_n, v_n, t_step);
k_two = funct((y_n + ((h / 2)*k_one)), (v_n + ((h / 2)*k_one)), (t_step + (h / 2)));
k_three = funct((y_n + ((h / 2)*k_two)), (v_n + ((h / 2)*k_two)), (t_step + (h / 2)));
k_four = funct((y_n + (h*k_three)), (v_n + (h*k_three)), (t_step + h));
double y_prime_next = v_n + ((h / 6)*(k_one + (2 * k_two) + (2 * k_three) + k_four));
v_n = y_prime_next;
y_n = y_next;
t_step += h;
}
return y_n;         //this value can be replaced with v_n if you would like to know the derivative y' at t instead.
}


double convolve(double(*funct1)(double), double(*funct2)(double), double t){
double time_step = t / 900;     //900 is an arbitrary number of partitions
double h = 0;                   //h is the time step counter
double tau_step = 0.00001;      //tau_step is also arbitrary in size
double tau = 0;                 //tau is tau increment counter
double answer = 0;
int counter = 0;
while (counter < 900)
{
tau = h;
while (tau < h + time_step)
{
answer += tau_step * funct1(tau) * funct2(t - tau);
tau += tau_step;
}
h += time_step;
counter++;
}
return answer;
}
*/