function [XYZ,M,M2]=fermat(N,k,Mode);
%function XYZ=fermat(N,k, Mode);
% Interesting is that this programm behaves differently 
% on Matlab5.1 and Matlab4.2.  Try fermat(100,2,1);
%
% Fermats Last Theorem (FLT) states that for k>0 no 
% non-zero integer  numbers X, Y, Z exist such that
%      X^k+Y^k=Z^k 
% fermat(N,k) is an algorithm aimed for calculating 
% all FLT tuples [X Y Z] with X,Y <= N.
% default k=2, N=100, Mode=0;
% [X Y Z] and [Y X Z] are included separably in XYZ.
% Therefore, XYZ consists of twice the number of tuples.
% XYZ is sorted
% XYZ(XYZ(:,1)>XYZ(:,2),:) give the unique tuples
%
% 0<Mode<5 is used for testing the accuracy of floating point operations 
% The following table shows the result for fermat(100,2,1)
% version platform  #tuples    Algorithm  
% 4.2    PCWIN        63      fermat(100,2,1)
% 4.2    DecAlpha     63      fermat(100,2,1)
% 5.1    PCWIN        20      fermat(100,2,1)
%
% 4.2    PCWIN        63      fermat(100,2,2)
% 4.2    DecAlpha     63      fermat(100,2,2)
% 5.1    PCWIN        25      fermat(100,2,2)
%
% 2.0.10 Octave/Linux 63      fermat(100,2,2)

%	Version 1.06
%	Copyright (c) 1998-2002,2006,2012 by Alois Schloegl
%	e-mail: alois.schloegl@ist.ac.at	
%       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.
%    


if nargin<1 N=100; end;
if nargin<2 k=2; end;
if nargin<3 Mode=0; end;

% Calculate all Combinations, b is a matrix, the row and cols indicate X and Y
a=1:N;
b=a(ones(N,1),:).^k;
b=b+b';

%find the fermat numbers; calculate the k-th root and test on integer.
if Mode==0
        if k==2 
                b=sqrt(b); [y,x]=find(b==fix(b)); 
        else    % assuming the proof of FLT is correct
                x=[];y=[];    
        end
elseif Mode==1
        b=b.^(1/k); [y,x]=find(b==fix(b)); 
elseif Mode==2
        [y,x]=find(b==(fix(b.^(1/k)).^k));  
elseif Mode==3
        b=b.^(1/k); [y,x]=find(abs(b - round(b)) <= 4*eps*b);
elseif Mode==4
        b=b.^(1/k); [y,x]=find(abs(b - round(b)) <= 4*eps);
elseif Mode==5
        b=exp(log(b)./k); [y,x]=find(b==fix(b)); 
elseif Mode==6
        b=2.^(log2(b)./k); [y,x]=find(b==fix(b)); 
elseif Mode==7
        [y,x]=find(b==(exp(log(fix(exp(log(b)./k))).*k)));  
elseif Mode==8
        [y,x]=find(b==(2.^(log2(fix(2.^(log2(b)./k))).*k)));  
end;   

%construct Fermat-tuple
XYZ=[x,y,(x.^k+y.^k).^(1/k)];

%output
M = length(x)/2;
if ~isempty(XYZ),
	M2 = sum(log2(XYZ(:,3).^k)<=53)/2;
else
	M2 = 0;        
end;
%fprintf(1,'Number of different FLT-tuples found: %i\t(%i)\n',M,M2);

%return;
