/* Solution of the S-d quantile equation for the case g<1.  
   Last modified on 10/11/01 by Sergej V. Aksenov. */

/* Return codes: 0 - ok, 1 - 6 lerchphi() codes. */

#include <math.h>

int sdqnt(double *sdpar, double *f0, double *f, double *x, double *prec)
{
  double alpha, g, h, x0, lambda, gamma, z1g, z1l, z2g, z2l, s, v, lerch1, lerch2, frac;
  int status, status1, status2, iter1, iter2;

  /* Set up parameters. */

  g = sdpar[0];
  h = sdpar[1];
  alpha = sdpar[2];
  x0 = sdpar[3];
  
  lambda = 1.0 - g;
  gamma = h - g;
  frac = lambda/gamma;
  v = 1.0 + frac;
  s = 1.0;
  
  if (*f == 0.0) z2g = z2l = 0.0;
        else
	  {
	    z2g = pow(*f,gamma);
	    z2l = pow(*f,lambda);
	  }
        
  if (*f0 == 0.0) z1g = z1l = 0.0;
        else
	  {
	    z1g = pow(*f0,gamma);
	    z1l = pow(*f0,lambda);
	  }

  /* Call lerchphi() function to compute the Lerch Phi. */

  status1 = lerchphi(&z1g, &s, &v, prec, &lerch1, &iter1);
  status2 = lerchphi(&z2g, &s, &v, prec, &lerch2, &iter2);
     
  *x = x0 + (z2l*(1.0+frac*z2g*lerch2) - z1l*(1.0+frac*z1g*lerch1)) /alpha/lambda;

  if (status1) status = status1;
  else status = 0;
  
  if (status2) status = status2;
  return status;
  
}