%FLT Test - Calculates the accuracy of floating point operations
%   based on calculating Numbers of Fermat's Last Theorem 
%   More information is available at
%       http://pub.ist.ac.at/~schloegl/matlab/fermat

%	Copyright (c) 1998,2006,2012 by Alois Schloegl
%	e-mail: alois.schloegl@gmail.com	
%       http://pub.ist.ac.at/~schloegl/matlab/fermat
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the  License, or (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

N = 1000;K=100;maxM=8;
O = [N*N/2 1034 zeros(1,K-2)];
M = zeros(maxM,K);
%for k=1:K,
for mode=1:maxM,
        %[XYZ,M(mode,k)]=fermat(N,k,mode);
	fn=lower(char([ computer, version,'.lx100',int2str(mode)]));
        fid=fopen(fn,'w+');
        fprintf(fid,'%s %s\n', computer, version);
        for k=1:K,
                [XYZ,M(mode,k)]=fermat(N,k,mode);
                if k>2,
                        fprintf(fid,'==== n=%i ==========\n',k);
                        if ~isempty(XYZ),
                                %[mode,k],
                                tmp=log2(XYZ(:,3).^k)<53;
                                tmp=XYZ(find(tmp),:);
                                fprintf(fid,'%i\t%i\t%i\n',tmp');
                        end;
                end;
        end;
        fclose(fid);
end;
%end;
VER = version;

E = abs(M-O(ones(maxM,1),:))';
%E=[E;sum(E)];

fprintf(1,'\nFTL benchmark of floating point accuracy: %i\n',sum(E));
fprintf(1,'The smaller the more accurate is the calculation\n'); 

fn=char([computer, version,'.txt1000']);
%fid=fopen('flt_rep.txt','w+');
fid=fopen(fn,'w+');
fprintf(fid,'FLTTEST%i Report \t%s\n',N,date);

if exist('OCTAVE_VERSION')        % OCTAVE
        fprintf(fid,'%s\t%s\tOCTAVE\n',computer,version);
        fprintf(fid,'\n');
else                            % MATLAB
        fprintf(fid,'%s\t%s\tMATLAB\n',computer,version);
        fprintf(fid,'%bx\t%bx\t%bx\n',realmax,realmin,eps);
end; 
fprintf(fid,'%80.70e\n',realmax);
fprintf(fid,'%80.70e\n',realmin);
fprintf(fid,'%80.70e\n',eps);

fprintf(fid,'E');
fprintf(fid,'\t%i\t%i\t%i\t%i\t%i\t%i\t%i\t%i\n',E');

for k=[2 10:10:K],
	fprintf(fid,'\n\nS%i\t',k);
	fprintf(fid,'%i\t',sum(E(1:k,:)));
end;
fclose(fid);

fprintf(1,'If the report of Your platform is not listed at\n'); 
fprintf(1,'\thttp://www.dpmi.tu-graz.ac.at/~schloegl/matlab/fermat\n'); 
fprintf(1,'send the file %s to a.schloegl@ieee.org\n',fn); 
type(fn);
