1 (bit ) ( ) PC WS CPU IEEE754 standard ( 24bit) ( 53bit)

Size: px

Start display at page:

Download "1 (bit ) ( ) PC WS CPU IEEE754 standard ( 24bit) ( 53bit)"

みそらいいはた
7 years ago
Views:

1 GNU MP BNCpack ( ) Linux Conference

2 1 (bit ) ( ) PC WS CPU IEEE754 standard ( 24bit) ( 53bit)

3 :9:; <>= R )TSVU!" # %$ & " # %$ & " # %' (*) " # %+-, " # " #./ " # " # 0 1 " # %' (*) " # 2, " # 1: IEEE754 IEEE754 3

4 CC "!$#&%('*)+,.- / :93; )+,.-</1032 = C E8@ EB@ 2: " #$%'& (*) +%'&!! 3: 4

5 PC IEEE754 5

6 open source software proprietary software GNU MP[6]( GMP ) GMP BNCpack[5] 6

7 2 (integer) (floating-point number) CPU (bit ) bit N 2 N+1 (overflow) 7

8 (base) (mantissa) (exponent) (round) 8

9 IEEE754 standard 2 IEEE754 IEEE (23) (52) bit 1bit bit p ε M ε M = 2 (p 1) (1) 0 ε M IEEE754 standard , 16 ε M 1 9

10 2 IEEE754 standard 2 10

11 2.1 ( ) 3 (ill-conditioned) n Ax = b A cond(a) = A A 1 x

12 Logistic x 0 = x n+1 := 4x n (1 x n ) {x n } 1 x n x n x 100 = IEEE [2] 12

13 2.2 A x A = , x = Conjugate-Gradient IEEE754 (float), (double) 128bit, 256bit, 512bit, 1024bit r k 2 = b Ax k 2 r

14 b-axi 2 1 IEEE754 float 1e-05 IEEE754 double 1e-10 1e-15 gmp 128bits gmp256bits 1e-20 gmp 1024bits gmp 512bits : CG IEEE754 CG 512bit IEEE754 IEEE754 14

15 10 15 Pentium IV IEEE GMP 4.1 (bits)

16 3 GNU MP GNU MP( GMP ) 1991 ANSI C (mpz t) (mpq t) Version 2 (mpf t) BNCpack GMP C 2048bit #include <stdio.h> #include "gmp.h" main() { int i; mpf_t x[102]; mpf_set_default_prec(2048); mpf_init_set_str(x[0], "0.7501", 10); for(i = 1; i <= 101; i++) mpf_init(x[i]); } for(i = 0; i <= 100; i++) { mpf_ui_sub(x[i+1], 1UL, x[i]); mpf_mul(x[i+1], x[i], x[i+1]); mpf_mul_ui(x[i+1], x[i+1], 4UL); if(i%10 == 0) gmp_printf("%5d %58.50Fe\n", i, x[i]); } 16

17 2 17

18 Version 4 C++ GMP Project MPFR [7] 4 C++ C #include <iostream.h> #include "gmpxx.h" main() { mpf_set_default_prec(2048); int i; mpf_class x[102]; x[0] = "0.7501"; } for(i = 0; i <= 100; i++) { x[i+1] = 4 * x[i] * (1 - x[i]); if(i%10 == 0) gmp_printf("%5d %58.50Fe\n", i, x[i].get_mpf_t()); } 4 MPFR GMP mpf t IEEE754 standard 18

19 GMP MPFR 1CPU PC GMP BNCpack GMP CPU Intel Pentium III 750MHz RAM 384MB OS Windows XP Professional C compiler GCC 2.95 GMP Ver (./configure; make; make install ) Software MuPAD Pro 2.0, Mathematica 4.1 MPFR st = cputime(); for (i=0;i<n;i++) mpf_add(z, x, y); printf("x+y took %1.2ems\n", (double)(cputime()-st)/n);

20 CPU GMP order GMP GMP GMP Knuth[9] 20

21 1: IEEE754 IEEE x + y x y x y x/y x : MuPAD Pro Mathematica GMP x + y x y x y x/y x x + y x y x y x/y x x + y x y x y x/y x

22 CPU Intel Pentium IV 2.2GHz RAM 512MB OS RedHat Linux 7.3 C compiler GCC 2.96 GMP & MPFR Ver.3.1.1, Ver.4.1 (./configure; make; make install ) 22

23 GMP 3.1.1, 4.1, MPFR GMP GMP 4.1 GMP 4.1& MPFR x + y x y x y x/y x x + y x y x y x/y x x + y x y x y x/y x

24 4 BNCpack BNCpack Basic Numerical Computation PACKage GMP Version 2 GMP GMP Fortran GMP Version 3 ANSI C C(not C++) 24

25 buggy 1. GMP (Ver.0.3 MPFR ) 2. (C++ Complex GSL[8] ) 3. ( ( ) ) 4. ( CG ) 5. ( ) 6. (Newton Newton, Regula-Falsi ) 7. (DKA ) 8. ( Romberg ) 9. ( Runge-Kutta ) 25

26 BNCpack GMP typdef MPFCmplx MPFPoly MPFArray, CMPFArray( ) MPFStack MPFVector MPFMatrix typedef struct{ unsigned long int prec; mpf_t re; mpf_t im; } mpfcmplx; typedef mpfcmplx *MPFCmplx; prec (bit ) BNCpack GMP mpf t BNCpack 26

27 4.1 #ifdef USE_GMP /* GMP */ /* */ MPFCmplx mpfca, mpfcb, mpfcc; /* */ /* Default: 128bit */ set_bnc_default_prec(128); /* Default */ mpfca = init_mpfcmplx(); mpfcb = init_mpfcmplx(); /* 256bit */ mpfcc = init2_mpfcmplx(256); /* mpfa := 1 + 2i */ /* mpfb := 3 + 4i */ /* mpfc := rand() + rand() i */ set_real_mpfcmplx_ui(mpfca, 1); set_image_mpfcmplx_ui(mpfca, 2); set_real_mpfcmplx_ui(mpfcb, 3); set_image_mpfcmplx_ui(mpfcb, 4); set_real_mpfcmplx_ui(mpfcc, rand()); set_image_mpfcmplx_ui(mpfcc, rand()); /* mpfcc := mpfca + mpfcb */ add_mpfcmplx(mpfcc, mpfca, mpfcb); /* */ printf("mpfca + mpfcb = "); print_mpfcmplx(mpfcc); /* */ free_mpfcmplx(mpfca); free_mpfcmplx(mpfcb); free_mpfcmplx(mpfcc); #endif 27

28 4.2 mpfmat a, mpfmat b #ifdef USE_GMP /* */ MPFMatrix mpfmat, mpfmat_a, mpfmat_b; mpf_t tmp, sqrt2; long int dim, row_dim, col_dim, i, j; /* Default: 128bit */ set_bnc_default_prec(128); /* */ mpfmat = init_mpfmatrix(row_dim, col_dim); mpfmat_a = init_mpfmatrix(row_dim, col_dim); mpfmat_b = init_mpfmatrix(row_dim, col_dim); mpf_init(tmp); mpf_init(sqrt2); mpf_sqrt_ui(sqrt2, 2UL); /* mpfmat[i][j] := sqrt(2) * (i+1) / (j+1) */ for(i = 0; i < row_dim; i++) { for(j = 0; j < col_dim; j++) { mpf_mul_ui(tmp, sqrt2, i+1); mpf_div_ui(tmp, tmp, j+1); set_mpfmatrix_ij(mpfmat, i, j, tmp); } } /* mpfmat_a := mpfmat */ subst_mpfmatrix(mpfmat_a, mpf_mat); /* mpfmat_b := mpfmat * mpfmat_a */ mul_mpfmatrix(mpfmat_b, mpfmat, mpfmat_a); /* */ printf("mpfmat_b:\n"); print_mpfmatrix(mpfmat_b); /* */ mpf_clear(tmp); mpf_clear(sqrt2); free_mpfmatrix(mpfmat); free_mpfmatrix(mpfmat_a); free_mpfmatrix(mpfmat_b); #endif 28

29 4.3 LU #ifdef USE_GMP /* Default: 256bit */ set_bnc_default_prec(256); /* */ mpf_init(reps); mpf_init(aeps); mpfa = init_mpfmatrix(dim, DIM); mpfb = init_mpfvector(dim); mpfx = init_mpfvector(dim); mpfans = init_mpfvector(dim); /* A b, x */ get_mpfproblem(mpfa, mpfb, mpfans); /* A */ print_mpfmatrix(mpfa); /* LU */ ret_mpf = MPFLUdecomp(mpfa); ret_mpf = SolveMPFLS(mpfx, mpfa, mpfb); /* */ for(i = 0; i < DIM; i++) { printf("%5ld ", i); mpf_out_str(stdout, 10, 0, \ get_mpfvector_i(mpfx, i)); printf(" "); mpf_out_str(stdout, 10, 0, \ get_mpfvector_i(mpfans, i)); printf("\n"); } /* */ mpf_clear(reps); mpf_clear(aeps); free_mpfmatrix(mpfa); free_mpfvector(mpfb); free_mpfvector(mpfx); free_mpfvector(mpfans); #endif 29

30 5 BNCpack BNCpack 1. ANSI C ANSI C GSL[8] C++ template BLAS LAPACK[10] 4. 30

31 5.1 BNCpack C GMP 2 C++ C++ STL 31

32 5.2 GSL(GNU Scientific Library)[8] BNCpack Bessel IEEE754 standard LAPACK(+BLAS) LAPACK 32

33 5.3 ( ) Mathematica Mathematica n := 50; digits := 100; a := Table[N[Sqrt[2]*i/j, digits], {i, 1, n}, {j, 1, n}]; b := Table[N[Sqrt[2]*i/j, digits], {i, 1, n}, {j, 1, n}]; time := Timing[a. b;]; Print[time]; Clear[a, b, time, n, digits]; 33

34 Mathematica BNCpack Mathematica BNCpack Mathematica BNCpack Mathematica BNCpack Mathematica 2 3 Mathematica Mathematica BNCpack Math./BNC

35 ( 3) 1 35

36 6 CPU Mathematica, MuPAD GMP GMP GMP BNCpack 36

37 BNCpack bug ANSI C 4. BNCpack free open 37

38 [1],,, [2], [3],,, [4],,, [5] BNCpack, [6] GNU MP, [7] The MPFR Library, [8] The GNU Scientific Library, [9] D.E.Knuth/, /,,1986. [10] LAPACK, [11] MuPAD, [12] Wolfram Research, 38

Krylov (b) x k+1 := x k + α k p k (c) r k+1 := r k α k Ap k ( := b Ax k+1 ) (d) β k := r k r k 2 2 (e) : r k 2 / r 0 2 < ε R (f) p k+1 :=

Krylov (b) x k+1 := x k + α k p k (c) r k+1 := r k α k Ap k ( := b Ax k+1 ) (d) β k := r k r k 2 2 (e) : r k 2 / r 0 2 < ε R (f) p k+1 := 127 10 Krylov Krylov (Conjugate-Gradient (CG ), Krylov ) MPIBNCpack 10.1 CG (Conjugate-Gradient CG ) A R n n a 11 a 12 a 1n a 21 a 22 a 2n A T = =... a n1 a n2 a nn n a 11 a 21 a n1 a 12 a 22 a n2 = A...