%
% MATLAB code to generate single precision approximations for erfinv
% following procedure described in GPU Gems article
%

close all
clear all

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

fig = 0;

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

  for pass = 1:min(2,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 = 8;
        wlo    = 0;
        whi    = 5;
        wmd    = 2.5;
      else           % tail region
        degree = 8;
        wlo    = sqrt(5);
        whi    = 4;
        wmd    = 3;
      end

    elseif (alt==2)
      if (pass==1)   % central region
        degree = 9;
        wlo    = 0;
        whi    = 6.25;
        wmd    = 3.125;
      else           % tail region
        degree = 7;
        wlo    = 2.5;
        whi    = 4;
        wmd    = 3.25;
      end

    else
      degree = 14;
      wlo    = 0;
      whi    = 16;
      wmd    = 8;
    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));
    else
      x = sqrt(1-exp(-w.^2));
    end

    y = erfinv(x);
    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))';

%
% plot results
%

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

    p1 = single(p1);
    fig = fig+1;  figure(fig);
    plot(w,(single(u)-polyval(p1,(w-wmd)/scal))./eps(single(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(' %12.9g',scal^(deg-degree)*p1(deg+1)))
%    end

%
% output to file
%

    if (pass==1)
      file = strvcat(' ',...
                 '__inline__ __device__ float MBG_erfinv(float x) ',...
                 '{                                               ',...
                 '  float w, p;                                   ',...
                 '                                                ',...
                 '  w = - __logf((1.0f-x)*(1.0f+x));              ',...
                 '                                                ');

       if (alt==3)
         file = strvcat(file,...
         sprintf('  w = w - %ff;                                  ',wmd) );
       else
         file = strvcat(file,...
         sprintf('  if ( w < %ff ) {                              ',whi),...
         sprintf('    w = w - %ff;                                ',wmd) );
       end
    else
      file = strvcat(file,'  else {',...
         sprintf('    w = sqrtf(w) - %ff;',wmd) );
    end

    if (alt==3)
      file = strvcat(file,sprintf('  p = %16.9gf;',scal^(-degree)*p1(1)));
    else
      file = strvcat(file,sprintf('    p = %16.9gf;',scal^(-degree)*p1(1)));
    end
    
    for deg = 1:degree
      if (alt==3)
        file = strvcat(file,sprintf('  p = %16.9gf + p*w;',scal^(deg-degree)*p1(deg+1)));
      else
        file = strvcat(file,sprintf('    p = %16.9gf + p*w;',scal^(deg-degree)*p1(deg+1)));
      end
    end

    if (alt<3)
      file = strvcat(file,'  }');
    end

    if (pass==2 || alt==3)
      file = strvcat(file,' ','  return p*x;','}');

      fid = fopen(sprintf('erfinv_SP_%d.cu',alt),'wt');
      fprintf(fid,'// \n// auto-generated by erfinv_SP_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

