main.dvi

Size: px

Start display at page:

Download "main.dvi"

えつみえいさか
5 years ago
Views:

1 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 GBYTE=1024MBYTE= BYTE 1

2 1 (1) (2) ( ' ) \ " (compute) \" (+,0, 2, 4 ) (computer) \" computer ( ) ( ) AIR TEX Graphics make

supercomputer supercomputer computer computer =computer 2 Computation Computation on the computer Algebraic Computation - computing without errors - ring of integers or field of rational numbers -

3 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

4 3 Computation Algebraic Computation algebraic computation Error-Free computation Algebraic computation : = 4 9 exact, sin 1 11, log = [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

5 exact 7 15 Numerical Computation Veried Inclusions Veried Inclusions (Veried computation) Algebraic Computation Numerical Computation (interval arithmetic) Numerical Computation p 2 p 2=1: [1:414; 1:415] fx 2 IR; 1:414 x 1:415g Veried computation exact 1:0 x=[ , ] ( ) 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

6 3 = a 0 + c c c n 10 n a 0, c n 0 0 a 0 < 1, 0 c n < [1] 10 1 t t = a 0 + c 1 t + c 2 t ; (0 c n <t) (1) m (1) t m 0 = 6( x 1 t + x 2 t 2 + x 3 t )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 ( )+1 = M1800/20

7 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 [2] 2 SSL II NUMPAC

8 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 = 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= x2=0.0 <--- x1= x2= <--- SSL II 8

x 2 0 SSL II (cf.[5]) 2 IF ([2]) 4 x1= 99999.

9 x 2 0 SSL II (cf.[5]) 2 IF ([2]) 4 x1= x2= x1= x2= =) =) 0! time ( / / ) 9

10 5 vs. computer Fortran x 4 0 4y 2 0 4y 4 for x =665857:0 and y = :0 z = x 4 0 4y 2 0 4y z = x 4 0 4y 2 0 4y 4 Fortran Fortran program ex1 real real*8 xr/ e0/,yr/ e0/ xd/ d0/,yd/ d0/ 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 : e+16 double precision :

11 program ex2 real xr/ e0/,yr/ e0/ real*8 xd/ d0/,yd/ d0/ 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 : e+16 double precision : x 4 : = 1: : = 1: : CALC exact (algebraic computation) exact Fortran CALC c CALC {C-style arbitrary precision calculator. calc version: k1.0 Copyright (c) 1992 David I. Bell Modied 1993 by Masahide Kashiwagi CALC UXP, MSP 11

12 > **4-4*470832**4-4*470832**2 1 3 <--- 02> ** < > 4*470832** > 4*470832** > < ( ) ( ) 1 x 4 0 4y 2 0 4y 4 for x = and y = x 4, y 2, y 4 x 4 x 2 2 x 2 12

13 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 ) x

14 program ex3 real real*8 xr/ e0/,yr/ e0/ xd/ d0/,yd/ d0/ real*16 xq/ q0/,yq/ q0/ print*,'single precision print*,'double precision print*,'quadruple precision end single precision : e+16 double precision : : ',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 :

15 :75b 6 + a 2 (11a 2 b 2 0 b b 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 y 4 0 2) + 5:5y 8 + x=(2y) x>0, y>0 \" z =333:75y 6 + x 2 (11x 2 y 2 0 y y 4 0 2) + 5:5y 8 + x=(2y) Fortran Fortran Fujitsu M1800/20 (4) f 1 (a; b) = 333:75b 6 +11a 4 b 2 0 b 6 a b 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 b 4 0 2) + a=(2b) f 3 (a; b) = 333:75b 6 + a 2 (11a 2 b 2 0 b b 4 0 2) + 5:5b 8 + a=(2b) 4 f 1 = f 2 = f 3 15

16 f 1, f 2, f 3 C C C & & & program ex3 real real*8 a1/ e0/,b1/ e0/,f1 a2/ d0/,b2/ d0/,f2 real*16 a3/ q0/,b3/ q0/,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** e0*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.75d E0*b**2)*b**6 + a**2*(11.0d0*a**2*b**2) + a**2*( -b** d0*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** q0*b**4-2.0q0) +5.5Q0*b**8 + a/(2.0q0*b) return end () single precision : double precision : quadruple precision : :

17 Mathematica algebraic computation exact 11 ( ) CALC Mathematica 12 Mathematica Mathematica CALC In[1]:= a= Out[1]= In[2]:= b= Out[2]= 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) Out[3]= -(-----) In[4]:= N[%,16] <--- exact value <--- Out[4]= <--- 6 Mathematica 00: :75 = 33375=100, 5:5 =55=10 12 Stephen Wolfram Fortran, C, TEX PostScript 17

18 7 1. LU 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

19 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: : Fortran + NUMPAC 3: : gcc + utility 7: : gcc + utility Fortran [7] ( ) n =2; 8; 10 0 B@ n CA B@ x 1 x 2 x 3 x CA = B@ (5) n Algebraic Computation Mathematica n doronpa% math 3 0 Mathematica 2.0 for SPARC <--- Mathematica Copyright 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 (5) 13 error 19

20 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] <--- n n n n n Out[2]= f10, , (-1-10 ), g In[3]:= Factor[%] <--- n n n n Out[3]= f10, , , g In[4]:= Quit <--- Mathematica (5) Numerical Computation 0 B@ x 1 x 2 x 3 x 4 1 CA = 10n 10 n n n +3 numerical computation Fortran n SSL II DLAX n= n =2 n n= n 20 1 CA

21 n= 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 [ , ] [ , ] [ , ] [ , ] 21

22 n=8 [ , ] [ , ] [ , ] [ , ] n=10 [ , ] [ , ] [ , ] [ , ] n =10 veried computation algebraic computation Numerical Computation [10] n 2 n A A = 0 B@ (6) NUMPAC HEQRVD n =17 Eigenvalue i i i 22 1 CA (6)

23 i i i i i 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[%] : <-- Out[2]= f1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1g Rump [10] W. M. Kahan \" + \ " \ " \ " 23

24 (prime number) ( ) 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).

25 [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

II

II 2016 7 21 computer-assisted proof 1 / 64 1. 2. 3. Siegfried M. Rump : [1] I,, 14:3 (2004), pp. 214 223. [2] II,, 14:4 (2004), pp. 346 359. 2 / 64 Risch 18 3 / 64 M n = 2 n 1 (n = 1, 2,... ) 2 2 1 1