// recursivequantile.cpp
// Normal quantile using 
// Steinbrecher's non-linear recursion
// This test program works on 0.004<u<0.996
// with max relative error < 3e-15
// using Bloodshed Dev C++
// (1.1e-15 with double -> long double everywhere)
// William Shaw
// Version 0.2 - provisional version
// for education and benchmarking only
// February 2007 
// william.shaw@kcl.ac.uk

#include <cmath>
#include <iostream>
#include <fstream>
#include "nrutil_nr.h" //public domain Numerical Recipes vectors
using namespace std;

double RecursiveQuantile(double u, NRVec<double> &rat)
{   int P=rat.size();
    int i;
    double b,s,t,q;
    const double rootpi = 1.77245385090551602729816748334;
    const double roottwo = 1.41421356237309504880168872421;
    b = rootpi*(u-0.5);
    i=-1;
    t=0;
    s=b;
    q=b*b;
    while(fabs(b)>1.0e-17 && i<= P) s=(t=s)+(b*=q*rat[i+=1]);
    if (i==P) 
    {cout << "recursion limit reached \n";
     cout << "switch to tail series \n ";
     cout << "last term = " << b  << "\n";}
    return roottwo*s;
}


int main()
{
    NRVec<double> qarr(2701);
    NRVec<double> ratios(2700);
    double u;
    double quantile;
    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] = Sum[q[m]*q[k - 1 - m]/(m + 1)/(2m + 1), {m, 0, k - 1}]
    qarr[0] = 1.0;
    qarr[1] = 1.0;
    for (k=2; k <= 2700; k++)
    {qarr[k]=0.0;
     for (m=0; m<k; m++)
    {qarr[k]+=qarr[m]*qarr[k-1-m]/(m+1)/(2*m+1);}
    };
    
    for (k = 1; k <= 2700; k++)
    {qarr[k]=qarr[k]/(2*k+1);
    ratios[k-1] = qarr[k]/qarr[k-1];}    

    cout << "Series ratios computed to order 2700 \n";
    cout << qarr[2700] << "\n";
    cout << ratios[2699]<< "\n";
    cout << "Working on test array \n";
    ofstream out("quantiles.txt");
    for (k=4; k<=996; k++)
    {u = k/1000.; 
    quantile = RecursiveQuantile(u, ratios);
    out.precision(25);
    out << u <<","<< quantile << "\n";}
    cout << "Output written to quantiles.txt \n";
    cout << "Hit any key to quit \n";
    cin >> name;
    return(0);
}