function z = cplxpair_Mod (x, tol, dim)
%
% Modified Code wrt Release : R14 (now made suitable for R13 and R12 also)
% Modified portions by : Sundar Krishnan (sundark100@yahoo.com)
% June-July 2005
%
% Modifications' Description :
% ----------------------------
% cplxpair_Mod is a modified version of the Standard R14 Matlab's cplxpair.
% The modification was found necessary to tackle genuine cases of
% Conjugate-Paired roots of Polynomials having a reasonably high degree,
% or Polynomials with multiple roots.
% The errors creep in from Matlab's standard roots and conv.
% With this modification, I hope that we can avoid getting unncecessary
% Error Messages for cplxpair operation of roots of such
% genuine cases of Polynomials.
% Essentially, I have differentiated between tol which is used for
% identifying (almost) Pure Real Nos and a newly introduced tol_Real
% which is used to identify if the pairs are Conjuagetes.
% I have adapted both tol and tol_Real by increasing the respective
% tolerances as the degree of the Polynomial increases.
% Also, I have assigned explicit Error Messages.
% The rest of the code essentially remains intact as it exists
% in Matlab's R14.
%
% &&&&&&&&&&
%
% A more comprehensive way to deal with Complex and Real Roots
% than cplxpair or cplxpair_Mod, is to instead use my Programme
% cmplx_roots_sort_Det.m.
% This Programme counts the number of Real and Complex Roots,
% and also calculates the multiplicity of each root ;
% cmplx_roots_sort_Det.m can also deal with cases where
% the Complex Roots may not be in Conjugate Pairs.
% cmplx_roots_sort_Det.m is very useful when we need to calculate
% Fractional Powers of Polynomials, as I have done in Poly_POWER.m
%
% cmplx_roots_sort_Det.m makes two calls to cplxpair_Mod.
% One call is via cmplx_roots_sort.m, and one is a direct call,
% both of which involve calls with inputs as Complex Roots only.
%
% If this suite of Programmes (in this zip file) don't work for you
% directly when downloaded into an independent directory, you could
% try loading these files in the Matlab dir : ...\toolbox\matlab\specfun
%
% &&&&&&&&&&&&&&&&&
%
% Author's Legal Disclaimer :
% ---------------------------
% I share this modified code, developed with lots of efforts,
% in good faith.
% Everyone is free to use this code at their own risk provided
% they include my name for this part in their works.
% I am NOT responsible for any technical or legal complications
% directly or indirectly, arising out of any direct or indirect use
% of my codes or my interpretations, including my modified codes.
% Any apparent TECHNICAL fault in interpreting any external references,
% or others' replies in any Public User Group Forums, is my own,
% not theirs. However, I am NOT LEGALLY responsible for any of my
% technical or other interpretation of others' comments in ANY Forum.
%
% Since this is a modified code, I have retained the original initial
% comments of Matlab (see below after Usage Egs).
%
% &&&&&&&&&&&&&&&&&
%
% % Usage Egs for testing the Modifications :
% % Case 1 :
% clc
% Poly_Rts_Cmplx = [ 3 + 4.02j, 4 + 3j, 4 - 3j, ...
% 12 - 1j, 12 + 1j, 3 - 4.02j, ...
% 6 + 8j, 6 - 8j, 6 + 8.05j, 6 - 8.05j, 6 - 8j, ...
% 6 + 8j, 3 - 4.020003j, 3 + 4.020003j, ...
% 1 - 12j, 1 + 12j, 8 + 6j, 8 - 6j ] ;
%
% len_Poly_Rts_Cmplx = length ( Poly_Rts_Cmplx ) ;
%
% Poly_T_1 = [1] ;
% for i = 1 : 14 % len_Poly_Rts_Cmplx
% Poly_T_1 = conv ( Poly_T_1, [1, -Poly_Rts_Cmplx(i)] ) ;
% end
% Poly_T_1 % deg = 14
% roots_Poly_T_1 = roots ( Poly_T_1 )
%
% Poly_T_2 = Poly_POWER (Poly_T_1, 2) % deg = 28
% roots_Poly_T_2 = roots ( Poly_T_2 )
%
% cplxpair_Mod (roots_Poly_T_1)
% cplxpair_Mod (roots_Poly_T_2)
%
%
% Poly_T_3 = [1] ;
% for i = 1 : 16 % len_Poly_Rts_Cmplx
% Poly_T_3 = conv ( Poly_T_3, [1, -Poly_Rts_Cmplx(i)] ) ;
% end
% Poly_T_3 % deg = 16
% roots_Poly_T_3 = roots ( Poly_T_3 )
%
% Poly_T_4 = Poly_POWER (Poly_T_3, 2) % deg = 32
% roots_Poly_T_4 = roots ( Poly_T_4 )
%
% cplxpair_Mod (roots_Poly_T_3)
% cplxpair_Mod (roots_Poly_T_4)
%
% % &&&&&&&&&&
%
% % Case 2 :
% % Repeat above with lower magnitudes :
% clc
% Poly_Rts_Cmplx = 0.1 * Poly_Rts_Cmplx ;
%
% % Proceed as above for tests at 4 degrees.
%
% % &&&&&&&&&&
%
% % Case 3 :
% clc
% Poly_Rts_Cmplx = [ -13 + 4.02j, 4 + 13j, 4 - 13j, ...
% 12 - 1j, 12 + 1j, -13 - 4.02j, ...
% 6 + 8j, 6 - 8j, 6 + 18.05j, 6 - 18.05j, 6 - 8j, ...
% 6 + 8j, 13 - 4.020003j, 13 + 4.020003j, ...
% 1 - 12j, 1 + 12j, 8 + 6j, 8 - 6j ] ;
% % Proceed as above for tests at 4 degrees.
%
% % &&&&&&&&&&
%
% % Case 4 : Case 3 Plus Multiple Real Roots
% % This will "fail" because the Std Matlab counts only by Real Parts !
% clc
% Poly_Rts_Cmplx = [ -6, -6, 7, 7, -6, 8, ...
% -13 + 4.02j, 4 + 13j, 4 - 13j, ...
% 12 - 1j, 12 + 1j, -13 - 4.02j, ...
% 6 + 8j, 6 - 8j, 6 + 18.05j, 6 - 18.05j, 6 - 8j, ...
% 6 + 8j, 13 - 4.020003j, 13 + 4.020003j, ...
% 1 - 12j, 1 + 12j, 8 + 6j, 8 - 6j ] ;
% % Proceed as above for tests at 4 degrees.
%
% % &&&&&&&&&&&&&&&&&
%
% Matlab's Original Initial Comments :
% CPLXPAIR Sort numbers into complex conjugate pairs.
% Y = CPLXPAIR(X) takes a vector of complex conjugate pairs and/or
% real numbers. CPLXPAIR rearranges the elements of X so that
% complex numbers are collected into matched pairs of complex
% conjugates. The pairs are ordered by increasing real part.
% Any purely real elements are placed after all the complex pairs.
% Y = CPLXPAIR(X,TOL) uses a relative tolerance of TOL for
% comparison purposes. The default is TOL = 100*EPS.
%
% For X an N-D array, CPLXPAIR(X) and CPLXPAIR(X,TOL) rearranges
% the elements along the first non-singleton dimension of X.
% CPLXPAIR(X,[],DIM) and CPLXPAIR(X,TOL,DIM) sorts X along
% dimension DIM.
% L. Shure 1-27-88
% Revised D. Orofino 4-30-96
% $Revision: 5.14.4.1 $ $Date: 2003/05/01 20:41:23 $
% *******************
% 1) Inits :
if isempty(x)
z = x ;
return % Quick exit if input is empty
end
if nargin == 3
nshifts = 0 ;
perm = [ dim : max(ndims(x),dim) 1 : dim-1 ] ;
% Note that it is something like having 2 horizontal sections above ! :
% Eg 1 : x = [ 25 : 12 9 : 14 ]
% x = [ 9 10 11 12 13 14 ]
% Eg 2 : x = [ 25 : 27 9 : 10 ]
% x = [ 25 26 27 9 10 ]
x = permute(x, perm) ;
else
[x, nshifts] = shiftdim(x) ;
perm = [ ] ;
end
% Supply initial values for tolerances :
% (SK : see my modifications further down.)
% Often the need is more stringent on tol_Real than on tol.
if nargin < 2 | isempty(tol) % | for R12
tol = 1e-8 ; % 100*eps for both tol and tol_Real is toooo small.
% Additional tol = tol_Real for Real and Imag Parts added by SK :
tol_Real = 9e-5 ;
end
% Reshape x to a 2-D matrix:
xsiz = size(x) ; % original shape of input
x = x(:, :) ; % reshape to a 2-D matrix
z = zeros(size(x)) ; % preallocate temp storage
errmsg = 'Complex numbers can''t be paired.' ;
errmsg_prefix = 'MATLAB:cplxpair_Mod:ComplexValuesPaired' ; % by SK
% No space as in : 'MATLAB : cplxpair_Mod : ComplexValuesPaired' ;
% See Message Identifiers in Help.
% Foll by SK :
[r_x, c_x] = size(x) ;
deg = max (r_x, c_x) ;
% Usually, we might give a row vector as input.
% But it should work even if the input is a col vector.
% SK : June 2005 :
% If Matlab's default tol of 100 * eps is retained for both
% tol and tol_Real, we would not be able to make any serious headway
% with the results obtained from Matlab's standard roots and conv !
%
% So, tol and tol_Real are modified below as a function of the degree
% of the Polynomial ie, allow for more tol as the degree goes higher.
%
% Often the need is more stringent on tol_Real than on tol.
% I have modified tol proportional to the degree of the Polynomial,
% but the values of tol below are however, still significantly lower than
% pure_real_fctr that I have used in cmplx_roots_sort_Det.m
% to eliminate the "superfluous" negligible imag parts created
% by Matlab's roots.
% The main reason is that I did not want to trample the original code
% drastically ; secondly, tol does not affect much in the working of
% cmplx_roots_sort_Det.m or cmplx_roots_sort.m
% So, the real tough task is that of adjusting tol_Real.
if nargin < 2 | isempty(tol) % | for R12
if deg > 6
% ie, for deg = 1 to 6 :
% tol = 1e-8 ; tol_Real = 9e-5
if deg <= 10 % ie, for deg = 7, 8, 9, 10
tol = 5e-7 ;
tol_Real = 9e-5 * (deg/6)^0.25 ; % 9.35e-5 to 1.02e-4
elseif deg <= 15
tol = 3e-6 ;
tol_Real = 1.5e-4 * (deg/10)^3 ; % 2e-4 to 5.06e-4
elseif deg <= 20
tol = 5e-6 * (deg /16)^1.5 ; % 5e-6 to 6.99e-6
tol_Real = 6e-4 * (deg/16)^3 ; % 6e-4 to 1.17e-3
elseif deg <= 25
tol = 1e-5 * (deg /21)^12 ; % 1e-5 to 8.1e-5
tol_Real = 1.5e-3 * (deg/21)^14 ; % 0.0015 to 0.0172
elseif deg <= 30
tol = 1e-4 * (deg /26)^24 ; % 1e-4 to 3.1e-3
tol_Real = 0.021 * (deg/26)^18 ; % 0.021 to 0.276
elseif deg <= 35
tol = 4e-3 * (deg /31)^8 ; % 4e-3 to 1.06e-2
tol_Real = 0.35 * (deg/31)^9 ; % 0.35 to 1.04
else
tol = 0.015 * (deg /36)^5 ; % 0.015 to 0.025 (40)
tol_Real = 1.5 * (deg/36)^6 ; % 1.5 to 2.82 (40)
end
end % if deg > 6
end % if nargin < 2 || isempty(tol)
% SK : The foll (commented code) is just an Average Check.
% In the Main Code below, the abs is taken for EACH element.
% deg, tol, tol_Real
% mean_abs = mean ( mean ( abs( x ) ) )
% abs_tol = tol * mean_abs
%
% pause
% &&&&&&&&&&&&&&&&&
% 2) The original cplxpair code with more clear Error Messages
for k = 1 : size(x, 2)
% Get next column of x:
xc = x(:, k) ;
% Find purely-real entries:
idx = find ( abs(imag(xc)) <= tol * abs(xc) ) ;
nr = length(idx) ; % Number of purely-real entries
if ~isempty(idx)
% Store sorted real's at end of column:
z(end-nr+1 : end, k) = sort( real(xc(idx)) ) ;
xc(idx) = [ ] ; % Remove values from current column
end
nc = length(xc) ; % Number of entries remaining in input column
if nc > 0
% Complex values in list:
if rem(nc, 2)==1
% Odd number of entries remaining
fprintf('\nProblem in cplxpair_Mod : wrt Odd no of entries\n');
x, deg, tol, tol_Real
errmsg_1 = strcat ( errmsg, ...
' : Odd number of entries remaining' ) ;
error ( errmsg_prefix, errmsg_1 ) ;
end
% Sort complex column-vector xc, based on its real part:
[xtemp, idx] = sort ( real(xc) ) ;
xc = xc(idx) ; % Sort complex numbers based on real part
% Check if real parts occur in pairs:
% [ SK : June 2005 :
% This is often the culprit code if Matlab's default tol is used!
% when the degree of the Polynomial is high like for eg, > 15. ]
% Compare to xc() so imag part is considered
% (in case real part is nearly 0).
% Arbitrary choice of using abs(xc(1:2:nc))
% or abs(xc(2:2:nc)) for tolerance
if any ( abs ( xtemp(1 : 2 : nc) - xtemp(2 : 2 : nc) ) > ...
tol_Real .* abs(xc(1 : 2 : nc)) )
fprintf('\nProblem in cplxpair_Mod : wrt Real parts Test\n');
x, deg, tol, tol_Real
errmsg_2 = strcat ( errmsg, ' : Real parts not in pairs' ) ;
error ( errmsg_prefix, errmsg_2 ) ;
end
% Check real part pairs to see if imag parts are conjugates:
nxt_row = 1 ; % next row in z(:, k) for results
while ~isempty(xc)
% Find all real parts "identical" to real(xc(1)):
idx = find ( abs ( real(xc) - real(xc(1)) ) ...
<= tol_Real .* abs(xc) ) ;
nn = length(idx) ; % # of values with identical real parts
if nn <= 1
% Only 1 value found - certainly not a pair!
fprintf('\nProblem in cplxpair_Mod : Imag Parts Test 1\n');
x, deg, tol, tol_Real
errmsg_3 = strcat (errmsg, ...
' : Imag parts not conjugates (Test 1)') ;
error ( errmsg_prefix, errmsg_3 ) ;
end
% There could be multiple pairs with "identical" real parts.
% Sort the imaginary parts of those values with identical real
% parts - these SHOULD be the next N entries, with N even.
[xtemp, idx] = sort ( imag(xc(idx)) ) ;
xq = xc(idx) ; % Get complex-values with identical real parts,
% which are now sorted by imaginary component.
% Verify conjugate-pairing of imaginary parts :
if any ( abs ( xtemp + xtemp(nn : -1 : 1) ) ...
> tol_Real .* abs(xq) )
fprintf('\nProblem in cplxpair_Mod : Imag Parts Test 2\n');
x, deg, tol, tol_Real, xtemp
errmsg_4 = strcat (errmsg, ...
' : Imag parts not conjugates (Test 2)') ;
error ( errmsg_prefix, errmsg_4 ) ;
end
% Keep value with pos imag part,
% and compute conjugate for pair.
% List value with most-neg imag first, then its conjugate.
z(nxt_row : nxt_row+nn-1, k) = reshape ( ...
[ conj( xq(end : -1 : nn/2+1) ), ...
xq(end : -1 : nn/2+1) ].' , ...
nn, 1 ) ;
nxt_row = nxt_row + nn ; % Bump next-row pointer
xc(idx) = [ ] ; % Remove entries from xc
end % while ~isempty(xc)
end % of complex-values check : if nc > 0
end % of column loop : for k = 1 : size(x, 2)
% &&&&&&&&&&&&&&&&&
% 3) Reshape Z to appropriate form
z = reshape(z, xsiz) ;
if ~isempty(perm)
z = ipermute(z, perm) ; % Inverse permute
end
if nshifts ~= 0
z = shiftdim(z, -nshifts) ;
end
% *******************
