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