//Utilities. Custom Functions.
//Erfh. We were not able to run Apache... erf reliably, we have to
// rapidly implement approximate erf.
// Usage must be strictly prepared before use.
// For example, put prepareErfc in Applet init code.
public class UCFun {
public static double erfh(double x, boolean HeavisideMode){
if(HeavisideMode)return x>0 ? 0.0 :1.0;
return erfc(x);
}
//=========================================================================
//Brute force (asymmetric) erfc calculation.
//-------------------------------------------------------------------------
//Enable this function if you don't have Internet and
//don't have decent Java Math. package.
public static double erfc(double x){
if(x<0)return 1.0-erfc(-x);
double fraction=x%step;
int left=Math.max(0,(int)((x-fraction)*step1));
int right=left+1;
if(right>=accur)return 0.5;
//Interpolate:
return (mErfc[right]-mErfc[left])*(x*step1-left)+mErfc[left];
}
private static double[] mErfc;
private static int accur;
private static double norm;
private static double step;
private static double step1;
//Input: Combination 1000000, 100 of accuracy,subaccuracy
// gave an overall integration accuracy 2*10^(-7)
// (About 4 seconds of preparation on 2GHz PC).
//
public static void prepareErfc(int accuracy,int subaccuracy){
accur=accuracy;
mErfc=new double[accur];
norm=1.0/Math.sqrt(2.0*Math.PI);
step=10.0/accur;
step1=1.0/step;
int lim=accur-1;
double factor=step*norm;
mErfc[0]=0.5;
double sum=0;
mErfc[0]=0.5;
for(int i=0; i<lim; i++){
double x=i*step;
//We set beginning point of an interval.
//Now, do integrate carefully through this interval:
double substep=step/subaccuracy;
double subfactor=factor/subaccuracy;
double subx=x;
double subsum=0.0;
for(int ii=0; ii<subaccuracy; ii++){
subx+=substep;
subsum+=Math.exp(-(subx*subx*0.5))*subfactor;
}
sum+=subsum;
mErfc[i+1]=0.5-sum;
//UTil.con("i"+i+" x"+x+" sum="+sum);
}
UTil.con(" erfc accuracy reached="+mErfc[lim-1]);
for(int k=lim-1; k>-1; k--){
if(mErfc[k]<0){
mErfc[k]=0.0; //Fix accuracy defect.
}else{
UTil.con(""+k+"=k, point number where integration overruns barrier of 0.5.");
break;
}
}
}
//-------------------------------------------------------------------------
//Brute force (asymmetric) erfc calculation.
//=========================================================================
public static double erfIntegrand(double x){
return (erfc(x)-erfc(x+step))*step1;
}
public static double erfDerivative(double x){
return Math.exp(-(x*x*0.5))*norm;
}
}
Copyright (C) 2009 Konstantin Kirillov