// Normal quantile: Steinbrecher's recursion
// This version works on 0.0005<u<0.9995 with max abs relative error 1.7< E-15,
// and on 0.0035 < u < 0.9965 with max abs relative error < 2.5 E-16, 
// using Bloodshed Dev C++, Win XP, Pentium M 1.7GHz
// Results on other hardware, compilers WILL vary, esp if poor long double.
// William Shaw
// Version 0.6 - provisional version for education & precision benchmarking only
// April 2007 
// william.shaw@kcl.ac.uk

#include <cmath>
#include <iostream>
#include <fstream>
using namespace std;

long double RecursiveQuantile(long double u, long double rat[19999])
{   int P=19999;
    int i;
    long double b,s,t,q;
    const long double rootpiotwo = 1.25331413731550025120788264241;
    b =2*u-1.0;
    i=-1;
    t=0;
    s=b;
    q=b*b;
    while((fabs(b)>1.0e-18 && i<= P-2) || i<5) s=(t=s)+(b*=q*rat[i+=1]);
    if (i==P-1) 
    {cout << "recursion limit reached \n";
     cout << "switch to tail series \n ";}
    return rootpiotwo*s;
}


int main()
{   long double qarr[20000];
    long double ratios[19999];
    long double u;
    long double quantile;
    const long double pioverfour = 0.785398163397448309615660845820;
    int maxsize = 20000;

    char name[5];
    int k,m;
    cout << "Steinbrecher's recursive quantile \n ";
    cout << "Initializing series..... \n";

//  Initialize the array of coefficients and ratios
//  using Steinbrecher recursion:
//  q[k] = (Pi/4) Sum[q[m]*q[k - 1 - m]/(m + 1)/(2m + 1), {m, 0, k - 1}]
    qarr[0] = 1.;
    for (k=1; k <= maxsize-1; k++)
    {qarr[k]=0.0;
     for (m=0; m<k/2; m++)
     {qarr[k]+=qarr[m]*qarr[k-1-m]*(1.0/(m+1)/(2*m+1)+1.0/(k-m)/(2*k-2*m-1));}
    if ((k-1)%2==0) {qarr[k]+= 2*qarr[(k-1)/2]*qarr[(k-1)/2]/(k+1)/k;}
    qarr[k]*=pioverfour;};
 
//  Now work out what to feed to power series
    for (k = 1; k <= maxsize-1; k++)
    {qarr[k]=qarr[k]/(2*k+1);
    ratios[k-1] = qarr[k]/qarr[k-1];}    

    cout << "Series ratios computed to order " << maxsize-1 << "\n";
    ofstream out("quantiles.txt");
    for (k=5; k<=9995; k++)
    {u = k/10000.; 
    quantile = RecursiveQuantile(u, ratios);
    out.precision(25);
    out << quantile << "\n";}
    cout << "Output written to quantiles.txt \n";
    cout << "Hit any key to quit \n";
    cin >> name;
    return(0);
}