//	file test.S uses this
//	gcc -m64 -c test.c
// Note that location changed to -location in tailslogMSS (two changes)

#include <math.h>
#include <stdio.h>

#define ERRORLOG "errorlog.txt"

#define MIN(a,b) ((a)<(b)? (a):(b))
#define MAX(a,b) ((a)<(b)? (b):(a))

#define OI 16
#define HOI 8
#define OIM1 15
/*      OI is the order of interpolation to be used
HOI is half the order of interpolation to be used
OIM1 is OI-1    */

#define FALSE 0
#define TRUE 1

/*=====================================================================*/
double expoffset(double x)
{
/* exp(x)-1 using continued fraction 4.2.40, page 70 of Abramowitz and Stegun */
/* Number of terms has been checked as just adequate for 8-byte precision */

double xsq,val;
if(fabs(x) > 0.1) val=exp(x)-1.;
else {
	xsq=x*x;
	val=x/(1.-x/2.+xsq/(12.*(1.+xsq/(60.*(1.+xsq/(140.*(1.+xsq/(252.))))))));
}
return val;
}
/*=====================================================================*/
double logoffset(double x)
{
/* log(1+x) using fraction 4.1.40, page 68 of Abramowitz and Stegun */
/*	with (1+z)/(1-z)=1+x, so that z=x/(2+x). */
/* Number of terms has been checked as adequate for 8-byte precision  */

double val,z,zsq;
if(fabs(x) > 0.1) val=log(1.+x);
else {
	z=x/(2.+x);
	zsq=z*z;
       	val=2.*z/(1.-zsq/(3.-4.*zsq/(5.-9.*zsq/(7.-16.*zsq/(9.-25.*zsq/(11.))))));
}
return val;
}
/*====================================================================== */
double millsratio(double x)
{
/* Computes Mill's ratio, which is the ratio of the right tail */
/*	probability to the density for a standard normal distribution. */
/* Note that this function should not be called for negative arguments. */

FILE * out;
double val,t;
double rhodiff[26]={1.,3.94766755064487339,
	.262699373075432297,-4.79530994886257864,-.472387694632788044,
	-4.64178612251027508,.888029797897972282E-1,3.34506626993363627,
	-.253018663557510120,-3.22271919537673200,.367907414060725404,
	.515135649489030960,-.782205852853700406,-1.66658903030166903,
	-.317783299853388699,3.44266858016113165,.110601063267646032,
	-113.675951225046565,-.475674678528123486E-2,891.504346794291090,
	.155839071476740257E-3,-52392.5264048076612,-.145796580254323180E-4,
	5331522.33187226848,.169353718097630593E-6,-83341339.7482781260};
double xdata[25]={0.,1.,.2704,.5776,.0784,.04,.0576,.0256,.1024,
	.1296,.16,.1936,.2304,.0064,.3136,.36,.4096,.4624,
	.5184,.0144,.64,.7056,.7744,.8464,.9216};

if(x<0.){
	out=fopen(ERRORLOG,"a");
	fprintf(out,"(millsratio) Argument x should not be negative");
	fclose(out);
	return 0.;
}

/* Thiele interpolation by Abramowitz and Stegun 25.2.50*/
t=1./(x+1.);
val=t*(rhodiff[0]+(t-xdata[0])/(rhodiff[1]+(t-xdata[1])/(rhodiff[2]+
	(t-xdata[2])/(rhodiff[3]+(t-xdata[3])/(rhodiff[4]+
	(t-xdata[4])/(rhodiff[5]+(t-xdata[5])/(rhodiff[6]+
	(t-xdata[6])/(rhodiff[7]+(t-xdata[7])/(rhodiff[8]+
	(t-xdata[8])/(rhodiff[9]+(t-xdata[9])/(rhodiff[10]+
	(t-xdata[10])/(rhodiff[11]+(t-xdata[11])/(rhodiff[12]+
	(t-xdata[12])/(rhodiff[13]+(t-xdata[13])/(rhodiff[14]+
	(t-xdata[14])/(rhodiff[15]+(t-xdata[15])/(rhodiff[16]+
	(t-xdata[16])/(rhodiff[17]+(t-xdata[17])/(rhodiff[18]+
	(t-xdata[18])/(rhodiff[19]+(t-xdata[19])/(rhodiff[20]+
	(t-xdata[20])/(rhodiff[21]+(t-xdata[21])/(rhodiff[22]+
	(t-xdata[22])/(rhodiff[23]+(t-xdata[23])/(rhodiff[24]+
	(t-xdata[24])/(rhodiff[25]))))))))))))))))))))))))));
return val;
}
/*=====================================================================*/
double lognormaltail(double z)
{
/* Finds the logarithm of right hand tail probability for a standard normal distribution. */
/* Underflow will affect the result if the argument z is large negative. */

double val;
const double dnorm0=.39894228040143267793994605993438187;
if(z<0.) val=logoffset(-dnorm0*exp(-.5*z*z)*millsratio(-z));
else val=log(dnorm0*millsratio(z))-.5*z*z;
return val;
}
/*=====================================================================*/

double normaltail(double z)
{
/* Finds the right hand tail probability for a standard normal distribution */
/* Note that precision is poor if the argument z is negative. */
/* Will return zero if z is large. */

double val;
double const dnorm0=.39894228040143267793994605993438187;
if(z<0.) val=1.-dnorm0*exp(-.5*z*z)*millsratio(-z);
else val=dnorm0*exp(-.5*z*z)*millsratio(z);
return val;
}

/*=====================================================================*/

double LogGamma(double x)
{
/*  Abramowitz & Stegun equation 6.1.48	 */

double Z,val,c0=0.9189385332046727417803297364056177,
	a0=1/12.,a1=1/30.,a2=53/210.,a3=195/371.,a4=22999/22737.,
	a5=29944523./19733142.,a6=109535241009./48264275462.;
if(x>10) val=(x-.5)*log(x)+c0-x+a0/(x+a1/(x+a2/(x+a3/(x+a4/(x+a5/(x+a6/x))))));
else {
	Z=x+9;
	val=(Z-.5)*log(Z)+c0-Z+a0/(Z+a1/(Z+a2/(Z+a3/(Z+a4/(Z+a5/(Z+a6/Z
		))))))-log(x*(x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6)*(x+7)*(x+8));
}
return val;
}
/*=====================================================================*/
double expectationlogMSS(double alpha,double location,double logscale,double power)
/* Calculates expectation of a specified power of a log MSS distribution */
// This should not give exceptions for the values likely to be input..???
{
static const double hpi=1.57079632679489661923132169163975144209858;
double res,ps,eps,s,logps;

/* When alpha < 0.5: log MSS variable is
	exp{location-exp(logscale)*(parametrization C standard)} */
/* When alpha >=0.5: log MSS variable is
	exp{location-exp(logscale)*(parametrization M=S0 standard)} */

logps=log(power)+logscale;
if(alpha<0.5)
	res=exp(power*location-exp(logps*alpha));
else {
	ps=exp(logps);
	if(alpha == 1)
		res=exp(power*location+logps*ps/hpi);
	else if(alpha == 2)
		res=exp(power*location+ps*ps);
	else{
		/* Simpler formula; less accurate if alpha is very near to unity
		res=exp(power*location-pow(ps,alpha)/cos(hpi*alpha)+ps*tan(hpi*alpha));
		*/
		eps=alpha-1;
		s=sin(.5*hpi*eps);
		res=exp(power*location+ps*(expoffset(eps*logps)+2*s*s)/sin(hpi*eps));
	}
}
return res;
}
/*==================================================================== */
void iidcombine(double alpha,double location,double logscale,
	double number,double *location_product,double *logscale_product)
{
/* Given parameters of a log maximally skew stable distribution,
	find the location and logscale for the product of a specified
	number of such distributions. (The number may be less than unity.) */
double calpha,scale;
static const double hpi=1.57079632679489661923132169163975144209858;

*logscale_product = logscale +log(number)/alpha;
/* For A=S1 parametrization, location parameters add and scales^(alpha) add */
if(alpha<.5){
		*location_product=number*location;
}
else{
	calpha=1-alpha;
	scale=exp(logscale);
	if (calpha == 0)*location_product=number*location-
		scale*number*log(number)/hpi;
	else if (fabs(calpha) < 0.1)*location_product=number*location-
		scale*number*expoffset(calpha*log(number)/alpha)/tan(hpi*calpha);
	else *location_product=number*location-
		scale*tan(hpi*alpha)*(pow(number,1/alpha)-number);
}
}
/*==================================================================== */
void match2logMSS(double mean,double standard_deviation,double alpha,double *location,double *logscale)
{
/* Given alpha, find parameters of a log maximally skew stable distribution */
/*	to fit given mean and standard deviation */
double logmr,s,alphaminus1,scale_to_alpha,sinalpha,cv,scale;
static const double log2=.6931471805599453094;
static const double pi=3.141592653589793238462643383279502884197;
static const double hpi=1.57079632679489661923132169163975144209858;
FILE *out;

/* When alpha < 0.5:
	MSS variable is exp(logscale)*(parametrization C standard)-location
	log MSS variable is exp{location-exp(logscale)*(parametrization C standard)} */
/* When alpha >=0.5:
	MSS variable is exp(logscale)*(parametrization M=S0 standard)-location
	log MSS variable is exp{location-exp(logscale)*(parametrization M=S0 standard)} */

if(mean <= 0 | standard_deviation <= 0){
	out=fopen(ERRORLOG,"a");
	fprintf(out,"(match2logMSS) Mean and standard deviation must be strictly positive\n");
	fprintf(out,"(match2logMSS) Mean and standard deviation given were %f %f",mean, standard_deviation);
	fclose(out);
	return;
}

cv=standard_deviation/mean;
logmr=logoffset(cv*cv);
/*printf("logmr= %f\n",logmr); */
alphaminus1=alpha-1;
/*printf("\nalpha= %f\n",alpha); */

/* Case when alpha is less than 0.5 */
if(alpha<0.5){
	scale_to_alpha=logmr/(2-pow(2,alpha));
/*	printf("scale_to_alpha= %f\n",scale_to_alpha); */
	*logscale=log(scale_to_alpha)/alpha;
/*	printf("low: logscale= %f\n",*logscale); */
	*location=log(mean)+scale_to_alpha;
/*	printf("location= %f\n",*location); */
/* Check by printing calculated moments */
/*	printf("Calculated first moment: %f\n",exp(*location-scale_to_alpha));*/
/*	printf("Calculated second moment: %f\n",exp(2* *location-pow(2,alpha)*scale_to_alpha)); */
}

/* Case when alpha is precisely 1 */
else if(alphaminus1 == 0){
	scale=pi*logmr/(4*log2);
/*	printf("unity: scale= %f\n",scale); */
	*logscale=log(scale);
/*	printf("unity: logscale= %f\n",*logscale); */
	*location=log(mean)-*logscale*2*(scale)/pi;
/*	printf("location= %f\n",*location); */
/* Check by printing calculated moments */
/*	printf("Calculated first moment: %f\n",exp(*location)*pow(scale,(2*(scale)/pi))); */
/*	printf("Calculated second moment: %f\n",exp(2**location)*pow((2*(scale)),(4*(scale)/pi))); */
}

/* Case when alpha exceeds 0.5 but is not precisely 1 */
else{
	s=sin(alphaminus1*hpi);
	sinalpha=sin(alpha*hpi);
/*	printf("s= %f\n",s); */
	scale_to_alpha=.5*logmr*s/expoffset(alphaminus1*log2);
/*	printf("scale_to_alpha= %f\n",scale_to_alpha); */
	*logscale=log(scale_to_alpha)/alpha;
/*	printf("high: logscale= %f\n",*logscale); */
	*location=log(mean)-(scale_to_alpha - exp(*logscale)*sinalpha)/s;
/*	printf("location= %f\n",*location); */

/* Check by printing calculated moments */
	scale=pow(scale_to_alpha,(1/alpha));
/*	printf("high: scale= %f\n",scale); */
/*	printf("Calculated first moment: %f\n",exp(*location+scale_to_alpha/s-(scale)*sinalpha/s)); */
/*	printf("Calculated second moment: %f\n",exp(2**location+pow(2,alpha)*scale_to_alpha/s-2*(scale)*sinalpha/s)); */
}
}
/*=====================================================================*/
double mid3(double x,double xa,double ya,double y1a,double y2a,
	double xb,double yb,double y1b,double y2b)
{
/* Find interpolated value for y at x given values of y and its first
two derivatives at a and at b. */

double h,t,h2,a0,a1,a2,a3,a4,a5,r0,r1,r2;

//Change x-axis so that a=0 and b=1
h=xb-xa;
y1a=y1a*h;
y1b=y1b*h;

h2=h*h;
y2a=y2a*h2;
y2b=y2b*h2;

a0=ya;
a1=y1a;
a2=.5*y2a;
r0=yb-ya-a1-a2;
r1=y1b-a1-2*a2;
r2=y2b-2*a2;
a5=.5*r2-3*r1+6*r0;
a4=r1-3*r0-2*a5;
a3=yb-a0-a1-a2-a4-a5;

t=(x-xa)/h;
return a0+t*(a1+t*(a2+t*(a3+t*(a4+t*a5))));
}

/* ============================================================= */
/* -----------------          zbrentf           ------------------*/
/* ====      Based on functions from Numerical Recipes      ==== */
/* ====      by Press, Flannery, Teukolsy & Vetterling      ==== */
/* ============================================================= */

double zbrentf(double (*func)(double),double xa,double xb,
	double fa,double fb,double tol,int *errorcode)
/* To find a zero of a univariate function func.
A maximum of 100 iterations is set.
On entry, x values xa and xb must bracket a zero.
Function values must be passed as fa and fb.
Errorcode will be zero on succesful return from zbrentf.*/

{
    int iter;
    double a,b,c,d,e,min1,min2;
    double fc,p,q,r,s,tol1,xm;
    FILE * out;

    *errorcode = 0;
    if (fb*fa > 0.){
	*errorcode = -1; /* Root must be bracketed */
	out=fopen(ERRORLOG,"a");
	fprintf(out,"(zbrentf) Root not bracketed by values at xa and xb");
	fclose(out);
	}
    else {
	a=xa;
	b=xb;
	fc=fb;
	for (iter=1;iter<=100;iter++) {
	    if (fb*fc > 0.) {
		c = a;
		fc = fa;
		e = d = b-a;
	    }
	    if (fabs(fc) < fabs(fb)) {
		a = b;
		b = c;
		c = a;
		fa = fb;
		fb = fc;
		fc = fa;
	    }
	    tol1 = 1.e-15*fabs(b)+0.5*tol;
	    xm = 0.5*(c-b);
	    /************/
	    /* printf("Iteration %d: Current estimate is %.15g\n",iter + 1,b);*/
	    /************/
	    if (fabs(xm) <= tol1 || fb == 0.) break; /* Done */
	    if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
		s = fb/fa;
		if (a == c) {
		    p = 2.0*xm*s;
		    q = 1.0-s;
		}
		else {
		    q = fa/fc;
		    r = fb/fc;
		    p = s*(2.*xm*q*(q-r)-(b-a)*(r-1.));
		    q = (q-1.)*(r-1.)*(s-1.);
		}
		if (p > 0.)  q = -q;
		p = fabs(p);
		min1 = 3.*xm*q-fabs(tol1*q);
		min2 = fabs(e*q);
		if (2.*p < (min1 < min2 ? min1 : min2)) {
		    e = d;
		    d = p/q;
		}
		else {
		    d = xm;
		    e = d;
		}
	    }
	    else {
		d = xm;
		e = d;
	    }
	    a = b;
	    fa = fb;
	    if (fabs(d) > tol1) b += d;
	    else b += (xm > 0. ? fabs(tol1) : -fabs(tol1));
	    fb = (*func)(b);
	}
	if (iter == 100){
	    *errorcode = -2;
	    out=fopen(ERRORLOG,"a");
	    fprintf(out,"(zbrentf) Number of iterations reached maximum (100)");
	    fclose(out);
	    }
    }
    return (b);
}
/*====================================================================== */
double approx_quantile(double v,double xa,double va,double da,
	double xb,double vb,double db,double *estimated_precision)
{
/* Find approximate quantile of a function, given that it has value va
and derivative da at a and value vb and derivative db at b. */

double t2,t3,val,h,t2pt3,t22pt33,Q,rootQ,tp1,tp2,first;

//Change y-axis so that require a zero.
va=va-v;
vb=vb-v;

//Check that function values are one each side of the required value.
if(va*vb >0){printf("Root apparently not within range\n");}

//Change x-axis for that a=0 and b=1
h=xb-xa;
//Check that a and b are different.
if(fabs(h) == 0){printf("Ends of range are identical???/n");return 0.;}
da=da*h;
db=db*h;

//Find quadratic and cubic terms.
t2pt3=vb-va-da;
t22pt33=db-da;
t3=t22pt33-2*t2pt3;
t2=t2pt3-t3;

//Check that fitted derivative does not have a root in the range.
Q=t2*t2-4*da*t3;
if(Q>0){
	rootQ=sqrt(Q);
	tp1=(rootQ-t2)/(2*t3);
	tp2=(-rootQ-t2)/(2*t3);
	if((tp1 >0 & tp1 <1) |(tp2 >0 & tp2 <1)){
		printf("No good\n");
		return 0.;
	}
}

//Do three Newton-Raphson iterations on the cubic approximation.
if(fabs(va)<fabs(vb))val=-va/da;
	else val=1-vb/db;
first=val;
val=val-(va+val*(da+val*(t2+val*t3)))/(da+val*(2*t2+val*3*t3));
val=val-(va+val*(da+val*(t2+val*t3)))/(da+val*(2*t2+val*3*t3));
*estimated_precision=fabs(val-first);

return xa+h*val;
}
/*====================================================================== */
double x4coeff(double f0,double d0,double h,double fh,double dh,double a,double fa)
/* Compute the coefficient of x^4 in a polynomial through
	f(0)=f0,	f'(0)=d0
	f(h)=fh,	f'(h)=dh
	f(a)=fa			*/
{
double rh,ra,D,E,F,r;

rh=1/h;
D=rh*(rh*(fh-f0)-d0);
E=rh*(dh-d0);
ra=1/a;
F=ra*(ra*(fa-f0)-d0);
r=1/(a-h);
return ((F-D)*r+(2*D-E)*rh)*r;
}
/*====================================================================== */
void calc_recip_denom(int nx,double x[],double denom[nx-OIM1][OI])
/* Calculates reciprocals of denominators for use in later interpolation */
{
double product;
int i,j,k,start,offset;
for (i=0; i<nx-OIM1; i++){
	offset=i;
	for (j=0; j<OI; j++){
		product=1;
		for (k=0; k<OI; k++){
			if(k != j)product=product*(x[j+offset]-x[k+offset]);
		}
		denom[i][j]=1./product;
	}
}
}
/*========================================================================= */
void interpolate(double x,double *f,double *d,int nxn,double xn[nxn],
	double fn[nxn],double dn[nxn],double xdenomn[nxn-OIM1][OI])
{
/* Find interpolated values for x using vectors xn, fn and dn with */
/* reciprocals of denominators specified by array xdenomn */
double difference[OI],product,weight;
int low,high,mid,offset,start,k;
FILE *out;

/* First find the smallest index of a larger value of xn by binary chop */
low=0;
high=nxn-1;
if(xn[high]<x){
	out=fopen(ERRORLOG,"a");
	fprintf(out,"(interpolate) x outside range of data\n");
	fclose(out);
	printf("Problem in subroutine interpolate\n");
	printf("x,nxn %f %d\n",x,nxn);
	printf("first and last xn: %f %f\n",xn[0],xn[nxn-1]);
	return;
}
do{
	mid=(low+high)/2;
/* printf("low,mid,high= %f %f %f\n",low,mid,high); */
	if(xn[mid] >= x) high=mid;
	else{
		if(low == mid)break;
		low=mid;
	}
} while(TRUE);

/*printf("high= %d\n",high); */
start=MIN(MAX(0,high-HOI),nxn-OI);
offset=start;
/*printf("start= %d\n",start); */
product=1;
for (k=0; k<OI; k++){
	difference[k]=x-xn[k+offset];
	product=product*difference[k];
/*	printf("difference[%d]= %g\n",k,difference[k]); */
}
/*printf("product= %g\n",product); */
if(product == 0){
	for (k=0; k<OI; k++){
/*Use appropriate alpha value */
		if(x == xn[k+offset]){
			*f=fn[k+offset];
			*d=dn[k+offset];
			break;
		}
	}
}
else{
	*f=0;
	*d=0;
	for (k=0; k<OI; k++){
		weight=product*xdenomn[start][k]/difference[k];
		/*printf("k,weight= %d %g\n",k+offset,weight); */
		*f+=weight*fn[k+offset];
		/*printf("fn[%d]= %g\n",k+offset,fn[k+offset]); */
		/*printf("f= %g\n",*f); */
		*d+=weight*dn[k+offset];
		/*printf("d= %g\n",*d); */
		}
	}
}
/*=================================================================== */
void interpolate_over_alpha(int nx,int nalpha,double alphalist[nalpha],
	double thisalpha,double tablef[nx][nalpha],double tabled[nx][nalpha],
	double thisf[nx],double thisd[nx],double denom[nalpha-OIM1][OI])
{
/* To interpolate the tables tablef and tabled over alpha */
double weight,product,difference[OI];
int i,j,k,start,offset;

/* First find the smallest index of a larger value of alpha */
/* (This could be made faster with a binary chop algorithm) */
for (j=0; j<nalpha; j++){
	if(alphalist[j] > thisalpha)break;
}
start=MIN(MAX(0,j-HOI),nalpha-OI);
offset=start;
/* printf("start= %d\n",start); */
product=1;
for (k=0; k<OI; k++){
	difference[k]=thisalpha-alphalist[k+offset];
	product=product*difference[k];
}
/*printf("product:%g\n",product); */
/*printf("difference:%g\n",difference) */

/* Minor option: when thisalpha is a tabulated value */
if(product == 0){
	for (k=0; k<OI; k++){
		if(thisalpha == alphalist[k+offset]){
			for (i=0; i<nx; i++){
				thisf[i]=tablef[i][k+offset];
				thisd[i]=tabled[i][k+offset];
			}
			break;
		}
	}
}

/* Major option: need to interpolate */
else{
	for (i=0; i<nx; i++){
		thisf[i]=0;
		thisd[i]=0;
	}
	for (k=0; k< OI; k++){
		weight=product*denom[start][k]/difference[k];
/*printf("difference[%d] %g\n",k,difference[k]); */
/*printf("denom[%d][%d]: %g\n",start,k,denom[start][k]); */
/*printf("weight=  %g\n",weight); */
		for (i=0; i<nx; i++){
			thisf[i]+=weight*tablef[i][k+offset];
			thisd[i]+=weight*tabled[i][k+offset];
		}
/*printf("thisf %f\n",thisf) */
/*printf("thisd %f\n",thisd) */
	}
}
}
/*====================================================================== */
static const double pi=3.141592653589793238462643383279502884197;
static const double hpi=1.57079632679489661923132169163975144209858;
static const double EulersC=.5772156649015328606065120900824024310422;
static const double neglarge=-1.7E308;
static const double poslarge=1.7E308;
static const double logposlarge=709.7;

static double previous_alpha=-999.;

static double absam1,a1,sa2,nu,Calpha_C,Calpha_M,midpoint,xi,ximid,xlowlimit;
static double Clogd,alphastar,calpha,calphastar,ratio,angle,C1,C2,C3,C4,eta;
static double ffound,dfound,logapprox,logmidpoint,logscalef;

/* Boolean variables stored as integers. */
static int distributiontabulated, optionstabulated, alphasmall, alphalarge;

/* Interpolation table i gives values for density as di and for distribution */
/* function (or its complement) as fi for nxi values of x and nyi values of y*/

/* First tables are for alpha < 1, alpha*xi > 1/5 in C parametrization */
/* Vy1 is alpha, Vx1 is proportional to 1/(alpha*xi)	*/ 
#define nx1 70
#define ny1 20
static double Vx1[nx1],Vy1[ny1],f1[nx1],d1[nx1];
static double tablef1[nx1][ny1],tabled1[nx1][ny1];
static double xdenom1[nx1-OIM1][OI],ydenom1[ny1-OIM1][OI];

/* Second tables are for alpha < .5, alpha*xi < 1/5 and x<1 in C parametrization */
/* Vy2 is alpha, Vx2 is proportional to x**(-1/alpha) */
#define nx2 20
#define ny2 20
static double Vx2[nx2],Vy2[ny2],f2[nx2],d2[nx2];
static double tablef2[nx2][ny2],tabled2[nx2][ny2];
static double xdenom2[nx2-OIM1][OI],ydenom2[ny2-OIM1][OI];

/* Third tables are for alpha < .5, x > 1 in C parametrization */
/* Vy3 is alpha, Vx3 is proportional to x**(-1/alpha) */
#define nx3 20
#define ny3 20
static double Vx3[nx3],Vy3[ny3],f3[nx3],d3[nx3];
static double tablef3[nx3][ny3],tabled3[nx3][ny3];
static double xdenom3[nx3-OIM1][OI],ydenom3[ny3-OIM1][OI];

/* Fourth tables are for 1.7 <alpha < 2, -1.3 < x (M=S0) < 20 */
/* Use the difference from the alpha=2 (normal) distribution for all larger x, */
/* scaling so that the asymptotic behaviour as x tends to infinity is handled accurately. */
/* Vy4 is alpha, Vx4 is based on x in the S0 parametrization */
#define nx4 100
#define ny4 17
static double Vx4[nx4],Vy4[ny4],f4[nx4],d4[nx4];
static double tablef4[nx4][ny4],tabled4[nx4][ny4];
static double xdenom4[nx4-OIM1][OI],ydenom4[ny4-OIM1][OI];
static double f4_alpha2[nx4],d4_alpha2[nx4];

/* Fifth tables are for 1.7 <alpha < 2, x (M=S0) > 20 */
/* Use Zolotarev 2.5.6 */
/* Vy5 is alpha, Vx5 is proportional to y**(-1/alpha) where x=y+eta y**(1-alpha) */
#define nx5 20
#define ny5 17
static double Vx5[nx5],Vy5[ny5],f5[nx5],d5[nx5];
static double tablef5[nx5][ny5],tabled5[nx5][ny5];
static double xdenom5[nx5-OIM1][OI],ydenom5[ny5-OIM1][OI];

/* Sixth tables are for 0.5 <alpha < 1.7, x (M=S0) > 5 */
/* Use Zolotarev 2.5.6 */
/* Vy6 is alpha, Vx6 is proportional to y**(-1/alpha) where x=y+eta y**(1-alpha) */
#define nx6 20
#define ny6 40
static double Vx6[nx6],Vy6[ny6],f6[nx6],d6[nx6];
static double tablef6[nx6][ny6],tabled6[nx6][ny6];
static double xdenom6[nx6-OIM1][OI],ydenom6[ny6-OIM1][OI];

/* Seventh tables are for 0.5 <alpha < 1.7, from xi=2/5 to x=7.3 */
/* 	Table 1 goes from xi=.2/alpha to xi=infinity, so regions overlap. */
/* Vy7 is alpha, Vx7 is x (S0 parametrization) */
#define nx7 60
#define ny7 40
static double Vx7[nx7],Vy7[ny7],f7[nx7],d7[nx7];
static double tablef7[nx7][ny7],tabled7[nx7][ny7];
static double xdenom7[nx7-OIM1][OI],ydenom7[ny7-OIM1][OI];

/* Tabulation of x, density, distribution function and right tail */
int ntabs;
#define MAX_NTABS 200
static double tab_x[MAX_NTABS],tab_density[MAX_NTABS],tab_F[MAX_NTABS],tab_cF[MAX_NTABS],
	tab_lx[MAX_NTABS],tab_ld[MAX_NTABS],tab_lF[MAX_NTABS],tab_lcF[MAX_NTABS];
/* MAX_NTABS	Largest number of tabulated points for any alpha */
/* ntabs	Number of tabulated points for current alpha */
/* tab_x		x coordinates corresponding to tabulated values (ascending order) */
/* tab_lx	logarithms of tabx values (used for alpha < 0.5) */
/* tab_density		density of maximally-skew stable distribution */
/* tab_ld	logarithm of density of maximally-skew stable distribution */
/* tab_F		distribution function of maximally-skew stable distribution */
/* tab_lF	logarithm of distribution function of maximally-skew stable distribution */
/* tab_cF	right tail of maximally-skew stable distribution */
/* tab_lcF	logarithm of right tail of maximally-skew stable distribution */

/* F, cF, density, call option value for various strike prices */
int nstrikes;
static double previous_location, previous_logscale, previous_precision=1.e-5;
double meanprice,explogscale;
#define MAX_NSTRIKES 500
static double l_strike[MAX_NSTRIKES], l_call[MAX_NSTRIKES],
	l_cF[MAX_NSTRIKES], l_density[MAX_NSTRIKES];

/*========================================================================= */
static void setalpha(double alpha)
{
/* functions setalpha might be made faster by initializing arrays */
double weight,product,y,difference[OI],sinangle;
int i,j,k,start,offset;
FILE *out;
FILE *in;

/* If alpha is not in permissible range then print message but do nothing */
if(alpha>2. | alpha <= 0.){
	out=fopen(ERRORLOG,"a");
	fprintf(out,"(SETALPHA) Alpha not in range 0 < alpha <= 2");
	fclose(out);
	return;
}

/* If alpha is 2 or the same as the last time, then do nothing. */
if (alpha == 2 | alpha == previous_alpha)return;

/* On first call, read array margins as vectors and tabulated arrays */
if (previous_alpha == -999.){

/* First check that MAX_NTABS is large enough */
	if(MAX_NTABS < nx1+MAX(MAX(nx2+nx3,nx4+nx5),nx6+nx7)){
		out=fopen(ERRORLOG,"a");
		fprintf(out,"(SETALPHA) Need to increase parameter MAX_NTABS");
		fclose(out);
		return;
	}

	in=fopen("interpolationdata.keep","r");
	for (i=0; i<nx1; i++) fscanf(in,"%lf",&Vx1[i]);
	for (i=0; i<ny1; i++) fscanf(in,"%lf",&Vy1[i]);
	for (i=0; i<nx2; i++) fscanf(in,"%lf",&Vx2[i]);
	for (i=0; i<ny2; i++) fscanf(in,"%lf",&Vy2[i]);
	for (i=0; i<nx3; i++) fscanf(in,"%lf",&Vx3[i]);
	for (i=0; i<ny3; i++) fscanf(in,"%lf",&Vy3[i]);
	for (i=0; i<nx4; i++) fscanf(in,"%lf",&Vx4[i]);
	for (i=0; i<ny4; i++) fscanf(in,"%lf",&Vy4[i]);
	for (i=0; i<nx5; i++) fscanf(in,"%lf",&Vx5[i]);
	for (i=0; i<ny5; i++) fscanf(in,"%lf",&Vy5[i]);
	for (i=0; i<nx6; i++) fscanf(in,"%lf",&Vx6[i]);
	for (i=0; i<ny6; i++) fscanf(in,"%lf",&Vy6[i]);
	for (i=0; i<nx7; i++) fscanf(in,"%lf",&Vx7[i]);
	for (i=0; i<ny7; i++) fscanf(in,"%lf",&Vy7[i]);
	for (j=0;j<ny1;j++) for (i=0;i<nx1;i++) fscanf(in,"%lf",&tablef1[i][j]);
	for (j=0;j<ny1;j++) for (i=0;i<nx1;i++) fscanf(in,"%lf",&tabled1[i][j]);
	for (j=0;j<ny2;j++) for (i=0;i<nx2;i++) fscanf(in,"%lf",&tablef2[i][j]);
	for (j=0;j<ny2;j++) for (i=0;i<nx2;i++) fscanf(in,"%lf",&tabled2[i][j]);
	for (j=0;j<ny3;j++) for (i=0;i<nx3;i++) fscanf(in,"%lf",&tablef3[i][j]);
	for (j=0;j<ny3;j++) for (i=0;i<nx3;i++) fscanf(in,"%lf",&tabled3[i][j]);
	for (j=0;j<ny4;j++) for (i=0;i<nx4;i++) fscanf(in,"%lf",&tablef4[i][j]);
	for (j=0;j<ny4;j++) for (i=0;i<nx4;i++) fscanf(in,"%lf",&tabled4[i][j]);
	for (j=0;j<ny5;j++) for (i=0;i<nx5;i++) fscanf(in,"%lf",&tablef5[i][j]);
	for (j=0;j<ny5;j++) for (i=0;i<nx5;i++) fscanf(in,"%lf",&tabled5[i][j]);
	for (j=0;j<ny6;j++) for (i=0;i<nx6;i++) fscanf(in,"%lf",&tablef6[i][j]);
	for (j=0;j<ny6;j++) for (i=0;i<nx6;i++) fscanf(in,"%lf",&tabled6[i][j]);
	for (j=0;j<ny7;j++) for (i=0;i<nx7;i++) fscanf(in,"%lf",&tablef7[i][j]);
	for (j=0;j<ny7;j++) for (i=0;i<nx7;i++) fscanf(in,"%lf",&tabled7[i][j]);

printf("Reading of data from file interpolationdata.keep complete\n");
	fclose(in);

/* Compute reciprocals of denominators for later interpolation */
	calc_recip_denom(nx1,Vx1,xdenom1);
/* printf("corners of xdenom1 %g %g %g %g\n",xdenom1[0][0],xdenom1[0][15],xdenom1[nx1-16][0],xdenom1[nx1-16][15]); */
	calc_recip_denom(ny1,Vy1,ydenom1);
/* printf("corners of ydenom1 %g %g %g %g\n",ydenom1[0][0],ydenom1[0][15],ydenom1[ny1-16][0],ydenom1[ny1-16][15]); */
	calc_recip_denom(nx2,Vx2,xdenom2);
	calc_recip_denom(ny2,Vy2,ydenom2);
	calc_recip_denom(nx3,Vx3,xdenom3);
	calc_recip_denom(ny3,Vy3,ydenom3);
	calc_recip_denom(nx4,Vx4,xdenom4);
	calc_recip_denom(ny4,Vy4,ydenom4);
	calc_recip_denom(nx5,Vx5,xdenom5);
	calc_recip_denom(ny5,Vy5,ydenom5);
	calc_recip_denom(nx6,Vx6,xdenom6);
	calc_recip_denom(ny6,Vy6,ydenom6);
	calc_recip_denom(nx7,Vx7,xdenom7);
	calc_recip_denom(ny7,Vy7,ydenom7);

/* Also calculate Gaussian distribution for tabulated x's for use with table4 */
	for (i=0; i<nx4; i++){
		f4_alpha2[i]=normaltail(Vx4[i]/sqrt(2.));
		d4_alpha2[i]=1/(sqrt(4*pi))*exp(-Vx4[i]*Vx4[i]*.25);
	}
/*printf("x4,f4_alpha2,d4_alpha2 %g\n",x4,f4_alpha2,d4_alpha2) */

}			/* end of initialization */
/* ===================================================== */
/* Store standard numbers which vary with alpha */
previous_alpha=alpha;
distributiontabulated=FALSE;
/* printf("SETALPHA: alpha= %f\n",alpha) */

alphasmall= alpha < .5;
alphalarge= alpha > 1.7;

/* Interpolate appropriate table(s) and store the result(s) */
/* Case when alpha < .5 */
if(alphasmall){
	calpha=1-alpha;
	absam1=fabs(calpha);
	a1=alpha/absam1;
	sa2=(2-alpha)/(2*alpha);
	nu=pow(absam1,(-1/alpha));
	eta=tan(hpi*alpha);
	Calpha_C=exp(LogGamma(alpha))*sin(pi*alpha)/pi;
	ximid=.2/alpha;
/*	midpoint=pow((absam1/ximid),(1/a1))*alpha */
	logmidpoint=log(absam1/ximid)/a1+log(alpha);
/*	xi=absam1*(z/alpha)**a1) */
/*	printf("logmidpoint %g\n",logmidpoint); */
	Clogd=log(nu/sqrt(2*pi*alpha));
/*printf("eta,Calpha_C %g %g %g\n",eta,Calpha_C); */

	interpolate_over_alpha(nx1,ny1,Vy1,alpha,tablef1,tabled1,f1,d1,ydenom1);
	interpolate_over_alpha(nx2,ny2,Vy2,alpha,tablef2,tabled2,f2,d2,ydenom2);
	interpolate_over_alpha(nx3,ny3,Vy3,alpha,tablef3,tabled3,f3,d3,ydenom3);
	}
else{
/* Case when alpha > .5 */
	if(alpha>1) alphastar=1/alpha;
	else alphastar=alpha;
	ximid=.4;
	if(alphastar == 1.){
		midpoint=(-log(hpi*ximid)-1)/hpi;
		nu=1;
		eta=0;
		logscalef=log(hpi);
/* Lower limit where xi=10**30; take density to be zero below here */
		xlowlimit=-(1+log(hpi*1.E30))/hpi;
/*		printf("xlowlimit (alpha=1) %g\n",xlowlimit); */
		}
	else{
		calpha=1-alpha;
		ratio=calpha/alpha;
		angle=hpi*calpha;
		sinangle=sin(angle);
		C1=alpha/sinangle;
		C2=calpha/alpha;
		C3=calpha/sinangle;
		C4=(2*pow(sin(.25*pi*calpha),2)-calpha)/sinangle;
		midpoint=C1*expoffset(C2*log(C3/ximid))+C4;
		nu=pow(fabs(1-alpha),(-1/alpha));
		eta=1/tan(angle);
/* Lower limit where xi=10**30; take density to be zero below here */
		xlowlimit=C1*expoffset(C2*log(C3/1.E30))+C4;
/*		printf("xlowlimit %g\n",xlowlimit); */

/* Note: when nu is near 1, precision would be improved by avoiding the */
/*	subtraction implicit in use of logscalef and Clogd */

		logscalef=log(fabs(sinangle))/alpha;
/*printf("logscalef= %g\n",logscalef); */
/*printf("midpoint= %g\n",midpoint); */
/*printf("nu= %g\n",nu); */
/*printf("eta= %g\n",eta); */
/*printf("C1,C2,C3,C4= %g %g %g %g\n",C1,C2,C3,C4); */
	}
	sa2=(2-alpha)/(2*alpha);
	Clogd=log(nu/sqrt(2*pi*alpha));
/*printf("Clogd %g\n",Clogd); */
	Calpha_M=2*exp(LogGamma(alpha))*sin(hpi*alpha)/pi;
/*printf("Calpha_M %g\n",Calpha_M); */

	interpolate_over_alpha(nx1,ny1,Vy1,alphastar,
		tablef1,tabled1,f1,d1,ydenom1);
	if(alphalarge){
		interpolate_over_alpha(nx4,ny4,Vy4,alpha,
			tablef4,tabled4,f4,d4,ydenom4);
		interpolate_over_alpha(nx5,ny5,Vy5,alpha,
			tablef5,tabled5,f5,d5,ydenom5);
		}
	else{
		interpolate_over_alpha(nx6,ny6,Vy6,alpha,
			tablef6,tabled6,f6,d6,ydenom6);
		interpolate_over_alpha(nx7,ny7,Vy7,alpha,
			tablef7,tabled7,f7,d7,ydenom7);
	}
}
}
/*========================================================================= */
void tailsMSS(int n,double x[],double d[],double logd[],double F[],
	double logF[],double cF[],double logcF[],
	double alpha,double location,double logscale)

// Only need to return logd,F and cF.
// For left-skewed, need to swap F and cF.
{
/* Computes density, distribution function and complement for a maximally skew
	stable distribution skewed to the right */
/* When alpha < 0.5:
	MSS variable is exp(logscale)*(parametrization C standard)-location
	log MSS variable is exp{location-exp(logscale)*(parametrization C standard)} 
   When alpha >=0.5:
	MSS variable is exp(logscale)*(parametrization M=S0 standard)-location
	log MSS variable is exp{location-exp(logscale)*(parametrization M=S0 standard)}
   In both cases the log MSS variable = exp( - MSS variable) and
		MSS variable = -log( log MSS variable).		*/

static const double roothalf=.7071067811865475244;
static const double log_density_mode2=-1.2655121234846454;
double z,y,dy,difference,logz,t,temp,temp2,approx,scale;
int i;
/* printf("Checks (1/2) %f\n",roothalf*roothalf); */
/* printf("log_density_mode2 %.16f\n",log(.5/sqrt(pi))); */

/*printf("(tailsMSS) alpha %f\n",alpha); */
/* If appropriate, set up for new alpha */
setalpha(alpha);

/* Case when alpha is precisely 2 */
scale=exp(logscale);
if (alpha == 2){
	for(i=0; i<n; i++){
		z=(x[i]-location)/scale;
		logd[i]=log_density_mode2-.25*z*z-logscale;
		d[i]=exp(logd[i]);
		if(z<0){
			F[i]=normaltail(-z*roothalf);
			cF[i]=1-F[i];
		}
		else{
			cF[i]=normaltail(z*roothalf);
		F[i]=1-cF[i];
		}
	}
}
/* Case when alpha is less than 0.5 */
else if(alphasmall){
    for(i=0; i<n; i++){
	difference=x[i]-location;
	if(difference <= 0){
		F[i]=0.;
		logF[i]=neglarge;
		cF[i]=1.;
		logcF[i]=0.;
		d[i]=0.;
		logd[i]=neglarge;
	}
	else{
		logz=log(difference)-logscale;
/*printf("logz %g\n",logz) */

	 /*Case covered by table 1: small x */
		if(logz<logmidpoint){
		/*xi=absam1*(z/alpha)**(-a1) */
			xi=absam1*exp(-a1*(logz-log(alpha)));
			t=.2/(alpha*xi);
/*printf("xi,t %g %g\n",xi,t); */
			interpolate(t,&ffound,&dfound,nx1,Vx1,f1,d1,xdenom1);
/*printf("ffound,dfound %g %g\n",ffound,dfound); */
			logd[i]=Clogd+sa2*log(xi)-xi+log(dfound)-logscale;
			if(d[i]<logposlarge) d[i]=exp(logd[i]);
			else d[i]=poslarge;
			logF[i]=-.5*log(2*pi*alpha*xi)-xi+log(ffound);
			F[i]=exp(logF[i]);
			logcF[i]=logoffset(-F[i]);
			cF[i]=1.-F[i];
		}

	 /*Case covered by table 2: small alpha, middle range for x */
		else if(logz<0.){
			temp2=exp(-alpha*logmidpoint)-1;
/*independent of z */
			temp=exp(-alpha*logz); 	/*= z**(-alpha) */
			t=(temp-1)/temp2;
/*printf("t %g\n",t); */
			interpolate(t,&ffound,&dfound,nx2,Vx2,f2,d2,xdenom2);
/*printf("ffound,dfound %g %g\n",ffound,dfound); */
			logF[i]=log(ffound*(1+EulersC*alpha*temp))-temp;
/*printf("F[i] %g\n",F[i]); */
			F[i]=exp(logF[i]);
			logcF[i]=logoffset(-F[i]);
			cF[i]=1-F[i];
			logd[i]=log(dfound*(1-EulersC*alpha+EulersC*alpha*temp)*alpha*temp)-temp-logz-logscale;
			d[i]=exp(logd[i]);
		}
	 /*Case covered by table 3: small alpha, upper range for x */
		else{
			t=exp(-alpha*logz);
/*= z**(-alpha) */
/*printf("t %g\n",t); */
			interpolate(t,&ffound,&dfound,nx3,Vx3,f3,d3,xdenom3);
/*printf("ffound,dfound %g %g\n",ffound,dfound); */
			approx=Calpha_C*t;
			cF[i]=ffound*approx;
			logcF[i]=log(cF[i]);
			F[i]=1-cF[i];
			logF[i]=logoffset(-cF[i]);
			logd[i]=log(dfound*approx*alpha)-logz-logscale;
			d[i]=exp(logd[i]);
		}
	}
    }
}

/* Case when alpha is between 1.7 and 2 */
else if(alphalarge){
    scale=exp(logscale);
    for(i=0; i<n; i++){
	z=(x[i]-location)/scale;
/*printf("z= %g\n",z); */

/* Case when z is below limit where xi can be calculated */
	if(z<xlowlimit){
		F[i]=0.;
		logF[i]=neglarge;
		cF[i]=1.;
		logcF[i]=0.;
		d[i]=0.;
		logd[i]=neglarge;
	}
/* Case covered by table 1: low range for x */
	else if(z<midpoint){
		xi=C3/exp(logoffset((z-C4)/C1)/C2);
		t=.2/(alphastar*xi);
/*printf("xi,t %g %g\n",xi,t); */
		interpolate(t,&ffound,&dfound,nx1,Vx1,f1,d1,xdenom1);
/*printf("ffound,dfound %g %g\n",ffound,dfound); */
		logd[i]=Clogd+sa2*log(xi)-xi+log(dfound)-logscale+logscalef;
/*printf("logd[i]= %g\n",logd[i]); */
		d[i]=exp(logd[i]);
/*printf("d[i]= %g\n",d[i]); */
		logF[i]=-.5*log(2*pi*alpha*xi)-xi+log(ffound);
/*printf("logF[i]= %g\n",logF[i]); */
		F[i]=exp(logF[i]);
/*printf("F[i]= %g\n",F[i]); */
		logcF[i]=logoffset(-F[i]);
/*printf("logcF[i]= %g\n",logcF[i]); */
		cF[i]=1.-F[i];
/*printf("cF[i]= %g\n",cF[i]); */
	}

/* Case covered by table 4: large alpha, middle range for x */
	else if(z<20.){
/*printf("z %g\n",z); */
/*printf("z %g\n",z); */
		interpolate(z,&ffound,&dfound,nx4,Vx4,f4,d4,xdenom4);
/*printf("ffound,dfound %g %g\n",ffound,dfound); */
/*printf("ffound,dfound %g %g\n",ffound,dfound); */
		cF[i]=normaltail(z*sqrt(.5))+ffound*(2-alpha);
		logcF[i]=log(cF[i]);
		F[i]=1.-cF[i];
		logF[i]=logoffset(-cF[i]);
//		d[i]=exp(log_density_mode2-.25*z*z-logscale)+dfound*(2-alpha);
		d[i]=(exp(log_density_mode2-.25*z*z)+dfound*(2-alpha))/scale;
		logd[i]=log(d[i]);
	}

/* Case covered by table 5: large alpha, upper range for x */
	else{
/* Find y such that x=y+eta*(y**(1-alpha)-1) */
		y=z;
/* A good first approximation */
		do{
			dy=(z-y-expoffset(calpha*log(y))*eta)/(1+eta*calpha*pow(y,(-alpha)));
			y=y+dy;
/*printf("y,dy= %g %g\n",y,dy); */
		} while(fabs(dy)>1.e-10*y);
		t=pow((y/19.5),(-alpha));
/*printf("t %g\n",t); */
		interpolate(t,&ffound,&dfound,nx5,Vx5,f5,d5,xdenom5);
/*printf("ffound,dfound %g %g\n",ffound,dfound); */
		logapprox=log(Calpha_M)-alpha*log(y);
		logcF[i]=logapprox+log(ffound);
		cF[i]=exp(logcF[i]);
		F[i]=1.-cF[i];
		logF[i]=logoffset(-cF[i]);
		logd[i]=logapprox-logscale+log(alpha*dfound/y);
		d[i]=exp(logd[i]);
	}
    }
}
/*======================================== */
/*Case when alpha is between 0.5 and 1.7 */
else{
    scale=exp(logscale);
    for(i=0; i<n; i++){
	z=(x[i]-location)/scale;
/*printf("alpha mid range: z= %g\n",z); */

/* Case when z is below limit where xi can be calculated */
	if(z<xlowlimit){
		F[i]=0.;
		logF[i]=neglarge;
		cF[i]=1.;
		logcF[i]=0.;
		d[i]=0.;
		logd[i]=neglarge;
	}
/* Case covered by table 1: low range for x */
	else if(z<midpoint){
		if(alpha == 1) xi=exp(-1-hpi*z)/hpi;
/*	printf("Check z,xi %g %g\n",z,xi); */
		else xi=C3/exp(logoffset((z-C4)/C1)/C2);
		t=.2/(alphastar*xi);
/*printf("xi,t %g %g\n",xi,t); */
		interpolate(t,&ffound,&dfound,nx1,Vx1,f1,d1,xdenom1);
/*printf("ffound,dfound %g %g\n",ffound,dfound); */
		logd[i]=Clogd+sa2*log(xi)-xi+log(dfound)-logscale+logscalef;
		d[i]=exp(logd[i]);
		logF[i]=-.5*log(2*pi*alpha*xi)-xi+log(ffound);
		F[i]=exp(logF[i]);
		logcF[i]=logoffset(-F[i]);
		cF[i]=1.-F[i];
	}

/* Case covered by table 7: middle range for alpha, middle range for x */
	else if(z<7.3){
		t=(z-midpoint)/(7.3-midpoint);
/*printf("z %g\n",z); */
		interpolate(t,&ffound,&dfound,nx7,Vx7,f7,d7,xdenom7);
/*printf("ffound,dfound %g %g\n",ffound,dfound); */
		logcF[i]=ffound;
		cF[i]=exp(ffound);
		F[i]=1.-cF[i];
		logF[i]=logoffset(-cF[i]);
		logd[i]=dfound-logscale;
		d[i]=exp(logd[i]);
	}

/* Case covered by table 6: middle range for alpha, upper range for x */
	else{
		if(alpha == 1){
			y=z;
			do{
				dy=(z-y-log(y)/hpi)/(1+1/(y*hpi));
				y=y+dy;
/*printf("y,dy= %g %g\n",y,dy) */
			}
			while(fabs(dy)>1.e-10*y);
		}
		else{
		/* Find y such that x=y+eta*(y**calpha-1) */
		/*=y+expoffset(calpha*log(y))/tan(hpi*calpha) */
			temp=1/tan(hpi*calpha);
			y=z;
/* A good first approximation */
/*printf("y= %g\n",y); */
			do{
				dy=(z-y-expoffset(calpha*log(y))*eta)/(1+eta*calpha*pow(y,(-alpha)));
				y=y+dy;
			}
/*printf("y,dy= %g %g\n",y,dy); */
			while(fabs(dy)>1.e-10*y);
		}
		t=pow((y/5.),(-alpha));
/*printf("t %g\n",t); */
		interpolate(t,&ffound,&dfound,nx6,Vx6,f6,d6,xdenom6);
/*printf("ffound,dfound %g %g\n",ffound,dfound); */
		logapprox=log(Calpha_M)-alpha*log(y);
		logcF[i]=logapprox+log(ffound);
		cF[i]=exp(logcF[i]);
		F[i]=1.-cF[i];
		logF[i]=logoffset(-cF[i]);
		logd[i]=logapprox-logscale+log(alpha*dfound/y);
		d[i]=exp(logd[i]);
	}
    }
}
}
/*====================================================================== */
void tailslogMSS(int n,double x[],double d[],double F[],double cF[],
	double alpha,double location,double logscale)
/* Computes density, distribution function and complement for a
	log maximally skew stable distribution */
{
double nlogx[n],logd[n],logF[n],logcF[n]; 
int i;

if(alpha < 0.5){
	for (i=0; i<n; i++){nlogx[i]=-log(x[i]);}
	tailsMSS(n,nlogx,d,logd,cF,logcF,F,logF,alpha,-location,logscale);
	for (i=0; i<n; i++){d[i]=exp(logd[i]+nlogx[i]);}
}

/* Case when alpha > 0.5 */
else{
	for (i=0; i<n; i++){nlogx[i]=-log(x[i]);}
	tailsMSS(n,nlogx,d,logd,cF,logcF,F,logF,alpha,-location,logscale);
	for (i=0; i<n; i++){d[i]=exp(logd[i]+nlogx[i]);}
//Perhaps replace with?:	for (i=0; i<n; i++){d[i]=d[i]/x[i];}
}
}
/*=========================================================================*/
void tabtails(double alpha)
{
/* Computes density and tail areas for tabulated points for standard
stable distributions.  Does nothing if alpha=2. */

/* When alpha < 0.5: parametrization  C   standard */
/* When alpha >=0.5: parametrization M=S0 standard */

double x,y,logx,t,temp,temp2,approx;
int i,j;

/* If appropriate, set up for new alpha */
setalpha(alpha);
if(distributiontabulated)return;

/* Store tabulated points using running index j, starting at j=0. */
j=-1;

/* When alpha is precisely 2, do nothing. */
if (alpha == 2){
}
/* Case when alpha is less than 0.5 */
else if(alphasmall){
//printf("\nTable 1\n");
	 /*Case covered by table 1: small x */
	for (i=1; i<nx1; i++){
		t=Vx1[i];
		xi=.2/(alpha*t);
//printf("i=%d, t=%g, xi=%g\n",i,t,xi);
		j++;
		tab_lx[j]=-log(xi/absam1)/a1+log(alpha);
		tab_ld[j]=Clogd+sa2*log(xi)-xi+log(d1[i]);
		tab_lF[j]=-.5*log(2*pi*alpha*xi)-xi+log(f1[i]);
		tab_F[j]=exp(tab_lF[j]);
		tab_lcF[j]=logoffset(-tab_F[j]);
		tab_cF[j]=1.-tab_F[j];
//printf("t=%g log x=%g log den=%g F=%g, cF=%g\n",t,tab_lx[j],tab_ld[j],tab_F[j],tab_cF[j]);
	}

	/*Case covered by table 2: small alpha, middle range for x */
	temp2=exp(-alpha*logmidpoint)-1; /*independent of z */
//printf("\nTable 2\n");
	for (i=nx2-1; i>=0; i=i-1){
		t=Vx2[i];
		temp=t*temp2;
		logx=-logoffset(temp)/alpha;
		if(logx > tab_lx[j]){
#ifdef CHECK
if(i == nx2-1)printf("No overlap at start of table 2; tab_lx[j]=%g, logx=%g\n",
	tab_lx[j],logx);
#endif
			j++;
			tab_lx[j]=logx;
			temp=temp+1;
			tab_lF[j]=log(f2[i]*(1+EulersC*alpha*temp))-temp;
			tab_F[j]=exp(tab_lF[j]);
			tab_lcF[j]=logoffset(-tab_F[j]);
			tab_cF[j]=1-tab_F[j];
			tab_ld[j]=log(d2[i]*(1-EulersC*alpha+EulersC*alpha*temp)*alpha*temp)-temp-tab_lx[j];
//printf("t=%g log x=%g log den=%g F=%g, cF=%g\n",t,tab_lx[j],tab_ld[j],tab_F[j],tab_cF[j]);
		}
	}
	/*Case covered by table 3: small alpha, upper range for x */
//printf("\nTable 3\n");
	for (i=nx3-1; i>0; i=i-1){
		t=Vx2[i];
		logx=-log(t)/alpha;
		if(logx > tab_lx[j]){
#ifdef CHECK
if(i == nx3-1)printf("No overlap at start of table 3; tab_lx[j]=%g, logx=%g\n",
	tab_lx[j],logx);
#endif
			j++;
			tab_lx[j]=logx;
			approx=Calpha_C*t;
			tab_cF[j]=f3[i]*approx;
			tab_lcF[j]=log(tab_cF[j]);
			tab_F[j]=1-tab_cF[j];
			tab_lF[j]=logoffset(-tab_cF[j]);
			tab_ld[j]=log(d3[i]*approx*alpha)-tab_lx[j];
//printf("t=%g log x=%g log den=%g F=%g, cF=%g\n",t,tab_lx[j],tab_ld[j],tab_F[j],tab_cF[j]);
		}
	}
}

/* Case when alpha is between 1.7 and 2 */
else if(alphalarge){
//printf("\nTable 1(large alpha)\n");
	/* Case covered by table 1: low range for x */
	for (i=1; i<nx1; i++){
		t=Vx1[i];
		xi=.2/(alphastar*t);
		j++;
		tab_x[j]=C1*expoffset(C2*log(C3/xi))+C4;
		tab_ld[j]=Clogd+sa2*log(xi)-xi+log(d1[i])+logscalef;
		tab_density[j]=exp(tab_ld[j]);
		tab_lF[j]=-.5*log(2*pi*alpha*xi)-xi+log(f1[i]);
		tab_F[j]=exp(tab_lF[j]);
		tab_lcF[j]=logoffset(-tab_F[j]);
		tab_cF[j]=1.-tab_F[j];
//printf("t=%g x=%g log den=%g F=%g, cF=%g\n",t,tab_x[j],tab_ld[j],tab_F[j],tab_cF[j]);
	}

	/* Case covered by table 4: large alpha, middle range for x */
//printf("\nTable 4\n");
	for (i=0; i<nx4; i++){
		t=Vx4[i];
		if(t > tab_x[j]){
#ifdef CHECK
if(i == 0)printf("No overlap at start of table 4; tab_x[j]=%g, t=%g\n",
	tab_x[j],t);
#endif
			j++;
			tab_x[j]=t;
			tab_cF[j]=f4_alpha2[i]+f4[i]*(2-alpha);
			tab_lcF[j]=log(tab_cF[j]);
			tab_F[j]=1.-tab_cF[j];
			tab_lF[j]=logoffset(-tab_cF[j]);
			tab_density[j]=d4_alpha2[i]+d4[i]*(2-alpha);
//printf("t=%g x=%g log den=%g F=%g, cF=%g\n",t,tab_x[j],tab_ld[j],tab_F[j],tab_cF[j]);
		}
	}

/* Case covered by table 5: large alpha, upper range for x */
//printf("\nTable 5\n");
	for (i=nx5-1; i>0; i=i-1){
		t=Vx5[i];
		y=pow(t,-1/alpha)*19.5;
		x=y+eta*(pow(y,(1-alpha))-1);
		if(x > tab_x[j]){
#ifdef CHECK
if(i == nx5-1)printf("No overlap at start of table 5; tab_x[j]=%g, x=%g\n",
	tab_x[j],x);
#endif
			j++;
			tab_x[j]=x;
			logapprox=log(Calpha_M)-alpha*log(y);
			tab_lcF[j]=logapprox+log(f5[i]);
			tab_cF[j]=exp(tab_lcF[j]);
			tab_F[j]=1.-tab_cF[j];
			tab_lF[j]=logoffset(-tab_cF[j]);
			tab_ld[j]=logapprox+log(alpha*d5[i]/y);
			tab_density[j]=exp(tab_ld[j]);
//printf("t=%g x=%g log den=%g F=%g, cF=%g\n",t,tab_x[j],tab_ld[j],tab_F[j],tab_cF[j]);
		}
	}
}
/*======================================== */
/*Case when alpha is between 0.5 and 1.7 */
else{
//printf("\nTable 1(mid alpha)\n");
	/* Case covered by table 1: low range for x */
	for (i=1; i<nx1; i++){
		t=Vx1[i];
		xi=.2/(alphastar*t);
		j++;
		tab_x[j]=C1*expoffset(C2*log(C3/xi))+C4;
		if(alpha == 1) tab_x[j]=(-log(hpi*xi)-1)/hpi;
			else tab_x[j]=C1*expoffset(C2*log(C3/xi))+C4;
		tab_ld[j]=Clogd+sa2*log(xi)-xi+log(d1[i])+logscalef;
		tab_density[j]=exp(tab_ld[j]);
		tab_lF[j]=-.5*log(2*pi*alpha*xi)-xi+log(f1[i]);
		tab_F[j]=exp(tab_lF[j]);
		tab_lcF[j]=logoffset(-tab_F[j]);
		tab_cF[j]=1.-tab_F[j];
//printf("t=%g x=%g log den=%g F=%g, cF=%g\n",t,tab_x[j],tab_ld[j],tab_F[j],tab_cF[j]);
	}

/* Case covered by table 7: middle range for alpha, middle range for x */
//printf("\nTable 7\n");
	for (i=0; i<nx7; i++){
		t=Vx7[i];
		x=t*(7.3-midpoint)+midpoint;
		if(x > tab_x[j]){
#ifdef CHECK
if(i == 0)printf("No overlap at start of table 7; tab_x[j]=%g, x=%g\n",
	tab_x[j],x);
#endif
			j++;
			tab_x[j]=x;
			tab_lcF[j]=f7[i];
			tab_cF[j]=exp(f7[i]);
			tab_F[j]=1.-tab_cF[j];
			tab_lF[j]=logoffset(-tab_cF[j]);
			tab_ld[j]=d7[i];
			tab_density[j]=exp(tab_ld[j]);
//printf("t=%g x=%g log den=%g F=%g, cF=%g\n",t,tab_x[j],tab_ld[j],tab_F[j],tab_cF[j]);
		}
	}

/* Case covered by table 6: middle range for alpha, upper range for x */
//printf("\nTable 6\n");
	for (i=nx6-1; i>0; i=i-1){
		t=Vx6[i];
		y=pow(t,-1/alpha)*5.;
		if(alpha==1)x=y+log(y)/hpi;
		else x=y+eta*(pow(y,(1-alpha))-1);
		if(x > tab_x[j]){
#ifdef CHECK
if(i == nx6-1)printf("No overlap at start of table 6; tab_x[j]=%g, x=%g\n",
	tab_x[j],x);
#endif
			j++;
			tab_x[j]=x;
			logapprox=log(Calpha_M)-alpha*log(y);
			tab_lcF[j]=logapprox+log(f6[i]);
			tab_cF[j]=exp(tab_lcF[j]);
			tab_F[j]=1.-tab_cF[j];
			tab_lF[j]=logoffset(-tab_cF[j]);
			tab_ld[j]=logapprox+log(alpha*d6[i]/y);
			tab_density[j]=exp(tab_ld[j]);
//printf("t=%g x=%g log den=%g F=%g, cF=%g\n",t,tab_x[j],tab_ld[j],tab_F[j],tab_cF[j]);
		}
	}
}
ntabs=j+1;
distributiontabulated=TRUE;
}
/*====================================================================== */
void setuplogMSS(double alpha,double location,double logscale,
	double relative_precision)
/* Tabulates option prices based on a log maximally skew stable distribution. */
/* Does nothing if alpha=2. */
/* This procedure may fail for very small alpha when the probability
distribution is near to being a two-point discrete distribution with
nearly all of the probability at one points. */
//Remember to modify optionslogmss since list is now in ascending order. 
//Also check handling of tails of distributions.

/* When alpha is near to 2 (say, alpha > 1.7) then it would probably be
adequate to used the tabulated points without even considering what
precision would be achieved. */

{
int i,j,itab,ilist,istack,OK;
double abs_precision;
double h,yh,error,coeffx4_last,coeffx4_here,cutoff,lastcF;
double temp_density,temp_logdensity,temp_logcF,temp_F,temp_logF;
double approx,nlogstrike,hsq,sum,temp;
int ltab,rtab,offset;
double maxcF,leftstrike,rightstrike,leftdensity,rightdensity,leftcF,rightcF;
double maxx4,laststrike,lastdensity,minhy,maxhy,nlogleftstrike,nlogrightstrike;
double keepstrike,keepcF,keepdensity,tabstrike;
double factor,c,minh,maxh,maxstrike,est_right_tail,term;

/* ----------- Meanings of some variables ----------------
abs_precision is the absolute precision required for each subregion
coeffx4_here is a local estimate of the coefficient of x^4
coeffx4_right is an estimate of the coefficient of x^4 from the right
Note that densities are of the log MSS variables.
*/

if(alpha == 2.)return;
if(alpha == previous_alpha & location == previous_location &
	 logscale == previous_logscale & relative_precision == previous_precision)return;
tabtails(alpha);
meanprice=expectationlogMSS(alpha,location,logscale,1.);
//printf("meanprice=%g\n",meanprice);
abs_precision=meanprice*relative_precision/MAX_NSTRIKES;
//printf("abs_precision=%g\n",abs_precision);
//======================================= ==============================
/* Note that tabulated values are in order of increasing x which is decreasing
order of strike price, since strike=exp{location-exp(logscale)*x}. */

/* Use ilist as index of values being listed.  Start at end. */
ilist=MAX_NSTRIKES-1;

/* ========== Find 3 points near mean price for starting loops =========== */
/* Find itab for first tabulated strike below meanprice */
if(alphasmall){
	/* Condition is location-exp(logscale+tab_lx[itab]) <log(meanprice)
	i.e. logscale+tab_lx[itab] > log(location-log(meanprice))
	i.e. tab_lx[itab] > log(location-log(meanprice)) - logscale */
	cutoff=log(location-log(meanprice)) - logscale;
//	printf("cutoff: %g\n",cutoff);
	for(itab=0; itab<ntabs; itab++)	if(tab_lx[itab] > cutoff)break;
	if(itab==ntabs-1)printf("SHOULDN'T HAVE HAPPENED!\n");
	nlogstrike=location-exp(logscale+tab_lx[itab]);
	l_strike[ilist]=exp(nlogstrike);
	l_density[ilist] = exp(tab_ld[itab]-logscale-nlogstrike);
}
else{
	explogscale=exp(logscale);
	/* Condition is exp(location-explogscale*tab_x[itab]) < meanprice
	i.e. explogscale*tab_x[itab] > location-log(meanprice) */
	cutoff=(location-log(meanprice))/explogscale;
//	printf("cutoff: %g\n",cutoff);
	for(itab=0; itab<ntabs; itab++)	if(tab_x[itab] > cutoff)break;
	l_strike[ilist]=exp(location-explogscale*tab_x[itab]);
	l_density[ilist] = tab_density[itab]/(explogscale*l_strike[ilist]);

}
l_cF[ilist]=tab_F[itab];
//printf("itab: %d\n",itab);
//printf("strike: %g\n",l_strike[ilist]);
//printf("cF: %g\n",l_cF[ilist]);
//printf("density: %g\n",l_density[ilist]);

/* If this strike is smaller than .6*meanprice, use non-tabulated value. */
if(l_strike[ilist] < .6*meanprice){
	l_strike[ilist]=meanprice;
	nlogstrike=-log(l_strike[ilist]);
        tailsMSS(1,&nlogstrike,&temp_density,&temp_logdensity,&l_cF[ilist],&temp_logcF,&temp_F,&temp_logF,alpha,-location,logscale);
	l_density[ilist]=exp(temp_logdensity+nlogstrike);
//	printf("strike: %g\n",l_strike[ilist]);
//	printf("cF: %g\n",l_cF[ilist]);
//	printf("density: %g\n",l_density[ilist]);
	ltab=itab;
}
else{
	ltab=itab+1;
}
rtab=itab-1;
/*	ltab is nearest tabulated point to left and
	rtab is nearest tabulated point to right	*/

/* Find tabulated point to right with cF < .99*l_cF[ilist] */
cutoff=.99*l_cF[ilist];
//printf("cutoff: %g\n",cutoff);
for(itab=rtab; itab>=0; itab=itab-1)	if(tab_F[itab] < cutoff)break;
if(itab==0)printf("THIS SHOULDN'T HAVE HAPPENED!\n");
if(alphasmall){
	nlogstrike=exp(logscale+tab_lx[itab])-location;
	rightstrike=exp(-nlogstrike);
	rightdensity = exp(tab_ld[itab]-logscale-nlogstrike);
}
else{
	rightstrike=exp(location-explogscale*tab_x[itab]);
	rightdensity= tab_density[itab]/(explogscale*rightstrike);
}
rightcF=tab_F[itab];
//printf("rightstrike: %g\n",rightstrike);
//printf("rightcF: %g\n",rightcF);
//printf("rightdensity: %g\n",rightdensity);

/* If this strike is greater than 2*meanprice, use non-tabulated value. */
if(rightstrike > 2.*meanprice){
	rightstrike=2.*meanprice;
	nlogstrike=-log(rightstrike);
        tailsMSS(1,&nlogstrike,&temp_density,&temp_logdensity,&rightcF,&temp_logcF,&temp_F,&temp_logF,alpha,-location,logscale);
	rightdensity=exp(temp_logdensity+nlogstrike);
//	printf("rightstrike: %g\n",rightstrike);
//	printf("rightcF: %g\n",rightcF);
//	printf("rightdensity: %g\n",rightdensity);
}
/* Keep values needed for loop to right about 150 lines below */
keepstrike=rightstrike;
keepcF=rightcF;
keepdensity=rightdensity;

/* Find tabulated point to left with cF at least .01+.99*l_cF[ilist] */
cutoff=.01+cutoff;
//printf("cutoff: %g\n",cutoff);
for(itab=ltab; itab<ntabs; itab++)	if(tab_F[itab] > cutoff)break;
if(alphasmall){
	nlogstrike=exp(logscale+tab_lx[ltab])-location;
	leftstrike=exp(-nlogstrike);
	leftdensity = exp(tab_ld[ltab]-logscale-nlogstrike);
}
else{
	leftstrike=exp(location-explogscale*tab_x[ltab]);
	leftdensity= tab_density[ltab]/(explogscale*leftstrike);
}
leftcF=tab_F[ltab];
//printf("leftstrike: %g\n",leftstrike);
//printf("leftcF: %g\n",leftcF);
//printf("leftdensity: %g\n",leftdensity);

/* If this strike is less than .5*meanprice, use non-tabulated value. */
if(leftstrike < .5*meanprice){
	leftstrike=.5*meanprice;
	nlogstrike=-log(leftstrike);
        tailsMSS(1,&nlogstrike,&temp_density,&temp_logdensity,&leftcF,&temp_logcF,&temp_F,&temp_logF,alpha,-location,logscale);
	leftdensity=exp(temp_logdensity+nlogstrike);
//	printf("leftstrike: %g\n",leftstrike);
//	printf("leftcF: %g\n",leftcF);
//	printf("leftdensity: %g\n",leftdensity);
}

/* ======= Initialize loop to left =============== */
/* Integrate to left using formulae for x as function of y. */
/* One approximation to coefficient of x^4 disregards density on left */ 
coeffx4_last=fabs(x4coeff(rightstrike,-1/rightdensity,
	l_cF[ilist]-rightcF,l_strike[ilist],-1/l_density[ilist],
	leftcF-rightcF,leftstrike));
//printf("coeffx4_last=%g\n",coeffx4_last);
/* Another approximation to coefficient of x^4 disregards density on right */ 
coeffx4_here=fabs(x4coeff(leftstrike,-1/leftdensity,
	l_cF[ilist]-leftcF,l_strike[ilist],-1/l_density[ilist],
	rightcF-leftcF,rightstrike));
//printf("coeffx4_here=%g\n",coeffx4_here);
/* Find larger absolute estimate of coefficient of x^4 */
maxx4=MAX(coeffx4_here,coeffx4_last);

laststrike=rightstrike;
lastcF=rightcF;
rightstrike=l_strike[ilist];
nlogrightstrike=-log(rightstrike);
rightcF=l_cF[ilist];
rightdensity=l_density[ilist];
/* ========================= Loop to left ======================= */
/* On start of each loop, know rightstrike, nlogrightstrike, rightcF,
   rightdensity, laststrike, lastcF, coeffx4_here, maxx4 and
   ltab (index of tabulated point to left of rightstrike) */
do{
//	printf("rightstrike=%g\n",rightstrike);
//	printf("nlogrightstrike=%g\n",nlogrightstrike);
//	printf("rightcF=%g\n",rightcF);
//	printf("rightdensity=%g\n",rightdensity);
//	printf("laststrike=%g\n",laststrike);
//	printf("lastcF=%g\n",lastcF);
//	printf("coeffx4_here=%g\n",coeffx4_here);
//	printf("maxx4=%g\n",maxx4);
//	printf("ltab=%d\n",ltab);
	/* Find spacing of points so that rounding error in computing
	approximation to coefficient of x^4 is not excessive.
	If spacing in x is h then rounding error for an error of 1.e-10
	in a cF is approximately 1.e-10 * h**-4 * density**-5
	(For right, remove density factor) */
	
	minhy=pow(1.e-10/(MAX(l_density[ilist],rightdensity)*maxx4),.25);
//	printf("minhy=%g\n",minhy);

	/* Find distance such that it is likely that numerical
	integration over that distance will be sufficiently accurate. */
	maxhy=pow(30.*abs_precision/maxx4,.2);
	if(maxhy<minhy) printf("SURPRISING: maxhy=%g < minhy=%g\n",maxhy,minhy);
//	printf("maxhy=%g\n",maxhy);
	maxcF=maxhy+rightcF;
//	printf("maxcF=%g\n",maxcF);
	if(maxcF>= 1) break;	/* Integration to strike=0 should be OK */

	/* If there is a suitable tabulated point, integrate to there. */
	for(itab=ltab; itab<ntabs; itab++){if(tab_F[itab] > maxcF)break;}
	itab=itab-1;
//	printf("itab=%d\n",itab);
	if(itab>ltab & tab_F[itab] > minhy +l_cF[ilist]){
		ltab=itab+1;
		if(alphasmall){
			nlogleftstrike=exp(logscale+tab_lx[itab])-location;
			leftstrike=exp(-nlogleftstrike);
			leftdensity =exp(tab_ld[itab]-logscale+nlogleftstrike);
		}
		else{
			nlogleftstrike=explogscale*tab_x[itab]-location;
			leftstrike=exp(-nlogleftstrike);
			leftdensity=tab_density[itab]/(explogscale*leftstrike);
		}
		leftcF=tab_F[itab];
	}
	else {
	/* If no such point tabulated, find new point at distance approx maxhy*/
	/* Note that this might fail if run out of tabulated points to left */
		ltab=itab;
		itab=itab+1;	/* To interpolate using far point */
		if(alphasmall){nlogstrike=exp(logscale+tab_lx[itab])-location;}
		else{nlogstrike=explogscale*tab_x[itab]-location;}
		c=tab_F[itab]-maxcF;
//		printf("c=%g\n",c);
		nlogleftstrike=(maxhy*nlogstrike+c*nlogrightstrike)/(maxhy+c);
//		printf("nlogleftstrike=%g\n",nlogleftstrike);
		leftstrike=exp(-nlogleftstrike);
        	tailsMSS(1,&nlogleftstrike,&temp_density,&temp_logdensity,&leftcF,&temp_logcF,&temp_F,&temp_logF,alpha,-location,logscale);
		leftdensity=exp(temp_logdensity+nlogleftstrike);

	}
	/* Compute contribution to call option from this region */
	/* Attribute to right hand side end of region.	*/
	h=rightstrike - leftstrike;
	yh=leftcF-rightcF;
//	printf("h=%g\n",h);
	l_call[ilist]=h*rightcF +.5*yh*(h+(1/6.)*yh*(1/leftdensity-1/rightdensity));
//	printf("Contrib. to integral=%g\n",l_call[ilist]);

	ilist=ilist-1;
	if(ilist<0){printf("Too many strikes needed\n");ilist=1/0;}
	l_cF[ilist]=leftcF;
	l_strike[ilist]=leftstrike;
	l_density[ilist] =leftdensity;
//	printf("strike: %g\n",l_strike[ilist]);
//	printf("cF: %g\n",l_cF[ilist]);
//	printf("density: %g\n",l_density[ilist]);

	/* Compute new approximation to coefficient of x^4. */
	coeffx4_last=coeffx4_here;
	coeffx4_here=fabs(x4coeff(leftstrike,-1/leftdensity,
		rightcF-leftcF,rightstrike,-1/rightdensity,
		lastcF-leftcF,laststrike));
//	printf("coeffx4_here=%g\n",coeffx4_here);
	/* Find conservative (large) estimate of coefficient of x^4 */
	if(coeffx4_last>coeffx4_here) maxx4=coeffx4_last;
	else maxx4=(coeffx4_here/coeffx4_last)*coeffx4_here;

	/* Check this by making minor amendments to some calls to x4coeff. */
/*	temp=fabs(x4coeff(leftstrike+1.e-10,-1/leftdensity,
		rightcF-leftcF,rightstrike,-1/rightdensity,
		lastcF-leftcF,laststrike));
	printf("  coeffx4 change=%g\n",coeffx4_here-temp);
	temp=fabs(x4coeff(leftstrike,-1/leftdensity,
		1.e-10+rightcF-leftcF,rightstrike,-1/rightdensity,
		lastcF-leftcF,laststrike));
	printf("  coeffx4 change=%g\n",coeffx4_here-temp);
	temp=fabs(x4coeff(leftstrike,-1/leftdensity,
		rightcF-leftcF,rightstrike,-1/rightdensity,
		1.e-10+lastcF-leftcF,laststrike));
	printf("  coeffx4 change=%g\n",coeffx4_here-temp);
*/
	/* Move data as required at the start of the loop */
	laststrike=rightstrike;
	lastcF=rightcF;
	rightstrike=leftstrike;
	nlogrightstrike=nlogleftstrike;
	rightcF=leftcF;
	rightdensity=leftdensity;
	
	/* If contribution from tail is acceptably small, stop. */
} while ((1.-leftcF)*leftstrike>abs_precision);
/* Add zero to list of strikes */
l_call[ilist]=.5*(1.+leftcF)*leftstrike;
ilist=ilist-1;
if(ilist<0){printf("Too many strikes needed\n");ilist=1/0;}
l_strike[ilist]=0;
l_cF[ilist]=1;
l_density[ilist]=poslarge;
l_call[ilist]=0;
	
/* Now move stored values to left hand end of arrays */
offset=ilist;
//printf("offset=%d\n",offset);
for(i=0; i<MAX_NSTRIKES-ilist ; i++){
	j=i+offset;
	l_strike[i]=l_strike[j];
//	printf("strike of %g moved from location %d to %d\n",l_strike[i],j,i);
	l_cF[i]=l_cF[j];
	l_density[i]=l_density[j];
	l_call[i]=l_call[j];
}
ilist=i-1;
//printf("\nAt end of loop left: ilist=%d\n",ilist);

/* ======= Initialize loop to right =============== */
/* Integrate to right using formulae for y as function of x. */
laststrike=l_strike[ilist-1];
lastcF=l_cF[ilist-1];
leftstrike=l_strike[ilist];
leftcF=l_cF[ilist];
leftdensity=l_density[ilist];
rightstrike=keepstrike;
rightcF=keepcF;
rightdensity=keepdensity;
//printf("laststrike=%g\n",laststrike);
//printf("lastcF=%g\n",lastcF);
//printf("leftstrike=%g\n",leftstrike);
//printf("leftcF=%g\n",leftcF);
//printf("leftdensity=%g\n",leftdensity);
//printf("rightstrike=%g\n",rightstrike);
//printf("rightcF=%g\n",rightcF);
//printf("rightdensity=%g\n",rightdensity);

/* One approximation to coefficient of x^4 disregards density on right */ 
coeffx4_last=fabs(x4coeff(lastcF,-l_density[ilist-1],
	leftstrike-laststrike,leftcF,-leftdensity,
	rightstrike-laststrike,rightcF));
//printf("coeffx4_last=%g\n",coeffx4_last);
/* Another approximation to coefficient of x^4 disregards density on right */ 
coeffx4_here=fabs(x4coeff(rightcF,-rightdensity,
	leftstrike-rightstrike,leftcF,-leftdensity,
	laststrike-rightstrike,lastcF));
//printf("coeffx4_here=%g\n",coeffx4_here);
/* Find larger absolute estimate of coefficient of x^4 */
maxx4=MAX(coeffx4_here,coeffx4_last);

/* ========================= Loop to right ======================= */
/* On start of each loop, know leftstrike, leftcF,
   leftdensity, laststrike, lastcF, coeffx4_here, maxx4 and
   rtab (index of tabulated point to right of leftstrike) */
do{
//	printf("leftstrike=%g\n",leftstrike);
//	printf("leftcF=%g\n",leftcF);
//	printf("leftdensity=%g\n",leftdensity);
//	printf("laststrike=%g\n",laststrike);
//	printf("lastcF=%g\n",lastcF);
//	printf("coeffx4_here=%g\n",coeffx4_here);
//	printf("maxx4=%g\n",maxx4);
//	printf("rtab=%d\n",rtab);
	/* Find spacing of points so that rounding error in computing
	approximation to coefficient of x^4 is not excessive.
	If spacing in x is h then rounding error for an error of 1.e-10
	in a cF is approximately 1.e-10 * h**-4		*/
	
	minh=pow(1.e-10/maxx4,.25);
//	printf("minh=%g\n",minh);

	/* Find distance such that it is likely that numerical
	integration over that distance will be sufficiently accurate. */
	maxh=pow(30.*abs_precision/maxx4,.2);
	if(maxh<minh) printf("SURPRISING: maxh=%g < minh=%g\n",maxh,minh);
//	printf("maxh=%g\n",maxh);
	maxstrike=maxh+leftstrike;
//	printf("maxstrike=%g\n",maxstrike);

	/* If there is a suitable tabulated point, integrate to there. */
	for(itab=rtab; itab>0; itab=itab-1){
		if(alphasmall){tabstrike=exp(location-exp(logscale+tab_lx[itab]));}
		else{tabstrike=exp(location-explogscale*tab_x[itab]);}
//		printf("itab=%d\n",itab);
//		printf("tabstrike=%g\n",tabstrike);
		if(tabstrike>maxstrike)break;}
	itab=itab+1;
//	printf("itab=%d\n",itab);
	if(alphasmall){tabstrike=exp(location-exp(logscale+tab_lx[itab]));}
	else{tabstrike=exp(location-explogscale*tab_x[itab]);}
//	printf("tabstrike=%g\n",tabstrike);
	if(itab<rtab & tabstrike > minh+leftstrike){
		rtab=itab-1;
		if(alphasmall){
			nlogrightstrike=exp(logscale+tab_lx[itab])-location;
			rightstrike=exp(-nlogrightstrike);
			rightdensity =exp(tab_ld[itab]-logscale+nlogrightstrike);
		}
		else{
			nlogrightstrike=explogscale*tab_x[itab]-location;
			rightstrike=exp(-nlogrightstrike);
			rightdensity=tab_density[itab]/(explogscale*rightstrike);
		}
		rightcF=tab_F[itab];
	}
	else {
	/* If no such point tabulated, use maxstrike. */
	/* Note that this might fail if run out of tabulated points to right */
		rightstrike=maxstrike;
		nlogrightstrike=-log(rightstrike);
        	tailsMSS(1,&nlogrightstrike,&temp_density,&temp_logdensity,&rightcF,&temp_logcF,&temp_F,&temp_logF,alpha,-location,logscale);
		rightdensity=exp(temp_logdensity+nlogrightstrike);

	}

	ilist=ilist+1;
	if(ilist==MAX_NSTRIKES){printf("Too many strikes needed\n");ilist=1/0;}
	l_cF[ilist]=rightcF;
	l_strike[ilist]=rightstrike;
	l_density[ilist] =rightdensity;
//	printf("strike: %g\n",l_strike[ilist]);
//	printf("cF: %g\n",l_cF[ilist]);
//	printf("density: %g\n",l_density[ilist]);

	/* Compute contribution to call option from this region */
	/* Attribute to right hand side end of region.	*/
	h=rightstrike - leftstrike;
//	printf("h=%g\n",h);
	l_call[ilist]=.5*h*(leftcF+rightcF+(1/6.)*h*(rightdensity-leftdensity));
//	printf("Contrib. to integral=%g\n",l_call[ilist]);

	/* Compute new approximation to coefficient of x^4. */
	coeffx4_last=coeffx4_here;
	coeffx4_here=fabs(x4coeff(rightcF,-rightdensity,
		leftstrike-rightstrike,leftcF,-leftdensity,
		laststrike-rightstrike,lastcF));
//	printf("coeffx4_here=%g\n",coeffx4_here);
	/* Find conservative (large) estimate of coefficient of x^4 */
	if(coeffx4_last>coeffx4_here) maxx4=coeffx4_last;
	else maxx4=(coeffx4_here/coeffx4_last)*coeffx4_here;

	/* Check this by making minor amendments to some calls to x4coeff. */
/*	temp=fabs(x4coeff(rightcF+1.e-10,-rightdensity,
		leftstrike-rightstrike,leftcF,-leftdensity,
		laststrike-rightstrike,lastcF));
	printf("  coeffx4 change=%g\n",coeffx4_here-temp);
	temp=fabs(x4coeff(rightcF,-rightdensity,
		leftstrike-rightstrike,leftcF+1.e-10,-leftdensity,
		laststrike-rightstrike,lastcF));
	printf("  coeffx4 change=%g\n",coeffx4_here-temp);
	temp=fabs(x4coeff(rightcF,-rightdensity,
		leftstrike-rightstrike,leftcF,-leftdensity,
		laststrike-rightstrike,lastcF+1.e-10));
	printf("  coeffx4 change=%g\n",coeffx4_here-temp);
*/	
	/* Move data as required at the start of the loop */
	laststrike=leftstrike;
	lastcF=leftcF;
	leftstrike=rightstrike;
	leftcF=rightcF;
	leftdensity=rightdensity;
	
	/* If contribution from tail is acceptably small, stop. */
} while (leftcF*leftcF>leftdensity*abs_precision);
/* Retain estimated right tail contribution */
nstrikes=ilist+1;
if(leftdensity>0.) est_right_tail=leftcF*leftcF/leftdensity;
	else est_right_tail=0;
if(!(est_right_tail > 0.)) est_right_tail=0.;
//printf("est_right_tail=%g\n",est_right_tail);
	
//======================================= ==============================
/* Now accumulate call option contributions. */
sum=est_right_tail;
for(i=nstrikes-1; i>0; i=i-1 ){
	term=l_call[i];
	l_call[i]=sum;
	sum=sum+term;
}

/* Check that contributions add to approximately the meanprice. */
printf("sum=%g, ought to be %g, discrepancy: %g\n",sum,meanprice,sum-meanprice);
if(fabs(sum/meanprice-1) > 2*relative_precision) printf("Doesn't add up\n");

/* Scale contributions so that they add to the meanprice */
factor=meanprice/sum;
for(i=1; i<nstrikes; i++ ) l_call[i]=l_call[i]*factor;
l_call[0]=meanprice;

previous_alpha=alpha;
previous_location=location;
previous_logscale=logscale;
previous_precision=relative_precision;
}
/*====================================================================== */
void optionslogMSS(int n,double x[],double call[],
	double alpha,double location,double logscale)
/* Computes values of options for a log maximally skew stable distribution */
/* Note that	strike+call=mean+put. */
{
	//Use mid3.
/* If new alpha, location or scale then use default relative precision (1.e-5). */
int i,j,lmin,lmax,mid;
double strike,logstrike,density,logdensity,cF,logcF,F,logF;
double h,yh,leftd,rightd,scale,rscale,z,mean;

if(alpha == 2.){ /* Use Black-Scholes formula  for lognormal distribution */
	scale=exp(logscale)*sqrt(2.);
	rscale=1./scale;
	mean=exp(location+exp(2*logscale));
//	printf("optionslogMSS: mean=%g\n",mean);
//	printf("optionslogMSS: scale=%g\n",scale);
	for(i=0; i<n; i++){
		if(x[i] <= 0) call[i]=mean;
		else {
			z=(log(x[i])-location)*rscale;
//			printf("optionslogMSS: x[i]=%g\n",x[i]);
//			printf("optionslogMSS: z=%g\n",z);
			call[i]=mean*(normaltail(z-scale))-x[i]*(normaltail(z));
//			printf("optionslogMSS: call[i]=%g\n",call[i]);
		}
	}
	return;
}
if(alpha != previous_alpha | location != previous_location |
		logscale != previous_logscale)
	setuplogMSS(alpha,location,logscale,1.e-5);
for(i=0; i<n; i++){
	strike=x[i];
	if(strike <= 0) call[i]=meanprice;
	else if(strike > l_strike[nstrikes-1]) call[i]=0;
	else{
		/* Find position in list of strikes by binary chop */
		lmin=0;
		lmax=nstrikes-1;
		while(lmax-lmin>1) {
			mid=(lmin+lmax)/2;
			if(strike>l_strike[mid]) lmin=mid;
			else lmax=mid;
		}
//printf("\nlmin=%d\n",lmin);
//printf("lmax=%d\n",lmax);
		logstrike=-log(strike);
        	tailsMSS(1,&logstrike,&density,&logdensity,&cF,&logcF,&F,&logF,alpha,-location,logscale);
		density=exp(logdensity+logstrike);
//printf("density=%g\n",density);
//printf("cF=%g\n",cF);
		h=l_strike[lmax]-strike;
//printf("h=%g\n",h);
		if(strike > meanprice){
			call[i]=l_call[lmax]+.5*h*(l_cF[lmax]+cF+(1/6.)*h*(l_density[lmax]-density));
		}
		else if(lmin==0) call[i]=meanprice-.5*(1.+cF)*strike;
		else {
		/* Cubic rule for x(y) on logMSS distribution */
			yh=l_cF[lmax]-cF;
			call[i]=l_call[lmax]+h*cF+.5*yh*(h+(1/6.)*yh*(1/density-1/l_density[lmax]));
		}
	}
}
}
/*====================================================================== */
/*==================================================================== */
void Sexpoffset(double *n,double *in,double *out)
{
int i,nn;
nn = *n;
for (i=0;i<nn;i++) {out[i]=expoffset(in[i]);}
}
/*====================================================================== */
void Slogoffset(double *n,double *in,double *out)
{
int i,nn;
nn = *n;
for (i=0;i<nn;i++) {out[i]=logoffset(in[i]);}
}
/*====================================================================== */
void Siidcombine(double *Snumber,double *Salpha,double *Slocation,
	double *Slogscale, double *location_product,double *logscale_product)
/* Given parameters of a log maximally skew stable distribution,
	find the location and logscale for the product of a specified
	number of such distributions. (The number may be less than unity.) */
{
double alpha,location,logscale,number;
alpha = *Salpha;
location = *Slocation;
logscale = *Slogscale;
number = *Snumber;
iidcombine(alpha,location,logscale,number,location_product,logscale_product);
}
/*==================================================================== */
void Smatch2logMSS(double *Salpha,double *Smean,double *Sstandard_deviation,
	double *location,double *logscale)
/* Given alpha, find parameters of a log maximally skew stable distribution */
/*	to fit given mean and standard deviation */
{
double mean,standard_deviation,alpha;
mean = *Smean;
standard_deviation = *Sstandard_deviation;
alpha = *Salpha;
match2logMSS(mean,standard_deviation,alpha,location,logscale);
}
/*========================================================================= */
void StailsMSS(double *Salpha,double *Slocation,double *Slogscale,
	double *Sn,double *x,double *d,double *logd,double *F,
	double *logF,double *cF,double *logcF)
/* Computes density, distribution function and complement for a maximally skew
	stable distribution skewed to the right */
{
int n;
double alpha,location,logscale;
n = *Sn;
alpha = *Salpha;
location = *Slocation;
logscale = *Slogscale;
tailsMSS(n,x,d,logd,F,logF,cF,logcF,alpha,location,logscale);
}
/*====================================================================== */
void StailslogMSS(double *Salpha,double *Slocation,double *Slogscale,
	double *Sn,double *x,double *d,double *F,double *cF)
/* Computes density, distribution function and complement for a
	log maximally skew stable distribution */
{
int n;
double alpha,location,logscale;
n = *Sn;
alpha = *Salpha;
location = *Slocation;
logscale = *Slogscale;
tailslogMSS(n,x,d,F,cF,alpha,location,logscale);
}
/*====================================================================== */
void SoptionslogMSS(double *Salpha,double *Slocation,
	double *Slogscale,double *Sn,double *x,double *call)
/* Computes values of options for a log maximally skew stable distribution */
/* Note that	strike+call=mean+put. */
{
int n;
double alpha,location,logscale;
n = *Sn;
alpha = *Salpha;
location = *Slocation;
logscale = *Slogscale;
optionslogMSS(n,x,call,alpha,location,logscale);
}

<html>
<body>
<p align="center"><font size="-2">
© Copyright 2012, CSIRO Australia<br>
Use of this web site and information available from<br>
it is subject to our </font>
<a href="http://www.csiro.au/legalnotices/disclaimer.html">
<font size="-2">Legal Notice and Disclaimer </font></a> and <a href="http://www.csiro.au/csiro/channel/pchb1,,.html"><font size="-2">Privacy
Statement</font></a>
</p>
<!-- Note: this is sourced from 'footer2.inc' -->
<!-- Clicky code for visitation statistices follows -->
<a title="Clicky Web Analytics" href="http://getclicky.com/16636"><!-- <img alt="Clicky Web Analytics" src="http://static.getclicky.com/media/links/badge.gif" border="0" /> --></a>
<script src="http://static.getclicky.com/26219.js" type="text/javascript"></script>
<!-- <noscript><p><img alt="Clicky" src="http://in.getclicky.com/26219-db2.gif" /></p></noscript> -->
</body>
</html>