
close all
clear all

th  = pi*(0.5 + 0:199)'/200;

fig = 0;

for alt = 1:2  % two alternatives
  disp(sprintf('\n%%\n%% alternative %d:',alt))

  for pass = 1:4-alt
    if (pass==1)
      disp(sprintf('%%\n%% central region \n%%'))
    else
      disp(sprintf('%%\n%% tail region \n%%'))
    end

    if (alt==1)
      if (pass==1)   % central region
        degree = 22;
        wlo    = 0;
        whi    = 6.25;
        wmd    = 3.125;
      elseif (pass==2)   % tail region 1
        degree = 18;
        wlo    = 2.5;
        whi    = 4;
        wmd    = 3.25;
      else               % tail region 2
        degree = 16;
        wlo    = 4;
        whi    = 6;
        wmd    = 5;
      end

    else
      if (pass==1)   % central region
        degree = 22;
        wlo    = 0;
        whi    = 6.25;
        wmd    = 3.125;
      else           % tail region
        degree = 26;
        wlo    = 2.5;
        whi    = 6;
        wmd    = 4.25;
      end
    end

    scal = 0.5*(whi-wlo);

    w = 0.5*(whi+wlo) - 0.5*(whi-wlo)*cos(th);
    if (pass==1)
      x = sqrt(1-exp(-w));
      y = erfcinv(exp(-w)./(1+sqrt(1-exp(-w))));
    else
      x = sqrt(1-exp(-w.^2));
      y = erfcinv(exp(-w.^2)./(1+sqrt(1-exp(-w.^2))));
    end
    u = y ./ x;

%
% diagonal scaling improves least squares matrix condition number
%

    D = diag(1./u);
    A = zeros(length(u),degree+1);
    for deg = 0:degree
      A(:,degree+1-deg) = ((w-wmd)/scal).^deg;
    end

    p1 = ((D*A)\(D*u))';

%
% error correction, using vpa (variable precision arithmetic)
%                   which needs symbolic toolbox
%
    ye = vpa(y);
    xe = erf(ye);
    disp(' computed xe')
    if (pass==1)
      we = -log(1-xe.^2);
    else
      we = sqrt(-log(1-xe.^2));
    end
    ue = ye./xe;

    p1e = vpa(p1);

    ze  = p1e(1);
    for i = 1:degree;
      ze = ze.*(we-wmd)/scal + p1e(i+1);
    end

    err = double(ue-ze);

    p1e = p1e + vpa( ((D*A)\(D*err))');
    p1  = double(p1e);

    ze  = p1e(1);
    for i = 1:degree;
      ze = ze.*(we-wmd)/scal + p1e(i+1);
    end
    err = ue - ze;

    fig = fig+1;  figure(fig);
    plot(w,err./u);
    title(sprintf('Alternative %d',alt)); ylabel('relative error')
    if (pass==1)
      xlabel('w')
    else
      xlabel('sqrt(w)')
    end

    p1 = double(p1);
    fig = fig+1;  figure(fig);
    plot(w,(u-polyval(p1,(w-wmd)/scal))./eps(u));
    title(sprintf('Alternative %d',alt)); ylabel('ulp error')
    if (pass==1)
      xlabel('w')
    else
      xlabel('sqrt(w)')
    end

%    disp(sprintf(' mid = %12.9g \n p%d  = ',wmd,pass))
%
%    for deg = 0:degree
%      disp(sprintf(' %20.18g',scal^(deg-degree)*p1(deg+1)))
%    end

%
% output to file
%

    if (pass==1)
      file = strvcat(' ',...
                 '__inline__ __device__ double MBG_erfinv(double x) ',...
                 '{                                                ',...
                 '  double w, p;                                   ',...
                 '                                                 ',...
                 '  w = - log((1.0-x)*(1.0+x));                    ',...
                 '                                                 ',...
         sprintf('  if ( w < %f ) {                              ',whi),...
         sprintf('    w = w - %f;                                ',wmd) );

    elseif(alt==1 & pass==2)
      file = strvcat(file,...
         sprintf('  else if ( w < %f ) {',whi^2),...
         sprintf('    w = sqrt(w) - %f;',wmd) );

    else
      file = strvcat(file,'  else {',...
         sprintf('    w = sqrt(w) - %f;',wmd) );
    end

    file = strvcat(file,sprintf('    p = %27.20g;',scal^(-degree)*p1(1)));
    
    for deg = 1:degree
      file = strvcat(file,sprintf('    p = %27.20g + p*w;',scal^(deg-degree)*p1(deg+1)));
    end

    file = strvcat(file,'  }');

    if (pass==4-alt)
      file = strvcat(file,' ','  return p*x;','}');

      fid = fopen(sprintf('erfinv_DP_%d.cu',alt),'wt');
      fprintf(fid,'// \n// auto-generated by erfinv_DP_approx.m on %s \n//\n\n',datestr(now));
      for n=1:size(file,1)
        fprintf(fid,'%s\n',file(n,:));
      end
      fclose(fid);
    end

  end

end

