//Program weron.cpp to generate levy distribution
//following algorithm of Weron  #{2864}
//October 2 1999
#include <iostream.h>
#include <stdio.h>
#include<stdlib.h>
#include<math.h>
long idum;
double levy(int flag);
double levy_alpha;		//exponent of distribution
double levy_beta;			//skewness parameter of distribution
double levy_bab;			//B(alpha,beta)
double levy_sab;			//S(alpha,beta)
double levy_pi2=2.*atan(1.);  //pi/2
double ran1(void);
double main(void)
	{
	static char file[30];   			//name of output file
	double flag=1;	//flag=1 if alpha=1
	FILE *fp;
	double i;
	double ndata;	//number of data points
	double small=.001;
	double x;		//levy distributed random variable
	cout<<"input rng seed ";
	cin >> idum;
	cout<<"exponent 0<alpha<=2 ";
	cin>> levy_alpha;
	cout<<"-1<=beta<=1 ";
	cin>>levy_beta;
	if (abs(levy_alpha-1.)-small)
		{
		flag=0;
		double y=levy_beta*tan(levy_pi2*levy_alpha);
		levy_bab=atan(y)/levy_alpha;
		levy_sab=pow(1+y*y,1/(2.*levy_alpha));
		}
	cout<<"filename for output: ";
	cout.flush();
	scanf("%s",&file);
	fp=fopen(file,"w");
	cout<<"number of data points ";
	cin>> ndata;
	fprintf(fp,"n=%lf, alpha= %lf beta=%lf \n", ndata, levy_alpha,levy_beta);
	for(i=1;i<=ndata;i++)
		{
		x=levy(flag);
		fprintf(fp,"%lf   %lf\n",i,x);
		}
	fclose(fp);
	return 0;
	}
double levy(int flag)
	{
	double v;		//uniform -pi/2,pi/2 distributed variable
	double w;		//exponentially distributed variable mean 1
	double x,y;    	//help variables
	v=levy_pi2*(2.*ran1()-1);
	w=-log(ran1());
	if(flag)
		{
		y=levy_pi2+levy_beta*v;
		x=(y*tan(v)-levy_beta*log(levy_pi2*w*cos(v)/y))/levy_pi2;
		}
	else
		{
		y=levy_alpha*(v+levy_bab);
		x=levy_sab*sin(y)/pow(cos(v),1/levy_alpha)*pow((cos(v-y)/w),
		1-1/levy_alpha);
		}
	return x;
	}
//a positive global variable long idum must be declared
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
double ran1()
{
	int j;
	long k;
	static long iy=0;
	static long iv[NTAB];
	double temp;

	if (!iy) 
		{
		for (j=NTAB+7;j>=0;j--) 
			{
			k=idum/IQ;
			idum=IA*(idum-k*IQ)-IR*k;
			if (idum < 0) idum += IM;
			if (j < NTAB) iv[j] = idum;
		        }
		iy=iv[0];
		}
	k=idum/IQ;
	idum=IA*(idum-k*IQ)-IR*k;
	if (idum < 0) idum += IM;
	j=iy/NDIV;
	iy=iv[j];
	iv[j] = idum;
	if ((temp=AM*iy) > RNMX) return RNMX;
	else return temp;
}
#undef IA
#undef IM
#undef AM
#undef IQ
#undef IR
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX