3 Signicant discrepancies (between the computed and the true result) are very rare, too rare to worry about all the time, yet not rare enough to ignore. by W.M.Kahan supercomputer 1 supercomputer supercomputer ( ) () supercomputer \GFLOPS 2 " \GBYTE 3 " (Fortran, C, C++, Pascal, etc.) ( ) () 3 1 JIS F 2 FLOPS G( ) 5GFLPOS 50 3 1GBYTE=1024MBYTE=1024 3 BYTE 1
1 (1) (2) ( ' ) \ " 4 5 6 (compute) \" (+,0, 2, 4 ) (computer) \" computer ( ) ( ) AIR TEX Graphics make 4 5 6 2
supercomputer supercomputer computer computer =computer 2 Computation Computation on the computer Algebraic Computation - computing without errors - ring of integers or field of rational numbers - algebraic number fields Numerical Computation - approximation of the exact result - floating-point arithmetic Verified Inclusions - approximation with a guaranteed error bound - interval arithmetic - fixed point theorem Mathematica CALC Fortran C++ Pascal ACRITH Pascal/XSC PROFIL 3
3 Computation Algebraic Computation algebraic computation Error-Free computation Algebraic computation 3 10 7 10 1 3 0:3333333 1 3 + 1 9 = 4 9 exact, sin 1 11, log 3 13 8 9689 11213 + 9941 19937 = 304638026 223553581 [11] Veri- ed Computation a; b 0 Z 1 0 (e 0 a2 x 2 0 e 0 b2 x 2 )dx =(b 0 a) p Algebraic Computation Numerical Computation Numerical Computation exact 7 2 ( 0 1 ) 8 4
exact 7 15 Numerical Computation Veried Inclusions Veried Inclusions (Veried computation) Algebraic Computation Numerical Computation (interval arithmetic) Numerical Computation p 2 p 2=1:4142 1112[1:414; 1:415] fx 2 IR; 1:414 x 1:415g Veried computation exact 1:0 x=[0.9999999999999993, 1.000000000000001] ( ) Veried Inclusions [3] Veried Inclusions [11] 5
3 = a 0 + c 1 10 + c 2 10 2 + 111+ c n 10 n + 111 a 0, c n 0 0 a 0 < 1, 0 c n < 10 10 [1] 10 1 t t = a 0 + c 1 t + c 2 t 2 + 111; (0 c n <t) (1) m (1) t m 0 = 6( x 1 t + x 2 t 2 + x 3 t 3 + 111)tm (2) x i 0 <x 1 <t,0x i <t(i =2; 3; 111) x i (oating point number) (2) n = 6( x 1 t + x 2 t + 111x n 2 t n )tm (3) m m L, m U 0m L m m U 4 n, t, m L, m U (t 0 1)t n01 (m L + m U +1) 2(t 0 1)t n01 (m L + m U +1)+1 n, t, m L, m U Fujitsu M1800/20 UXP/Fortran77 EX {z} 16 t {z} 14 n 064 {z} {z} m 63 m L m U 2215216 13 2 (64+63+1)+1 = 17293822569102704641 6 M1800/20
4 (round-o error) n n +1 (round-o error) 10 1X c i 10 i i=0 {z } = nx c i 10 i i=0 {z } n + 1X c i 10 i i=n+1 {z } ( ) 9 n Fortran 10... [2] 2 SSL II NUMPAC 9 3 10 7
2 ax 2 + bx + c =0; a 6= 0 x = 0b 6 p b 2 0 4ac 2a b p b 2 0 4ac [4] n n =5 a =1; b = 0(10 n +10 0n ); c =1 x 1 = 10 n x 2 = 10 0n x 1 = 100000 x 2 = 0:00001 SSL II RQDR Fujitsu M1800/20 Fortran77 EX VS. SSL II program niji real complex z(2) a/1.0/,b,c/1.0/,x1,x2 b = -1.0E5-1.0E-5 <--- n=5 x1 = ( -b + sqrt(b**2-4*a*c) )/(2.0*a) x2 = ( -b - sqrt(b**2-4*a*c) )/(2.0*a) print*,x1,x2 call rqdr(a,b,c,z,icon) print*,z end <--- <--- SSL II x1=100000.0 x2=0.0 <--- x1=100000.0 x2=0.0000099999997 <--- SSL II 8
x 2 0 SSL II (cf.[5]) 2 IF ([2]) 4 x1= 99999.999999999985 x2= 0.000013385357559 x1= 99999.99999999999999999999999999980 4 x2= 0.00001000000000000000000000000626847810 =) =) 0! time ( / / ) 9
5 vs. computer Fortran x 4 0 4y 2 0 4y 4 for x =665857:0 and y = 470832:0 z = x 4 0 4y 2 0 4y 4 10 5 0-20000 -40000-10 -5 0 5 10-10 -5 0 5 10 0-5 -10-10 -5 0 5 10 z = x 4 0 4y 2 0 4y 4 Fortran Fortran program ex1 real real*8 xr/665857.0e0/,yr/470832.0e0/ xd/665857.0d0/,yd/470832.0d0/ print*,'single precision : ',xr**4-4*yr**2-4*yr**4 print*,'double precision : ',xd**4-4*yd**2-4*yd**4 end 5... single precision : 7.2057594e+16 double precision : -16777216.000000000 10
program ex2 real xr/665857.0e0/,yr/470832.0e0/ real*8 xd/665857.0d0/,yd/470832.0d0/ print*,'single precision : ',xr**4-4*yr**4-4*yr**2 print*,'double precision : ',xd**4-4*yd**4-4*yd**2 end... single precision : 7.2056709e+16 double precision : -4891648.0000000000 x 4 : 196572970000000000000000 = 1:9657297 2 10 23 : 196573006004558190000000 = 1:9657300600455819 2 10 23 : 196573006004558194713601 CALC exact (algebraic computation) exact Fortran CALC c CALC {C-style arbitrary precision calculator. calc version: 1.25.0 - k1.0 Copyright (c) 1992 David I. Bell Modied 1993 by Masahide Kashiwagi CALC UXP, MSP 11
3 2 0 1 01> 665857**4-4*470832**4-4*470832**2 1 3 <--- 02> 665857**4 0 2 1 <--- 3 2 0 1 3 0 196573006004558194713601 03> 4*470832**4 196573006003671463624704 04> 4*470832**2 2 1 886731088896 05> 196573006004558194713601-196573006003671463624704 -886731088896 1 <--- 1 3 2 0 1 ( ) ( ) 1 x 4 0 4y 2 0 4y 4 for x =665857 and y = 470832 x 4, y 2, y 4 x 4 x 2 2 x 2 12
x 4 0 4y 2 0 4y 4 1 x 4 0 4y 2 0 4y 4 0 1=(x 2 +2y 2 +1)(x 2 0 2y 2 0 1) x 2, y 2 x 2 +2y 2 +1> 0 x 2 0 2y 2 =1 2 (x 2, y 2 ) 3 0 2 1 x 4 196573006004558194713601 196573006004558194713601 13
10 023 10 08 10 16 10 022 10 07 10 20 10 021 10 06 10 24 10 020 10 05 10 28 10 019 10 04 10 32 10 018 10 03 10 36 10 017 10 02 10 40 10 016 10 01 10 44 10 015 1 10 48 10 014 10 1 10 56 10 013 10 2 10 64 10 012 10 3 10 72 10 011 10 4 10 80 10 010 10 8 10 88 10 09 10 12 4 program ex3 real real*8 xr/665857.0e0/,yr/470832.0e0/ xd/665857.0d0/,yd/470832.0d0/ real*16 xq/665857.0q0/,yq/470832.0q0/ print*,'single precision print*,'double precision print*,'quadruple precision end single precision : 7.2057594e+16 double precision : -16777216.000000000 : ',xr**4-4*yr**2-4*yr**4 : ',xd**4-4*yd**2-4*yd**4 : ',xq**4-4*yq**2-4*yq**4 quadruple precision : 1.000000000000000000000000000000000 4 4 14
6 6 333:75b 6 + a 2 (11a 2 b 2 0 b 6 0 121b 4 0 2) + 5:5b 8 + a 2b (4) for a = 77617:0 and b = 33096:0 b 8 a=(2b) z =333:75y 6 + x 2 (11x 2 y 2 0 y 6 0 121y 4 0 2) + 5:5y 8 + x=(2y) x>0, y>0 \" 3000 200 1 100 0.8 0 0.6 0.2 0.4 0.4 0.6 0.8 0.2 1 z =333:75y 6 + x 2 (11x 2 y 2 0 y 6 0 121y 4 0 2) + 5:5y 8 + x=(2y) Fortran 1994 10 Fortran Fujitsu M1800/20 (4) f 1 (a; b) = 333:75b 6 +11a 4 b 2 0 b 6 a 2 0 121b 4 a 2 0 2a 2 +5:5b 8 + a=(2b) f 2 (a; b) = (333:75 + 5:5b 2 )b 6 + a 2 (11a 2 b 2 )+a 2 (0b 6 0 121b 4 0 2) + a=(2b) f 3 (a; b) = 333:75b 6 + a 2 (11a 2 b 2 0 b 6 0 121b 4 0 2) + 5:5b 8 + a=(2b) 4 f 1 = f 2 = f 3 15
f 1, f 2, f 3 C C C & & & program ex3 real real*8 a1/77617.0e0/,b1/33096.0e0/,f1 a2/77617.0d0/,b2/33096.0d0/,f2 real*16 a3/77617.0q0/,b3/33096.0q0/,f3 external f1,f2,f3 print*,'single precision print*,'double precision print*,'quadruple precision function f1(a,b) real a,b : ',f1(a1,b1) : ',f2(a2,b2) : ',f3(a3,b3) f1=333.75e0*b**6+11.0e0*a**4*b**2 -b**6*a**2-121.0e0*b**4*a**2-2.0e0*a**2 +5.5E0*b**8 + a/(2.0e0*b) return end real*8 function f2(a,b) real*8 a,b f2=(333.75d0 + 5.5E0*b**2)*b**6 + a**2*(11.0d0*a**2*b**2) + a**2*( -b**6-121.0d0*b**4-2.0d0) + a/(2.0d0*b) return end real*16 function f3(a,b) real*16 a,b f3=333.75q0*b**6+a**2*(11.0q0*a**2*b**2-b**6-121.0q0*b**4-2.0q0) +5.5Q0*b**8 + a/(2.0q0*b) return end () single precision : 1.1726036 double precision : 1.1726039400531785 quadruple precision : 1.172603940053178631858834904520183 7 1:172603... 16
Mathematica algebraic computation exact 11 ( ) CALC Mathematica 12 Mathematica Mathematica CALC In[1]:= a=77617 3 2 0 1 Out[1]= 77617 In[2]:= b=33096 3 2 0 1 Out[2]= 33096 3 2 0 1 In[3]:= 33375/100 b^6 + a^2(11 a^2 b^2 - b^6-121 b^4-2) + 55/10 b^8 + a/(2 b) 3 2 0 1 54767 Out[3]= -(-----) 66192 In[4]:= N[%,16] <--- exact value 3 2 0 1 <--- Out[4]= -0.827396059946821 <--- 6 Mathematica 00:827396 11 333:75 = 33375=100, 5:5 =55=10 12 Stephen Wolfram Fortran, C, TEX PostScript 17
7 1. LU 2. 3. band matrix band matrix full matrix LU ( band matrix ) Fortran LU [7] n 2 n A =(a ij ) 2 a ij =( n +1 ) ij 1 2 sin( n +1 ) A orthogonal symmetric matrix A 01 = A 1. Fortran77 EX + SSL II SSL II [5] DALU DLUIV DALU LU DLUIV LU 2. Fortran77 EX + NUMPAC NUMPAC[6] MINVD LU 3. gcc + utility Sun Sparc 10 gcc(c++) PROFIL [8] utility 18
n =10 n =150 A =(a ij ) numerical computation B =(b ij )( A = B ) error 13 error max 1i;jn ja ij 0 b ij j n =10 n =150 Fortran + SSL II 3:2058 2 10 015 5:4484 2 10 014 Fortran + NUMPAC 3:1503 2 10 015 1:2779 2 10 014 gcc + utility 7:2165 2 10 016 9:4386 2 10 015 gcc + utility Fortran [7] ( ) n =2; 8; 10 0 B@ 1 0 10 0n 01 1 01 1 01 1 01 1 01 0 0 1 0 0 01 10 CA B@ x 1 x 2 x 3 x 4 1 0 CA = B@ (5) n Algebraic Computation Mathematica n doronpa% math 3 0 Mathematica 2.0 for SPARC 03 02 01 03 21 <--- Mathematica Copyright 1988-91 Wolfram Research, Inc. -- Terminal graphics initialized -- In[1]:= m=ff1-10^(-n),-1,1,-1g,f1,-1,1,-1g,f1,-1,0,0g,f1,0,0,-1gg 1 CA 3 2 0 1 (5) 13 error 19
Out[1]= ff1-10 -n, -1, 1, -1g, f1, -1, 1, -1g, f1, -1, 0, 0g, f1, 0, 0, -1gg In[2]:= LinearSolve[m,f-3,-2,-1,-3g] 3 0 21 <--- n n n n n Out[2]= f10, 1 + 10, 4 + 3 10 + 2 (-1-10 ), 3 + 10 g In[3]:= Factor[%] 3 0 21 <--- n n n n Out[3]= f10, 1 + 10, 2 + 10, 3 + 10 g In[4]:= Quit 3 0 21 <--- Mathematica (5) Numerical Computation 0 B@ x 1 x 2 x 3 x 4 1 CA = 0B @ 10n 10 n +1 10 n +2 10 n +3 numerical computation Fortran n SSL II DLAX n=2 100.00000000000005 100 101.00000000000005 101 102.00000000000005 102 103.00000000000005 103 n =2 n n=8 100000000.05263561 100000000 100000001.05263561 100000001 100000002.05263561 100000002 100000003.05263561 100000003 n 20 1 CA
n=10 10000000560.375102 10000000000 10000000561.375102 10000000001 10000000562.375102 10000000002 10000000563.375102 10000000003 Veried Computation Hamburg-Harburg Olaf Knuppel C++ PROFIL ([8]) PROFIL solver S. M. Rump (cf.[9]) Theorem Let A2IPIR n2n, B2IPIR n be given and let ~x 2 IR n, R 2 IR n2n, ; 6= X 2 IPIR n, X being compact. Dene Z = R 1 (B 0A~x) and C = I 0 R 1A; L(X) =Z + C1X; all operations being power set operations. If L(X) int(x) then R and every A 2 IR n2n, A 2Ais nonsingular and for every b 2Bthe unique solution ^x = A 01 b satises ^x =~x + L(X): Ax = b ~x ^x A L(X) int(x) ~x + L(X) n PROFIL veried computation n=2 [ 99.99999999999959, 100.000000000001] [100.9999999999996, 101.000000000001] [101.9999999999996, 102.000000000001] [102.9999999999996, 103.000000000001] 21
n=8 [ 99999999.11763082, 100000002.0978633] [100000000.1176306, 100000003.0978635] [100000001.1176306, 100000004.0978635] [100000002.1176308, 100000005.0978633] n=10 [9999991201.131073, 10000029348.54311] [9999991201.971003, 10000029349.70318] [9999991202.971003, 10000029350.70318] [9999991204.131073, 10000029351.54311] n =10 veried computation algebraic computation Numerical Computation [10] n 2 n A A = 0 B@ 1 01 01 1 1 01 0........ 1 1 01 0 1 0 0 1 2 (6) NUMPAC HEQRVD n =17 Eigenvalue 1.1322859860 1.1229217217 6 0.0484689972i 1.0963368694 6 0.0897980776i 1.0566992036 6 0.1181569137i 22 1 CA (6)
1.0099364473 6 0.1299194454i 0.9626498727 6 0.1240285995i 0.9210696786 6 0.1018650121i 0.8902963972 6 0.0667226311i 0.8739468165 6 0.0232084483i z =1 exact A n 1 n =17 A 1 17 Mathematica (6) In[1]:= ReadList[``matrix.data'', Number, RecordLists -> True] Out[1]= ff1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1g, : ( ) In[2]:= Eigenvalues[%] : 3 2 0 1 <-- Out[2]= f1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1g Rump [10] W. M. Kahan 3 0 21 8 \" + \ " \ " \ " 23
(prime number) 1 2 3 5 7 11 13 17 19 23 111 ( ) 697 17 2 41 19937 Mathematica 10 ( ) etc. Colin Wilson & Damon Wilson \The Encyclopedia of Unsolved Mysteries" Colin Wilson [1] :, (1939). [2] :, 7, (1991). 24
[3] :,, Vol.31, No.9, 1177{1190 (1990). [4] Forsythe, G.E., Malcolm, M. A. and Moler, C. B. : Computer methods for mathematical computations, Prentice-Hall, Inc. (1977). [5] SSL II 99SP{0050, (1980). [6] ( : NUMPAC ), (1989). [7] Gregory, R. T. and Karney, D. L.: A collection of matrices for testing computational algorithms, John Wiley & Sons, New York (1969). [8] Knuppel, O. : PROFIL Programmer's runtime optimized fast interval library, Berichte des Forschungsschwerpunktes Informations- und Kommunikationstechnik, TUHH (1993). [9] Rump, S. M. : Solving algebraic problems with high accuracy, in Kulish, U. W., and Miranker, W. L., editors, A New Approach to Scientic Computation, Academic Press, New York, 51{120 (1983). [10] Rump, S. M. : Verication methods for dense and sparse system of equations, Berichte des Forschungsschwerpunktes Informations- und Kommunikationstechnik, TUHH (1993). [11], :, 831,, 53{72 (1993). 25