%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
%

% 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.
%    

tic; %flops(0);
N=100;K=20;maxM=8;
O = [5000 63 zeros(1,K-2)];%O = zeros(1,K);
M=zeros(maxM,K);
for mode=1:maxM,
	fn=lower(char([ computer, version,'.log',int2str(mode)]));
        fid=fopen(fn,'w+');
	fn=lower(char([ computer, version,'.lx',int2str(mode)]));
        fid2=fopen(fn,'w+');
        fprintf(fid,'%s %s\n', computer, version);
        fprintf(fid2,'%s %s\n', computer, version);
        for k=1:K,
		tic;
                [XYZ,M(mode,k)]=fermat(N,k,mode);
		T(mode,k)=toc;
                if k>2,
                        fprintf(fid,'==== n=%i ==========\n',k);
                        fprintf(fid2,'==== n=%i ==========\n',k);
                        if ~isempty(XYZ),
				for k2 = 1:size(XYZ,1);
		                        fprintf(fid,'%i\t%i\t%i\n',XYZ(k2,1),XYZ(k2,2),XYZ(k2,3));
				end;	
                                %[mode,k],
                                tmp=log2(XYZ(:,3)).*k<=53;
                                tmp=XYZ(find(tmp),:);
				for k2 = 1:size(tmp,1);
	                                fprintf(fid2,'%i\t%i\t%i\n',tmp(k2,1),tmp(k2,2),tmp(k2,3));
				end;	
                        end;
                end;
        end;
        fclose(fid);
        fclose(fid2);
end;
VER=[version '           '];
%computer
%ret_val=[toc; flops];
ret_val = toc; 

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=lower(char([ computer, version,'.txt']));
fid=fopen(fn,'w+');
%fid=fopen('flt_rep.txt','w+');
fprintf(fid,'FLTTEST Report \t%s\n',date);
fclose(fid);

if VER(1:6)=='MIDEVA'
        fid=fopen(fn,'a');
        fprintf(fid,'\n %s\t%s \n ',computer,version);
        fprintf(fid,'\n %bx \t%bx \t%bx\n',realmax,realmin,eps);
elseif exist('OCTAVE_VERSION')        % OCTAVE
        eval(['!uname -a >> "',fn,'"'],'-1');
        fid=fopen(fn,'a');
        fprintf(fid,'%s\t%s\tOCTAVE\n',computer,version);
        fprintf(fid,'\n');
elseif strcmp(version,'3.5')
	%% FreeMat
        fid=fopen(fn,'a');
        fprintf(fid,'\n %s\t%s \n ',computer,version);
        fprintf(fid,'\n %x \tx \t%x\n',realmax,realmin,eps);
else                            % MATLAB
        eval(['!uname -a >> "',fn,'"'],'-1');
        fid=fopen(fn,'a');
        fprintf(fid,'%s\t%s\tMATLAB\n',computer,version);
        fprintf(fid,'%bx\t%bx\t%bx\n',realmax,realmin,eps);
end; 

fprintf(fid,'%60.50e\n',realmax);
fprintf(fid,'%60.50e\n',realmin);
fprintf(fid,'%60.50e\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');
fprintf(fid,'\n\nS2\t');
fprintf(fid,'%i\t',E(2,:));
fprintf(fid,'\n\nS10\t');
fprintf(fid,'%i\t',sum(E(1:10,:)));
fprintf(fid,'\n\nS20\t');
fprintf(fid,'%i\t',sum(E(1:20,:)));
fclose(fid);

fprintf(1,'If the report of Your platform is not listed at\n'); 
fprintf(1,'\thttp://pub.ist.ac.at/~schloegl/matlab/fermat\n'); 
fprintf(1,'send the file %s to alois.schloegl@ist.ac.at\n',fn); 
fprintf(1,'\n',fn); 
type(fn);