#include <float.h>
#include <math.h>

void inverseerrorfunction(float x, float oneminx, float *inverf)
{
	float chepolsum(int, float, float []);
	float absx,p,betax,a[24];

	absx=fabs(x);
	if (absx > 0.8 && oneminx > 0.2) oneminx=0.0;
	if (absx <= 0.8) {
		a[0]  = 0.992885376618941;  a[1]  = 0.120467516143104;
		a[2]  = 0.016078199342100;  a[3]  = 0.002686704437162;
		a[4]  = 0.000499634730236;  a[5]  = 0.000098898218599;
		a[6]  = 0.000020391812764;  a[7]  = 0.000004327271618;
		a[8]  = 0.000000938081413;  a[9]  = 0.000000206734720;
		a[10] = 0.000000046159699;  a[11] = 0.000000010416680;
		a[12] = 0.000000002371501;  a[13] = 0.000000000543928;
		a[14] = 0.000000000125549;  a[15] = 0.000000000029138;
		a[16] = 0.000000000006795;  a[17] = 0.000000000001591;
		a[18] = 0.000000000000374;  a[19] = 0.000000000000088;
		a[20] = 0.000000000000021;  a[21] = 0.000000000000005;
		*inverf = chepolsum(21,x*x/0.32-1.0,a)*x;
	} else if (oneminx >= 25.0e-4) {
		a[0]  =  0.912158803417554;  a[1]  = -0.016266281867664;
		a[2]  =  0.000433556472949;  a[3]  =  0.000214438570074;
		a[4]  =  0.000002625751076;  a[5]  = -0.000003021091050;
		a[6]  = -0.000000012406062;  a[7]  =  0.000000062406609;
		a[8]  = -0.000000000540125;  a[9]  = -0.000000001423208;
		a[10] =  0.000000000034384;  a[11] =  0.000000000033584;
		a[12] = -0.000000000001458;  a[13] = -0.000000000000810;
		a[14] =  0.000000000000053;  a[15] =  0.000000000000020;
		betax=sqrt(-log((1.0+absx)*oneminx));
		p = -1.54881304237326*betax+2.56549012314782;
		p=chepolsum(15,p,a);
		*inverf = (x < 0.0) ? -betax*p : betax*p;
	} else if (oneminx >= 5.0e-16) {
		a[0]  =  0.956679709020493;  a[1]  = -0.023107004309065;
		a[2]  = -0.004374236097508;  a[3]  = -0.000576503422651;
		a[4]  = -0.000010961022307;  a[5]  =  0.000025108547025;
		a[6]  =  0.000010562336068;  a[7]  =  0.000002754412330;
		a[8]  =  0.000000432484498;  a[9]  = -0.000000020530337;
		a[10] = -0.000000043891537;  a[11] = -0.000000017684010;
		a[12] = -0.000000003991289;  a[13] = -0.000000000186932;
		a[14] =  0.000000000272923;  a[15] =  0.000000000132817;
		a[16] =  0.000000000031834;  a[17] =  0.000000000001670;
		a[18] = -0.000000000002036;  a[19] = -0.000000000000965;
		a[20] = -0.000000000000220;  a[21] = -0.000000000000010;
		a[22] =  0.000000000000014;  a[23] =  0.000000000000006;
		betax=sqrt(-log((1.0+absx)*oneminx));
		p = -0.559457631329832*betax+2.28791571626336;
		p=chepolsum(23,p,a);
		*inverf = (x < 0.0) ? -betax*p : betax*p;
	} else if (oneminx >= FLT_MIN) {
		a[0]  =  0.988575064066189;  a[1]  =  0.010857705184599;
		a[2]  = -0.001751165102763;  a[3]  =  0.000021196993207;
		a[4]  =  0.000015664871404;  a[5]  = -0.000000519041687;
		a[6]  = -0.000000037135790;  a[7]  =  0.000000001217431;
		a[8]  = -0.000000000176812;  a[9]  = -0.000000000011937;
		a[10] =  0.000000000000380;  a[11] = -0.000000000000066;
		a[12] = -0.000000000000009;
		betax=sqrt(-log((1.0+absx)*oneminx));
		p = -9.19999235883015/sqrt(betax)+2.79499082012460;
		p=chepolsum(12,p,a);
		*inverf = (x < 0.0) ? -betax*p : betax*p;
	} else
		*inverf = (x > 0.0) ? 26.0 : -26.0;
}