function hdmrd(K)
% hdmrd(K) demonstrates how Interval Arithmetic's error-bounds
% grow exponentially rather than linearly, as they should, for
% N repeated multiplications by an orthogonal matrix M obtained
% from Matlab's hadamard(2*K) whose every element is +1 or -1 ;
% M = hadamard(2*K)/(2^K) == M' and M*M == eye(4^K) exactly.
% After interval vector x is initialized, no roundoff occurs;
% but M*M*...*M*x widens exponentially instead of not at all.
format compact, format short
n = 4^K ; disp([' Dimension n = 4^K = ', int2str(n)])
clck = clock ;
rand('seed', round(10^6*clck(6))) %... seed rand(..) @ seconds
x = rand ; %... throw away the first non-random rand(...)
x = 2048*rand(n,1) - 1024 ; % -1024 < random x < +1024
xhi = ceil(x) ; %... rounded UP to the next integer
xlo = floor(x) ; %... rounded DOWN to the next integer
disp( ' Initial interval vector [xlo, xhi]'' & its widths: ' )
disp( [xlo, xhi, xhi-xlo]' )
H = hadamard(2*K) ;
disp( ' Hadamard Matrix H = hadamard(2*K) = ... '), disp(H)
M = H*(0.5^K) ;
disp( [' M = H/2^K ; Norm(M - M'') = ', num2str(norm(M-M')) ] )
disp( [' Norm(M*M - eye(n)) = ', num2str(norm(M*M-eye(n))) ] )
Mp = M.*(M>0) ; Mn = M - Mp ; Mpn = [Mp, Mn] ;
for N = 1:4
disp(' ')
t = Mpn*[xlo; xhi] ; xhi = Mpn*[xhi; xlo] ; xlo = t ;
disp( ' new x = M*x ; [xlo, xhi, xhi-xlo]'' = ...' )
disp( [xlo, xhi, xhi-xlo]' )
end
format long e
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diary h
È clr
È hdmrd(1)
Dimension n = 4^K = 4
Initial interval vector [xlo, xhi]' & its widths:
-304 -989 -528 8
-303 -988 -527 9
1 1 1 1
Hadamard Matrix H = hadamard(2*K) = ...
1 1 1 1
1 -1 1 -1
1 1 -1 -1
1 -1 -1 1
M = H/2^K ; Norm(M - M') = 0
Norm(M*M - eye(n)) = 0
new x = M*x ; [xlo, xhi, xhi-xlo]' = ...
-906.5 73.5 -387.5 609.5
-904.5 75.5 -385.5 611.5
2.0 2.0 2.0 2.0
new x = M*x ; [xlo, xhi, xhi-xlo]' = ...
-305.5 -990.5 -529.5 6.5
-301.5 -986.5 -525.5 10.5
4.0 4.0 4.0 4.0
new x = M*x ; [xlo, xhi, xhi-xlo]' = ...
-909.5 70.5 -390.5 606.5
-901.5 78.5 -382.5 614.5
8.0 8.0 8.0 8.0
new x = M*x ; [xlo, xhi, xhi-xlo]' = ...
-311.5 -996.5 -535.5 0.5
-295.5 -980.5 -519.5 16.5
16.0 16.0 16.0 16.0
È hdmrd(1)
Dimension n = 4^K = 4
Initial interval vector [xlo, xhi]' & its widths:
780 232 -328 205
781 233 -327 206
1 1 1 1
Hadamard Matrix H = hadamard(2*K) = ...
1 1 1 1
1 -1 1 -1
1 1 -1 -1
1 -1 -1 1
M = H/2^K ; Norm(M - M') = 0
Norm(M*M - eye(n)) = 0
new x = M*x ; [xlo, xhi, xhi-xlo]' = ...
444.5 6.5 566.5 539.5
446.5 8.5 568.5 541.5
2.0 2.0 2.0 2.0
new x = M*x ; [xlo, xhi, xhi-xlo]' = ...
778.5 230.5 -329.5 203.5
782.5 234.5 -325.5 207.5
4.0 4.0 4.0 4.0
new x = M*x ; [xlo, xhi, xhi-xlo]' = ...
441.5 3.5 563.5 536.5
449.5 11.5 571.5 544.5
8.0 8.0 8.0 8.0
new x = M*x ; [xlo, xhi, xhi-xlo]' = ...
772.5 224.5 -335.5 197.5
788.5 240.5 -319.5 213.5
16.0 16.0 16.0 16.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
È hdmrd(2)
Dimension n = 4^K = 16
Initial interval vector [xlo, xhi]' & its widths:
Columns 1 through 12
-972 733 369 127 671 114 647 507 -202 924 -652 -102
-971 734 370 128 672 115 648 508 -201 925 651 -101
1 1 1 1 1 1 1 1 1 1 1 1
Columns 13 through 16
-461 -221 -500 -736
-460 -220 -499 -735
1 1 1 1
Hadamard Matrix H = hadamard(2*K) = ...
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1
1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1
1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1
1 1 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1
1 -1 1 -1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 1
1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1
1 -1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 1 1 -1
1 1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1
1 -1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1
1 1 -1 -1 1 1 -1 -1 -1 -1 1 1 -1 -1 1 1
1 -1 -1 1 1 -1 -1 1 -1 1 1 -1 -1 1 1 -1
1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1
1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1
1 1 -1 -1 -1 -1 1 1 -1 -1 1 1 1 1 -1 -1
1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 1 -1 -1 1
M = H/2^K ; Norm(M - M') = 0
Norm(M*M - eye(n)) = 0
new x = M*x ; [xlo, xhi, xhi-xlo]' = ...
Columns 1 through 7
61.5 -613.5 229.5 -647.5 49.0 -960.0 137.0
65.5 -609.5 233.5 -643.5 53.0 -956.0 141.0
4.0 4.0 4.0 4.0 4.0 4.0 4.0
Columns 8 through 14
-618.0 1034.5 226.5 -785.5 -121.5 -8.0 -124.0
-614.0 1038.5 230.5 -781.5 -117.5 -890.0 -120.0
4.0 4.0 4.0 4.0 4.0 4.0 4.0
Columns 15 through 16
-324.0 -568.0
-320.0 -564.0
4.0 4.0
new x = M*x ; [xlo, xhi, xhi-xlo]' = ...
Columns 1 through 7
-979.5 725.5 361.5 119.5 663.5 106.5 639.5
-963.5 741.5 377.5 135.5 679.5 122.5 655.5
16.0 16.0 16.0 16.0 16.0 16.0 16.0
Columns 8 through 14
499.5 -209.5 916.5 -659.5 -109.5 -468.5 -228.5
515.5 -193.5 932.5 -643.5 -93.5 -452.5 -212.5
16.0 16.0 16.0 16.0 16.0 16.0 16.0
Columns 15 through 16
-507.5 -743.5
-491.5 -727.5
16.0 16.0
new x = M*x ; [xlo, xhi, xhi-xlo]' = ...
Columns 1 through 7
31.5 -643.5 199.5 -677.5 19.0 -990.0 107.0
95.5 -579.5 263.5 -613.5 83.0 -926.0 171.0
64.0 64.0 64.0 64.0 64.0 64.0 64.0
Columns 8 through 14
-648.0 1004.5 196.5 -815.5 -151.5 -924.0 -154.0
-584.0 1068.5 260.5 -751.5 -87.5 -860.0 -90.0
64.0 64.0 64.0 64.0 64.0 64.0 64.0
Columns 15 through 16
-354.0 -598.0
-290.0 -534.0
64.0 64.0
new x = M*x ; [xlo, xhi, xhi-xlo]' = ...
Columns 1 through 7
-1099.5 605.5 241.5 -0.5 543.5 -13.5 519.5
-843.5 861.5 497.5 255.5 799.5 242.5 775.5
256.0 256.0 256.0 256.0 256.0 256.0 256.0
Columns 8 through 14
379.5 -329.5 796.5 -779.5 -229.5 -588.5 -348.5
635.5 -73.5 1052.5 -523.5 26.5 -332.5 -92.5
256.0 256.0 256.0 256.0 256.0 256.0 256.0
Columns 15 through 16
-627.5 -863.5
-371.5 -607.5
256.0 256.0
È hdmrd(2)
Dimension n = 4^K = 16
Initial interval vector [xlo, xhi]' & its widths:
Columns 1 through 12
45 509 -911 469 628 -973 -1000 255 -527 -14 -970 -142
46 510 -910 470 629 -972 -999 256 -526 -13 -969 -141
1 1 1 1 1 1 1 1 1 1 1 1
Columns 13 through 16
-649 362 -452 -537
-648 363 -451 -536
1 1 1 1
Hadamard Matrix H = hadamard(2*K) = ...
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1
1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1
1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1
1 1 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1
1 -1 1 -1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 1
1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1
1 -1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 1 1 -1
1 1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1
1 -1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1
1 1 -1 -1 1 1 -1 -1 -1 -1 1 1 -1 -1 1 1
1 -1 -1 1 1 -1 -1 1 -1 1 1 -1 -1 1 1 -1
1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1
1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1
1 1 -1 -1 -1 -1 1 1 -1 -1 1 1 1 1 -1 -1
1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 1 -1 -1 1
M = H/2^K ; Norm(M - M') = 0
Norm(M*M - eye(n)) = 0
new x = M*x ; [xlo, xhi, xhi-xlo]' = ...
Columns 1 through 7
-976.8 -943.3 665.3 745.8 204.3 -653.3 114.3
-972.8 -939.3 669.3 749.8 208.3 -649.3 118.3
4.0 4.0 4.0 4.0 4.0 4.0 4.0
Columns .8through 14
-134.3 485.8 190.3 28.8 1136.3 392.8 -445.8
-130.3 489.8 194.3 32.8 1140.3 396.8 -441.8
4.0 4.0 4.0 4.0 4.0 4.0 4.0
Columns 15 through 16
179.8 -839.8
183.8 -835.8
4.0 4.0
new x = M*x ; [xlo, xhi, xhi-xlo]' =
Columns 1 through 7
37.5 501.5 -918.5 461.5 620.5 -980.5 -1007.5
53.5 517.5 -902.5 477.5 636.5 -964.5 -991.5
16.0 16.0 16.0 16.0 16.0 16.0 16.0
Columns 8 through 14
247.5 -534.5 -21.5 -977.5 -149.5 -656.5 354.5
263.5 -518.5 -5.5 -961.5 -133.5 -640.5 370.5
16.0 16.0 16.0 16.0 16.0 16.0 16.0
Columns 15 through 16
-459.5 -544.5
-443.5 -528.5
16.0 16.0
new x = M*x ; [xlo, xhi, xhi-xlo]' =
Columns 1 through 7
-1006.8 -973.2 635.3 715.8 174.3 -683.3 84.3
-942.8 -909.3 699.3 779.8 238.3 -619.3 148.3
64.0 64.0 64.0 64.0 64.0 64.0 64.0
Columns 8 through 14
-164.3 455.8 160.3 -1.3 1106.3 362.8 -475.8
-100.3 519.8 224.3 62.8 1170.3 426.8 -411.8
64.0 64.0 64.0 64.0 64.0 64.0 64.0
Columns 15 through 16
149.8 -869.8
213.8 -805.8
64.0 64.0
new x = M*x ; [xlo, xhi, xhi-xlo]' =
Columns 1 through 7
-82.5 381.5 -1038.5 341.5 500.5 -1100.5 -1127.5
173.5 637.5 -782.5 597.5 756.5 -844.5 -871.5
256.0 256.0 256.0 256.0 256.0 256.0 256.0
Columns 8 through 14
127.5 -654.5 -141.5 -1097.5 -269.5 -776.5 234.5
383.5 -398.5 114.5 -841.5 -13.5 -520.5 490.5
256.0 256.0 256.0 256.0 256.0 256.0 256.0
Columns 15 through 16
-579.5 -664.5
-323.5 -408.5
256.0 256.0
È quit
41368 flop(s).