Chapter 24 - Chapter 24 Chapter Numerics Bjarne Stroustrup

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: Chapter 24 Chapter Numerics Bjarne Stroustrup www.stroustrup.com/Programming Abstract Abstract This lecture is an overview of the fundamental This tools and techniques for numeric computation. For some, numerics are everything. For many, numerics are occasionally essential. Here, we present the basic problems of size, precision, truncation, and error handling in standard mathematical functions. We present multidimensional matrices and the standard library complex numbers. library Stroustrup/Programming 2 Overview Overview Precision, overflow, sizes, errors Matrices Random numbers Complex numbers Please note: Again, we are exploring how to express the fundamental notions of an Again, application domain in code application Numerical computation, especially linear algebra We are showing the basic tools only. The computational techniques for We using these tools well, say in engineering calculations, you’ll have to learn elsewhere learn That’s “domain knowledge” beyond the scope of this presentation Stroustrup/Programming 3 Precision, etc. Precision, When we use the built-in types and usual computational When techniques, numbers are stored in fixed amounts of memory techniques, Floating-point numbers are (only) approximations of real numbers float x = 1.0/333; float sum = 0; for (int i=0; i<333; ++i) sum+=x; cout << sum << "\n"; 0.999999 Integer types represent (relatively) small integers only short y = 40000; int i = 1000000; cout << y << " " << i*i << "\n"; -25536 -727379968 There just aren’t enough bits to exactly represent every number There we need in a way that’s amenable to efficient computation we Stroustrup/Programming 4 Sizes Sizes char short int float long int* double The exact sizes of C++ built-in types depends on the hardware and The the compiler the These sizes are for Windows using Microsoft , GCC on Linux, and others sizeof(x) gives you the size of x in bytes By definition, sizeof(char)==1 By sizeof(char)==1 Unless you have a very good reason for something else, stick to Unless char, int, and double char int and double Stroustrup/Programming 5 Sizes, overflow, truncation Sizes, // when you calculate, you must be aware of possible overflow and truncation // when // Beware: C++ will not catch such problems for you // Beware: void f(char c, short s, int i, long lg, float fps, double fpd) { c = i; // yes: chars really are very small integers // yes: s = i; i = i+1; // what if i was the largest int? // what lg = i*i; // beware: a long may not be any larger than an int // beware: // and anyway, i*i is an int – possibly it already overflowed // and fps = fpd; i = fpd; // truncates: e.g. 5.7 -> 5 fps = i; // you can lose precision (for very large int values) // you char ch = 0; for (int i = 0; i<500; ++i) { cout << int(ch) << "\t"; ch++; } // try this for try } Stroustrup/Programming 6 If in doubt, you can check If The simplest way to test Assign, then compare void f(int i) { char c = i; // may throw away information you don’t want to loose // may if (c != i) { // oops! We lost information, figure out what to do // oops! } // … // } Stroustrup/Programming 7 Math function errors Math If a standard mathematical function can’t do its If job, it sets errno from <cerrno>, for example errno void f(double negative, double very_large) // primitive (1974 vintage, pre-C++) but ISO standard error handling // primitive { errno = 0; // no current errors // no sqrt(negative); // not a good idea // not if (errno) { /* … */ } // errno!=0 means ‘something went wrong” /* // errno!=0 if (errno == EDOM) // domain error // domain cerr << "sqrt() not defined for negative argument"; cerr pow(very_large,2); // not a good idea // not if (errno==ERANGE) // range error // range cerr << "pow(" << very_large << ",2) too large for a double"; cerr } Stroustrup/Programming 8 Matrices Matrices The standard vector and the built-in array are The vector one dimensional one What if we need 2 dimensions? (e.g. a matrix) N dimensions? A vector (e.g. Matrix<int> v(4) ), also called a 1 dimensional Matrix, or even a 1*N matrix A 3*4 matrix (e.g. Matrix<int> m(3,4) ), also called a 2 dimensional Matrix Stroustrup/Programming 9 Matrices Matrices Matrix<int> m(3,4); A column A row A 3*4 matrix, also called a 2 dimensional Matrix 3 rows 4 columns Stroustrup/Programming 10 C-style multidimensional Arrays C-style A built-in facility int ai[4]; double ad[3][4]; double char ac[3][4][5]; ai[1] =7; ad[2][3] = 7.2; ac[2][3][4] = 'c'; // 1 dimensional array // dimensional // 2 dimensional array // dimensional // 3 dimensional array // dimensional Basically, Arrays of Arrays Stroustrup/Programming 11 C-style multidimensional Arrays C-style Problems C-style multidimensional Arrays are Arrays of Arrays Fixed sizes (i.e. fixed at compile time) Can’t be passed cleanly Not even assignment (copy) A major source of bugs major As usual, an Array doesn’t know its own size No Array operations An Array turns into a pointer at the slightest provocation No range checking If you want to determine a size at run-time, you’ll have to use free store And, for most people, they are a serious pain to write Look them up only if you are forced to use them Look only E.g. TC++PL, Appendix C, pp 836-840 Stroustrup/Programming 12 C-style multidimensional Arrays C-style You can’t pass multidimensional Arrays cleanly void f1(int a[3][5]); Can’t read vector with size from input and then call f1 // useful for [3][5] matrices only // useful (unless the size happens to be 3*5) Can’t write a recursive/adaptive function void f2(int [ ][5], int dim1); // 1st dimension can be a variable // void f3(int[ ][ ], int dim1, int dim2); void // error (and wouldn’t work anyway) error void f4(int* m, int dim1, int dim2) // odd, but works // odd, { for (int i=0; i<dim1; ++i) for (int j = 0; j<dim2; ++j) m[i*dim2+j] = 0; for } Stroustrup/Programming 13 A Matrix library Matrix // on the course site Matrix.h // on #include “Matrix.h” void f(int n1, int n2, int n3) { Matrix<double> ad1(n1); // default: one dimension // default: Matrix<int,1> ai1(n1); Matrix<double,2> ad2(n1,n2); // 2 dimensional // dimensional Matrix<double,3> ad3(n1,n2,n3); // 3 dimensional // dimensional ad1(7) = 0; // subscript using ( ) – Fortran style // subscript ad1[7] = 8; // [ ] also works – C style // also style ad2(3,4) = 7.5; // true multidimensional subscripting // true ad3(3,4,5) = 9.2; } Stroustrup/Programming 14 A Matrix library Matrix “like your math/engineering textbook talks about Matrixs” Compile-time and run-time checked Matrices of any dimension You can pass them around Usual Matrix operations 1, 2, and 3 dimensions actually work (you can add more if/as needed) Matrices are proper variables/objects Or about vectors, matrices, tensors Subscripting: ( ) Slicing: [ ] Assignment: = Scaling operations (+=, -=, *=, %=, etc.) Fused vector operations (e.g., res[i] = a[i]*c+b[2]) Fused res[i] Dot product (res = sum of a[i]*b[i]) Dot a[i]*b[i] Performs equivalently to hand-written low-level code You can extend it yourself as needed (“no magic”) Stroustrup/Programming 15 A Matrix library Matrix // compile-time and run-time error checking // compile-time void f(int n1, int n2, int n3) { Matrix<double> ad1(5); // default: one dimension // default: Matrix<int> ai(5); Matrix<double> ad11(7); Matrix<double,2> ad2(n1); // error: length of 2nd dimension missing // error: Matrix<double,3> ad3(n1,n2,n3); Matrix<double,3> ad33(n1,n2,n3); ad1(7) = 0; // Matrix_error exception; 7 is out of range // Matrix_error ad1 = ai; // error: different element types // error: ad1 = ad11; // Matrix_error exception; different dimensions ad2(3) = 7.5; // error: wrong number of subscripts // error: ad3 = ad33; // ok: same element type, same dimensions, same lengths ok: } Stroustrup/Programming 16 A Matrix library Matrix • As we consider the matrix (row, column): a[1][2] Matrix<int> a(3,4); a(1,2) a[0]: 01 02 03 a[1]: 10 11 12 13 a[2]: 00 20 21 22 23 As the elements are laid out in memory: 00 01 02 03 10 11 12 13 Stroustrup/Programming 20 21 22 23 17 A Matrix library Matrix void init(Matrix<int,2>& a) // initialize each element to a characteristic value // initialize { for (int i=0; i<a.dim1(); ++i) for (int j = 0; j<a.dim2(); ++j) a(i,j) = 10*i+j; } void print(const Matrix<int,2>& a) // print the elements of a row by row // print { for (int i=0; i<a.dim1(); ++i) { for (int j = 0; j<a.dim2(); ++j) cout << a(i,j) << '\t'; cout << '\n'; } } Stroustrup/Programming 18 2D and 3D Matrices 2D // 2D space (e.g. a game board): // 2D enum Piece { none, pawn, knight, queen, king, bishop, rook }; Matrix<Piece,2> board(8,8); // a chessboard // chessboard Piece init_pos = { rook, knight, bishop, queen, king, bishop, knight, rook }; // 3D space (e.g. a physics simulation using a Cartesian grid): // 3D int grid_nx; // grid resolution; set at startup // grid int grid_ny; int grid_nz; Matrix<double,3> cube(grid_nx, grid_ny, grid_nz); Stroustrup/Programming 19 1D Matrix 1D Matrix<int> a(10); // means Matrix<int,1> a(10); // means a.size(); a.dim1(); int* p = a.data(); a(i); a[i]; Matrix<int> a2 = a; a = a2; a *= 7; a.apply(f); a.apply(f,7); a.apply(f,7); b =apply(f,a); b =apply(f,a,7); =apply(f,a,7); // number of elements // number // number of elements number // extract data as a pointer to a C-style array extract // ith element (Fortran style), but range checked // // ith element (C-style), but range checked // // copy initialization copy // copy asssignment // copy // scaling a(i)*=7 for each i (also +=, -=, /=, etc.) // scaling // a(i)=f(a(i)) for each element a(i) a(i)=f(a(i)) // a(i)=f(a(i),7) for each element a(i) // a(i)=f(a(i),7) // make a new Matrix with b(i)==f(a(i)) // make a new Matrix with b(i)==f(a(i),7) Matrix<int> a3 = scale_and_add(a,8,a2); // fused multiply and add Matrix<int> // fused int r = dot_product(a3,a); // dot product // dot Stroustrup/Programming 20 2D Matrix (very like 1D) (very Matrix<int,2> a(10,20); a.size(); a.dim1(); a.dim2(); int* p = a.data(); a(i,j); a[i]; a[i][j]; Matrix<int> a2 = a; a = a2; a *= 7; a.apply(f); a.apply(f,7); a.apply(f,7); b=apply(f,a); b=apply(f,a,7); b=apply(f,a,7); a.swap_rows(7,9); a.swap_rows(7,9); // number of elements // number // number of elements in a row number // number of elements in a column number // extract data as a pointer to a C-style array extract // (i,j)th element (Fortran style), but range checked // (i,j) // ith row (C-style), but range checked // // (i,j)th element C-style (i,j) // copy initialization copy // copy asssignment // copy // scaling (and +=, -=, /=, etc.) // scaling // a(i,j)=f(a(i,j)) for each element a(i,j) a(i,j)=f(a(i,j)) // a(i,j)=f(a(i,j),7) for each element a(i,j) // a(i,j)=f(a(i,j),7) // make a new Matrix with b(i,j)==f(a(i,j)) // // make a new Matrix with b(i,j)==f(a(i,j),7) // // swap rows a[7] a[9] // Stroustrup/Programming 21 3D Matrix (very like 1D and 2D) (very Matrix<int,3> a(10,20,30); a.size(); a.dim1(); a.dim2(); a.dim3(); int* p = a.data(); a(i,j,k); a[i]; a[i][j][k]; Matrix<int> a2 = a; a = a2; a *= 7; a.apply(f); a.apply(f,7); a.apply(f,7); b=apply(f,a); b=apply(f,a,7); b=apply(f,a,7); a.swap_rows(7,9); a.swap_rows(7,9); // number of elements // number // number of elements in dimension 1 number // number of elements in dimension 2 number // number of elements in dimension 3 number // extract data as a pointer to a C-style Matrix extract // (i,j,k)th element (Fortran style), but range checked // (i,j,k) // ith row (C-style), but range checked // // (i,j,k)th element C-style (i,j,k) // copy initialization copy // copy asssignment // copy // scaling (and +=, -=, /=, etc.) // scaling // a(i,j,k)=f(a(i,j)) for each element a(i,j,k) a(i,j,k)=f(a(i,j)) // a(i,j,k)=f(a(i,j),7) for each element a(i,j,k) // a(i,j,k)=f(a(i,j),7) // make a new Matrix with b(i,j,k)==f(a(i,j,k)) // make a new Matrix with b(i,j,k)==f(a(i,j,k),7) // swap rows a[7] a[9] // Stroustrup/Programming 22 Using Matrix Using See book Matrix I/O §24.5.4; it’s what you think it is Solving linear equations example §24.6; it’s just like your algebra textbook says it is Stroustrup/Programming 23 Random numbers Random A “random number” is a number from a sequence that matches “random a distribution, but where it is hard to predict the next number in the sequence the Uniformly distributed integers <cstdlib> <cstdlib> Uniformly // a value in [0:RAND_MAX] value // the largest possible value for rand() // the // seed the random number generator seed Use int rand() RAND_MAX void srand(unsigned); int rand_int(int max) { return()%max; } int rand_int(int min, int max) {return min+rand_int(maxmin); } If that’s not good enough (not “random enough”) or you If need a nonuniform distribution, use a professional library need E.g. boost::random (also C++0x) E.g. (also Stroustrup/Programming 24 Complex Complex Standard library complex types from <complex> Standard <complex> template<class T> class complex { T re, im; // a complex is a pair of scalar values, a coordinate pair re, complex public: public: complex(const T& r, const T& i) :re(r), im(i) { } complex(const complex(const T& r) :re(r),im(T()) { } complex() :re(T()), im(T()) { } T real() { return re; } T imag() { return im; } // operators: = += -= *= /= // operators: }; // operators: + - / * == != // operators: // whatever standard mathematical functions that apply to complex: // whatever // pow(), abs(), sqrt(), cos(), log(), etc. and also norm() (square of abs()) // pow(), Stroustrup/Programming 25 Complex Complex // use complex<T> exactly like a built-in type, such as double // use // just remember that not all operations are defined for a complex (e.g. <) // just typedef complex<double> dcmplx; // sometimes complex<double> gets verbose typedef sometimes void f( dcmplx z, vector< complex<float> >& vc) {; dcmplx z2 = pow(z,2); dcmplx z3 = z2*9+vc[3]; dcmplx sum = accumulate(vc.begin(), vc.end(), dcmplx); } Stroustrup/Programming 26 Numeric limits Numeric Each C++ implementation specifies properties of the built-in types built-in used to check against limits, set sentinels, etc. From <limits> From <limits> for each type // smallest value // largest value For floating-point types min() max() … Lots (look it up if you ever need it) E.g. max_exponent10() E.g. max_exponent10() From <limits.h> and <float.h> <float.h> INT_MAX DBL_MIN // largest int value // largest value // smallest double value // smallest value Stroustrup/Programming 27 Numeric limits Numeric They are important to low-level tool builders If you think need them, you are probably too close to hardware, If but there are a few other uses. For example, but void f(const vector<int>& vc) { // pedestrian (and has a bug): // pedestrian int smallest1 = v[0]; for (int i = 1; i<vc.size(); ++i) if (v[i]<smallest1) smallest1 = v[i]; // better: // better: int smallest2 = numeric_limits<int>::max(); for (int i = 0; i<vc.size(); ++i) if (v[i]<smallest2) smallest2 = v[i]; // or use standard library: // or vector<int>::iterator p = min_element(vc.begin(),vc.end()); // and check for p==vc.end() // and } Stroustrup/Programming 28 A great link great http://www-gap.dcs.st-and.ac.uk/~history/ A great link for anyone who likes math Famous mathematicians Biography, accomplishments, plus curio, for example, who is the only major mathematician to win an Olympic medal? Famous curves Famous problems Mathematical topics Or simply needs to use math Algebra Algebra Analysis Analysis Numbers and number theory Numbers Geometry and topology Geometry Mathematical physics Mathematical Mathematical astronomy The history of mathematics … Stroustrup/Programming 29 ...
View Full Document

This note was uploaded on 02/18/2012 for the course CSCE 121 taught by Professor Walter daugherity during the Fall '09 term at Texas A&M.

Ask a homework question - tutors are online