AA 2004/2005

Course of Scientific Programming in C++

Suggested lectures

• Bruce Eckel’s Free Electronic Books:

  • Here you can find some free books about C++

• C++ Annotations

  • Here you can find the C++ Annotations by Frank B. Brokken

Lessons

Lesson of 16/5/2005

1.cpp

File: 1.cpp

/*
   A multi line comment
   start with a slash + star
   end with a star + slash
 */

// single line comment

// include
#include <iostream>

/*
 A program consist in a procedure named "main"
 with o without arguments
 */
 
int // means integer is the
    // return type of the routine main
main () // main is the name "fixed" of the program
        // () contains eventual arguments
{ // open brace is the beginning of the body of the program

  // cout = character out is the object for the output
  // cout write to the standard output 
  std::cout << "hello world!\n" ;
  // << is an operator, that in the case of the object cout means
  // that the thing on the right is put in the standard output.
  // the character \n is the newline. There are others special character
  // \r = is the return character 
  // \t = is the tabulation character
  // \g = the bell sound
  
  // first possibility to avoid the use of std::
  using std::cout ;
  cout << "cout without std::\r"
       << "XYZ\n" ;
           
  // use of endl
  using namespace std ;
  cout << "ABC" << endl ; // endl = "\n"
  // but endl also empty the buffer
  // endl is equivalent to the compose command
  cout << "a string" << endl ;
  cout << "a string" << "\n" << flush ;
  
  return 0 ;
} // closing brace is the end of the program

2.cpp

File: 2.cpp

#include <iostream>

using namespace std ;

int
main () {
  int a, b ; // define the integer a and b
  
  /*
   
   Other numerical type
   
   char   = character (typically 8 bits) 
   int    = integer of the size of the word of the machine
            (typically 32 bits, 4 byte)
   float  = floting point number (normally 32, bits 4byte)
   double = double precision floating point number
            (typically 64 bits, 8 bytes). 
   
   Modificator
   
   usigned
   signed
   
   [ unsigned int  usigned char
     signed int    signed char ]

   short
   long

   short int = small integer (typically 16 bits )
   long int  = large integer (typically 64 bits )
   
   "short" is a synonimous of "short int"
   "long"  is a synonimous of "long int"
   
   long double = quadruple precision floating point number
                 (hardware dependent)
                 
   The only sure thing anbout the size

   short <= int <= long <= long long
   float <= double<= long double

  */

  cout << "Given me a and b " ;
  // cin is the object dedicated to the input
  cin >> a >> b ;
  
  int c = a + b ;
  
  cout << a << " + " << b << " = " << c << "\n" ;

  return 0 ;
}

3.cpp

File: 3.cpp

#include <iostream>
// I/O standard definitions

#include <iomanip>
// I/O manipulators

using namespace std ; // load the namespace std to ::

int
main () {
  double a, b ; // define the integer a and b
  
  cout << "Given me a and b " ;
  // cin is the object dedicated to the input
  cin >> a >> b ;
  
  double c = a + b ;
  
  // set(int n) = reserve n character for the next output 
  cout << setw(20) << a << " + " << b << " = " << c << "\n" ;
  
  // other manipulators
  /*
   setprecision( int ) = set precision in printing number
   setbase( int )      = set base to print integer 10, 8, 16
   setfill( char )     = set the filling character
   dec                 = print decimal
   hex                 = print hexadecimal
   oct                 = print octal
   
   flush               = flush buffer
   
   setw(int)           = reserve spaces for next output
   
   left                = align left
   right               = align right
   
   scientific          = print number scientific e.g. 1.234E23
   floatfield          = print number e.g. 123.343
   
   */
  
  cout << setw(20) << setfill('*') << a << " + " << b << " = " << c << "\n" ;
  cout << setprecision(10) << a << " + " << b << " = " << c << "\n" ;
  cout << a << " + " << b << " = " << c << "\n" ;
  
  cout << setbase(10) << 123 << " " << setbase(8) << 123 << "\n" ;
  cout << setbase(10) << 123 << " " << setbase(16) << 123 << "\n" ;

  cout << dec << 123 << " " << oct << 123 << "\n" ;
  cout << dec << 123 << " " << hex << 123 << "\n" ;
  
  return 0 ;
}

4.cpp

File: 4.cpp

#include <iostream>

using namespace std ;

/*
 Conditional and loop
 */

int
main () {
  /*
    if ( condition ) statement ;

    if ( condition ) statement ;
    else             statement ;    
    
   */
  int a = 1 ;
  int b = 2 ;
  if ( a < b ) cout << "a < b is true \n" ;
  else         cout << "a < b is true \n" ;
  
  if ( a < b ) {
    int c = a+b ;
    cout << "a < b is true \n" ;
  }
  else {
    cout << "a < b is true \n" ;
  }
  
  /*
    while ( condition ) statement ;
   */

  int i = 0 ;
  while ( ++i < 100 ) cout << i << endl ;
  
  /* ++i means i = i+1 ;
     there is also --i that means i = i-1 
     ++i and i++ are different in the sense that
     int j = i++ ; means that j = i next i = i+1 ;
     int j = ++i ; means that i = i+1 next j = i ;
  */
  
  /*
    do {
      instructions ;
    } while ( condition ) ;
  */
  
  i = 0 ;
  do {
    i++ ;
    cout << i << ", " ;
  } while ( i < 50 ) ;
  cout << endl ;
  
  /*
    for ( init ; condition ; next ) statement ;
    
    init      = is a statement executed ONE time at the beginning
    condition = is tested at EACH iteration
    next      = is a statement executed at the end of EACH iteration
  */
  // int i = 1 ; is an error
  int sum = 0 ;
  for ( int i = 1 ; // init
        i < 100 ;
        i++ ) {
    sum += i ; // sum = sum + i ;
  }
  cout << "sum = " << sum << endl ;
  /*
    short cut form of =
    
    a OP= b ; means   a = a OP b ;
    
    for example
   
       a /= b ; ==> a = a / b ;
       a *= b ; ==> a = a * b ;
       a |= b ; ==> a = a | b ;
    
   */
   
  /*
    The switch statement
    
    switch ( int ) {
      case number:
      
      break ; // interrupt the switch,
              // if not present continue to the next case.
      case number:
      
      case number:
      
      default:
    }
  */

  int j = -1 ;
  switch (j) {
    case 7:   cout << "j = 7" << endl ;
    case -1:  cout << "j = -1" << endl ;
    case 0:   cout << "j = 0" << endl ;
    case 3:   cout << "j = 3" << endl ;
    break ;
    default:  cout << "j = none" << endl ;
  }
  return 0 ;
  
  /*
    break    instruction
    continue instruction
  
    for ( ) {
      ...
      ...
      if ( cond ) break ; // interrupt the loop
      ..
      ..
      if ( cond ) continue ; // skip to the next iteration
      ..
      ..
    }

  */
}

Lesson of 17/5/2005

5.cpp

File: 5.cpp

#include <iostream>

using namespace std ;

int
main()
{
  // comparison operator
  
  double a = 1 ;
  double b = 2 ;

  // equality
  if ( a == b ) { // = is the assign not the comparison sign
    cout << "a is equal to b\n" ;
  } else {
    cout << "a is NOT equal to b\n" ;  
  }
  
  // not equal
  if ( a != b ) {
    cout << "a is NOT equal to b\n" ;
  } else {
    cout << "a is equal to b\n" ;  
  }

  // greater 
  if ( a > b ) {
    cout << "a is greater than b\n" ;
  } else {
    cout << "a is NOT greater than b\n" ;  
  }

  // lesser
  if ( a < b ) {
    cout << "a is less than b\n" ;
  } else {
    cout << "a is NOT less than b\n" ;	
  }

  // greater or equal
  if ( a >= b ) { // => is not a comparison operator
    cout << "a is greater or equal than b\n" ;
  } else {
    cout << "a is NOT greater or equal than b\n" ;	
  }

  // lesser or equal
  if ( a <= b ) { // =< is not a comparison operator
    cout << "a is less or equal than b\n" ;
  } else {
    cout << "a is NOT less or equal than b\n" ;	 
  }

  // BOOLEAN
  double c = 3 ;
  if ( a == b || b == c ) ; // || is the OR operator
  if ( a == b && b == c ) ; // && is the AND operator

  // using boolean type
  bool eq  = a == b ; // eq takes the values "true" or "false"
  bool neq = b != c ;
  if ( eq || neq ) ;
  
  bool bc = true ;
  bool ac = false ;
  
  // loop with the use of boolean
  bool cont = true ;
  int i = 1 ;
  while ( cont ) {
    cout << i << endl ;
    cont = i < 100 ;
    ++i ;
  }
  
  // boolean arithmetic
  {
	bool a, b, c, d, e ;
    a = true ;
    b = true ;
    c = false ;
    d = false ;
    e = false ;
	
    a = (b && c) || e ; // (b AND c) OR e
    b = ! c ;			// negation
  }
  
  return 0 ;
}

6.cpp

File: 6.cpp

#include <iostream>
#include <cmath>

/*
  most important include
  
  iostream
  fstream
  iomanip
  cmath
  cstblib

  from standard template library

  string
  vector
  
*/

using namespace std ;

typedef long double value_type ; 
// make value_type an alias of long double

typedef int index_type ; 
// make index_type an alias of int

int
main()
{
  // example the use the trapezoidal rule to 
  // compute the integral of f(x)=1+sin(x)
  // using 100 intervals
  
  value_type a = 1 ;
  value_type b = 10 ;
  index_type n = 100 ;
  
  value_type fa = (1+sin(a)) ;
  value_type fb = (1+sin(b)) ;
  value_type h  = (b-a) / n ;
  
  value_type res = (fa+fb)*h/2 ;
  for ( index_type i = 1 ; i < n ; ++i ) {
    value_type xi = a + i*h ;
    value_type fi = (1+sin(xi)) ;
    res += h * fi ; 
  }

  cout << "The approximate integral is = " << res << "\n" ;
  return 0 ;
}

7.cpp

File: 7.cpp

#include <iostream>
#include <cmath>

using namespace std ;

typedef long double value_type ; 
// make value_type an alias of long double

typedef int index_type ; 
// make index_type an alias of int

// declaration of fun
value_type fun( value_type x ) ;

int
main()
{
  // example the use the trapezoidal rule to 
  // compute the integral of f(x)=1+sin(x)
  // using 100 intervals
  
  value_type a = 1 ;
  value_type b = 10 ;
  index_type n = 100 ;
  
  value_type fa = fun(a) ;
  value_type fb = fun(b) ;
  value_type h  = (b-a) / n ;
  
  value_type res = (fa+fb)*h/2 ;
  for ( index_type i = 1 ; i < n ; ++i ) {
    value_type xi = a + i*h ;
    value_type fi = fun(xi) ;
    res += h * fi ; 
  }

  cout << "The approximate integral is = " << res << "\n" ;
  return 0 ;
}

// definition of fun
value_type // return value of the function
fun( value_type x ) { // fun is the name of the function
  return 1+sin(x);    // x is one argument
}

8.1.cpp

File: 8.1.cpp

#include <iostream>
#include <cmath>

using namespace std ;

typedef long double value_type ; 
// make value_type an alias of long double

typedef int index_type ; 
// make index_type an alias of int

// declaration of fun
value_type fun( value_type x ) ;

int
main()
{
  // example the use the trapezoidal rule to 
  // compute the integral of f(x)=1+sin(x)
  // using 100 intervals
  
  value_type a = 1 ;
  value_type b = 10 ;
  index_type n = 100 ;
  
  value_type fa = fun(a) ;
  value_type fb = fun(b) ;
  value_type h  = (b-a) / n ;
  
  value_type res = (fa+fb)*h/2 ;
  for ( index_type i = 1 ; i < n ; ++i ) {
    value_type xi = a + i*h ;
    value_type fi = fun(xi) ;
    res += h * fi ; 
  }

  cout << "The approximate integral is = " << res << "\n" ;
  return 0 ;
}

8.2.cpp

File: 8.2.cpp

#include <iostream>
#include <cmath>

using namespace std ;

typedef long double value_type ; 
// make value_type an alias of long double

typedef int index_type ; 
// make index_type an alias of int

// declaration of fun
value_type fun( value_type x ) ;

// definition of fun
value_type // return value of the function
fun( value_type x ) { // fun is the name of the function
  return 2+sin(x);    // x is one argument
}

9.1.cpp

File: 9.1.cpp

#include "9.h"

int
main()
{
  // example the use the trapezoidal rule to 
  // compute the integral of f(x)=1+sin(x)
  // using 100 intervals
  
  value_type a = 1 ;
  value_type b = 10 ;
  index_type n = 100 ;
  
  value_type fa = fun(a) ;
  value_type fb = fun(b) ;
  value_type h  = (b-a) / n ;
  
  value_type res = (fa+fb)*h/2 ;
  for ( index_type i = 1 ; i < n ; ++i ) {
    value_type xi = a + i*h ;
    value_type fi = fun(xi) ;
    res += h * fi ; 
  }

  cout << "The approximate integral is = " << res << "\n" ;
  return 0 ;
}

9.1.cpp

File: 9.2.cpp

#include "9.h"

// definition of fun
value_type // return value of the function
fun( value_type x ) { // fun is the name of the function
  return 2+sin(x);    // x is one argument
}

9.h

File: 9.h

#include <iostream>
#include <cmath>

using namespace std ;

typedef double value_type ; 
// make value_type an alias of double

typedef int index_type ; 
// make index_type an alias of int

// declaration of fun
value_type fun( value_type x ) ;

10.1.cpp

File: 10.1.cpp

#include "10.h"
#include <iomanip>

value_type valabs( value_type x )   ;

int
main() {
  //  solving f(x) = 0 by Newton method
  value_type x       = 0 ;
  value_type epsilon = 1E-6 ;
  for ( index_type i = 1 ; i < 100 ; ++i ) {
    value_type f = fun(x) ;
    if ( valabs(f) < epsilon ) break ; // interrupt the loop
    value_type f_1 = fun_1(x) ;
    x = x - f/f_1   ;
    cout << "iter = " << setw(2) << i
         << " x = " << setw(10) << x
         << " residual = " << f << "\n" ;
  }   
  cout << "result x = " << x << "\n" ;
}

value_type
valabs( value_type x ) {
  return x > 0 ? x : -x ;
}

10.2.cpp

File: 10.2.cpp

#include "10.h"

value_type
fun(value_type x) {
 return pow(x-1,2) ; // (x-1)*(x-1)
}

value_type
fun_1(value_type x) {
 return 2*(x-1) ;
}

10.h

File: 10.h

#include <iostream>
#include <cmath>

using namespace std ;

typedef double value_type  ; 
// make value_type an alias of double

typedef int index_type ;
// make index_type an alias of int

// declaration of fun
value_type fun( value_type x ) ;

// declaration of  fun_1
value_type fun_1( value_type x ) ;

1.cpp

File: 11.cpp

#include "11.h"
#include <iomanip>

// include again 11.h
#include "11.h"

value_type valabs( value_type x ) ;

int
main () {
  // solving f(x) = 0 by newton method
  value_type x       = 0 ;
  value_type epsilon = 1E-6 ;
  for ( index_type i = 1 ; i < 100 ; ++i ) {
    value_type f   = fun(x) ;
    if ( valabs( f ) < epsilon ) break ;
    value_type f_1 = fun_1(x) ;
    x = x - f/f_1 ;
    cout << "iter = " << setw(2) << i
         << " x = " << setw(10) << x
         << " residual = " << f << "\n" ;
  } 
  cout << " result x = " << x << "\n" ;
}

11.h

File: 11.h

#ifndef ELEVEN_H
#define ELEVEN_H

#include <iostream>
#include <cmath>

using namespace std ;

typedef double value_type ; 
// make value_type an alias of double

typedef int index_type ; 
// make index_type an alias of int

inline
value_type
fun(value_type x) {
  return pow(x-1,2) ; // (x-1)*(x-1)
}

inline
value_type
fun_1(value_type x) {
  return 2*(x-1) ;
}

inline
value_type
valabs( value_type x ) {
  return x > 0 ? x : -x ;
}

#endif

12.cpp

File: 12.cpp

#include <iomanip>
#include <string>

using namespace std ;

/*
  Introduction to struct

  want to collect
  
  first name
  second name
  age
  nationality
  weight 
  hair
  driving licence
  
*/

typedef double value_type ;
typedef int    index_type ;

struct record {
  string     first_name ;
  string     second_name ;
  index_type age ;
  string     nationality ;
  value_type weight ;
  string     hair_color ;
  bool       have_a_driving_licence ;
}  ;

int
main () {
  struct record a ;
  
  a . first_name = "Pippo" ;
  a . second_name = "Pluto" ;
  a . age = 1204 ;
  a . have_a_driving_licence = true ;
  
  // Define and initialize
  struct record b = {
    "Pluto",
    "Dog",
    1204,
    "Topolinia",
    256,
    "Red",
    true
  } ;
  
  typedef struct record Record ;
  
  Record c ;
  
}

Lesson of 18/5/2005

13.cpp

File: 13.cpp

#include <iostream>
#include <iomanip>

using namespace std ;

typedef double value_type ;
typedef int    index_type ;

// vec is a pointer to memory allocated
// with 100 value_type

value_type vec[100] ;

int
main() {
  /*
    to access the elements of vec we use
    the operator [ ], and the index run
    from 0 to dim-1 (here dim = 100)
   */
   
  vec[1]  ; // is the second element
  vec[99] ; // is the last element

  for ( index_type i = 0 ; i < 100 ; ++i ) 
    vec[i] = i * i ;

  value_type * p_vec ; // is a pointer to value_type things
  p_vec = vec ; // p_vec point to the first value of
                // the vector vec
                
  // NO vec = p_vec ; vec is not assignable
  
  p_vec = & vec[0] ; // is equivalent to p_vec = vec 

  cout << "& vec[0] = " << &vec[0] << " ==> " 
       << (long)&vec[0]<< "\n"
       << " p_vec   = " << p_vec << " ==> " 
       << (long)p_vec << "\n" ;
       
  // (long) is a "casting" operator, reinterprets
  // p_vec (address) as a long integer.
       
  cout << "vec[0] = " << vec[0] << "\n"
       << "*p_vec = " << *p_vec << "\n" ;
  // * is the deference operator:
  // p_vec contains the address of a value_type value 
  // *p_vec is the value_type value at the address p_vec
  
  
  /*
    pointer ca be assigned, incremented,
    subtracted, ...
    some example
   */
   
  p_vec = & vec[11] ; // & estract the address from an lvalue
  
  cout << "vec[11] = " << vec[11] << "\n"
       << "*p_vec  = " << *p_vec << "\n" ;

  cout << "& vec[11] = " << &vec[11] << "\n"
       << " p_vec    = " << p_vec << "\n" ;
       
 /*
   Equivalent access with the * operator
  */
  
  // vec[13] is equalent to *(vec+13)
  cout << "vec[13] = " << vec[13]
       << "  *(vec+13) = " << *(vec+13) << "\n" ;
       
  //the address vec+13 is:
  // the address of vec + 13 * sizeof(value_type) bytes

  // assign a pointer to a translated value
  p_vec = vec + 14 ;
  cout << "vec[14] = " << vec[14]
       << "  *p_vec = " << *p_vec << "\n" ;

  /*
    Operation with pointers
  */
  
  // loop with pointers

  // old loop
  //for ( index_type i = 0 ; i < 100 ; ++i ) 
  //  vec[i] = 0 ;
  for ( p_vec = vec ; p_vec < vec + 100 ; ++p_vec )
    *p_vec = 0 ;
  
  // another way
  for ( p_vec = vec ; p_vec < vec + 100 ; )
    *p_vec++ = 0 ;
  
  // another way again
  p_vec = vec ;
  while ( p_vec < vec + 100 ) *p_vec++ = 0 ;
  
  // distance between pointer
  p_vec = &vec[11] ;
  cout << "p_vec - vec = " << p_vec - vec << "\n" ;
  cout << "(short*) p_vec - (short*) vec = " 
       << (short*) p_vec - (short*) vec << "\n" ;
  
  // the sizeof
  cout << "sizeof(value_type) = " << sizeof(value_type) << "\n"
       << "sizeof(short)      = " << sizeof(short) << "\n" ;
}

14.cpp

File: 14.cpp

#include <iostream>
#include <iomanip>

using namespace std ;

typedef double value_type ;
typedef int    index_type ;

// vec is a pointer to memory allocated
// with 100 value_type

// declare and initialize
value_type avec[10] = { 1,2,3,4,5,6,7,8,9,10 } ;

// declare and initialize, infer the size by the 
// initalization list
value_type bvec[]   = { 1,2,3,4,5,6,7,8,9,10 } ;

struct point {
  value_type x ;
  value_type y ;
  value_type z ;
} ;

typedef struct point Point ;

// collapsing the previous two definition as

typedef struct {
  value_type x ;
  value_type y ;
  value_type z ;
} point3d ;

int
main() {
  value_type * cvec ;
  
  // use dynamic allocation
  cvec = new value_type [10] ;
  // new is the allocator operator
  // value_type is the type I am allocating
  // [10] is the number of cell allocated by new
  
  // YES cvec = bvec + 100 ;
  // NO  bvec = cvec + 100 ;
  
  struct point * old_pnts ;
  point3d * pnts ;
  
  old_pnts = new struct point [100] ;
  old_pnts = new Point [100] ;
  pnts     = new point3d [100] ;
  
  // access to the x of the second element:
  pnts[1] . x ; // first get the first element then 
                // access by . operator
  (*(pnts+1)) . x ;
  (pnts+1) -> x ;
  
  for ( index_type i = 0 ; i < 100 ; ++i ) {
    pnts[i] . x = 1 ;
    pnts[i] . y = i ;
    pnts[i] . z = -3 ;
    // another equivalent possibility
    point3d * p = pnts + i ;
    p -> x = 1 ;
    p -> y = i ;
    p -> z = -3 ;
  }
  
  for ( index_type i = 0 ; i < 10 ; ++i ) 
    cvec[i] = avec[i] + bvec[i] ;
  
  // freeing memory operator is delete:
  delete [] cvec ;
  // delete is the deletion operator
  // [] mean that I am deallocating a vector
  // cvec is the idetifier of the pointer allocated
}

15.cpp

File: 15.cpp

#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std ;

typedef double value_type ;
typedef int    index_type ;

typedef struct {
  value_type x, y, z ;
} point3d ;

value_type
norm( point3d const p ) { // p is passed by VALUE!
  return sqrt( p.x*p.x +
               p.y*p.y +
               p.z*p.z ) ;
}

value_type
norm_addr( point3d const * p_addr ) { // p is passed by ADDRESS!
  p_addr = p_addr + 1 ;
  p_addr = p_addr - 1 ;
  return sqrt( p_addr -> x * p_addr -> x +
               p_addr -> y * p_addr -> y +
               p_addr -> z * p_addr -> z ) ;
}

value_type
norm_addr_1( point3d const * const p_addr ) { // p is passed by ADDRESS!
  // NO p_addr = p_addr + 1 ; pointer cannot be modified
  // NO p_addr = p_addr - 1 ;
  return sqrt( p_addr -> x * p_addr -> x +
               p_addr -> y * p_addr -> y +
               p_addr -> z * p_addr -> z ) ;
}

value_type
norm_addr_2( point3d * const p_addr ) { // p is passed by ADDRESS!
  // NO  p_addr = p_addr + 1 ; pointer cannot be modified
  // NO  p_addr = p_addr - 1 ;
  // YES p_addr -> x = 2 ;     thing pointed can be modified
  return sqrt( p_addr -> x * p_addr -> x +
               p_addr -> y * p_addr -> y +
               p_addr -> z * p_addr -> z ) ;
}

value_type
norm_ref( point3d const & p_ref ) { // p is passed by REFERENCE!
  return sqrt( p_ref . x * p_ref . x +
               p_ref . y * p_ref . y +
               p_ref . z * p_ref . z ) ;
}

// modification of p_ref
value_type
norm_ref_modify( point3d & p_ref ) { // p is passed by REFERENCE!
  value_type res = sqrt( p_ref . x * p_ref . x +
                         p_ref . y * p_ref . y +
                         p_ref . z * p_ref . z ) ;
  p_ref . z = 0 ;
  return res ;
}


int
main() {
  point3d pt ;
  
  pt . x = 1 ; 
  pt . y = 2 ;
  pt . z = 3 ;
  
  cout << " norm_ref_modify  = " << norm_ref_modify( pt ) << "\n" ;
  cout << " norm      = " << norm( pt ) << "\n" ;
  cout << " norm_addr = " << norm_addr( & pt ) << "\n" ;
  cout << " norm_ref  = " << norm_ref( pt ) << "\n" ;
  
  point3d volatile * p_pt = & pt ;
  // some operation
  // you cannot assume that pt is unchanged
   
}

Lesson of 23/5/2005

16.1.cpp

File: 16.1.cpp

#include "16.h"

namespace newton_solver_2d {

  // F(x) = / x0 - x1           \
  //        |                   |
  //        \ x0**2 + x1**3 - 2 /
  void
  F( vec2 const & x, vec2 & Fx ) {
    Fx[0] = x[0] - x[1] ;
    Fx[1] = x[0]*x[0] - x[1]*x[1]*x[1] - 2 ;
    // Fx[1] = pow(x[0],2) - pow(x[1],3) -2 ;
  }
  
  // JFx[row][column]
  void JF( vec2 const & x, mat2 & JFx ) {
    JFx[0][0] = 1            ; // DF[0]/Dx[0]
    JFx[0][1] = -1           ; // DF[0]/Dx[1]
    JFx[1][0] = 2*x[0]       ; // DF[1]/Dx[0]
    JFx[1][1] = -3*x[1]*x[1] ; // DF[1]/Dx[1]
  }
  
  /*
   Small 2x2 linear solver
   */
  
  bool
  linear_solver( mat2 const & A, // coefficients matrix
                 vec2 const & b, // r.h.s
                 vec2       & x  // the solution of the linear system
               ) {
    // Using Cramer rule
    
    value_type det = A[0][0] * A[1][1] - A[1][0] * A[0][1] ;
    if ( det == 0 ) return false ;
    
    x[0] = ( b[0] * A[1][1] - b[1] * A[0][1] ) / det ;
    x[1] = ( A[0][0] * b[1] - A[1][0] * b[0] ) / det ;

    return true ;
  }

  /*
     Newton scheme
                         -1
     x    = x   - JF(x  )   F(x )
      n+1    n        n        n
   
     the implementation:
    
       SOLVE: JF(x ) h = F(x )
                  n         n

       UPDATE:  x   = x + h
                n+1    n
   */
  bool
  solver( vec2 const & x0,
          F_type     * pF,
          JF_type    * pJF,
          index_type const n,
          value_type const & tol,
          vec2 & x
          ) {
    
    // set x to the initial value
    x[0] = x0[0] ;
    x[1] = x0[1] ;
    
    for ( index_type i = 0 ; i < n ; ++i ) {
      vec2 f, h ;
      mat2 jf ;
      // perform a Newton step 
      // substep 1: compute F(x)
      (*pF)(x,f) ; // shortcut borrowed for C syntax ==> pF(x,f) ;
      // substep 2: compute JF(x)
      (*pJF)(x,jf) ; // shortcut borrowed for C syntax ==> pJF(x,jf) ;
      // substep 3: compute the correction h = JF(x)^(-1)F(x)
      bool ok = linear_solver(jf, f, h) ;
      if ( !ok ) return false ;
      // substep 4: update x 
      x[0] -= h[0] ;
      x[1] -= h[1] ;
      // check for the stopping criteria
      value_type res = sqrt( h[0]*h[0]+h[1]*h[1] ) ;
      if ( res <= tol ) return true ; 
    }
    return false ;
  }
}

16.2.cpp

File: 16.2.cpp

#include "16.h"

int
main() {
  using newton_solver_2d::vec2 ;
  using newton_solver_2d::F ;
  using newton_solver_2d::JF ;
  using newton_solver_2d::solver ;
  using namespace std ;

  vec2 x0, x ;
  
  x0[0] = 1 ;
  x0[1] = 1 ;
  solver( x0, F, JF, 100, 1E-6, x ) ; 
  
  cout << "The solution is x = [ " 
       << x[0] << ", " << x[1] << "]\n" ;
}

16.h

File: 16.h

// standard trick for only one time inclusion
#ifndef SIXTEEN_H
#define SIXTEEN_H

#include <iostream>
#include <iomanip>
#include <cmath>

namespace newton_solver_2d {
  using namespace std ;
  
  // set the value_type to double,
  // value_type is the generic name for "number"
  // (non integer)type
  typedef double value_type ;

  // set the index_type to int,
  // index_type is the generic name for "integer" type
  typedef int index_type ;
  
  // define vec2 as a name for the vector with 2 components
  typedef value_type vec2[2] ; // value_type a[2] ; ==> vec2 a ;
  
  // define mat2 as a name for the matrix 2x2
  typedef value_type mat2[2][2] ; // value_type a[2][2] ; ==> mat2 a ;
 
  // an equivalent definition
  // typedef vec2 mat2[2] ;
  
  /*
   now we set the definition (prototipe)
   of the functions
   */
  
  /*
    This methos is not correct:
    vec2
    fun( vec2 x ) {
      vec2 res ;
      res[0] = .... ;
      res[1] = .... ;
      return res ;
    }
   
    return res ; // return the pointer to a temporary
  */
  
  // prototype of the non-linear system 
  extern void F( vec2 const & x, vec2 & Fx ) ;

  // prototype of the Jacobian of the non-linear system
  extern void JF( vec2 const & x, mat2 & JFx ) ;
  
  /*
    extern = means that the bodies of the functions
             are somewhere else.
   */

  // prototype for the solver
  
  extern
  bool // return true = if the solution was found
  solver( vec2 const & x0,
          // the starting point of the Newton scheme
          void (*pF)  ( vec2 const & x, vec2 & Fx ),
          // pF is a pointer to the function
          void (*pJF) ( vec2 const & x, mat2 & JFx ),
          // pF is a pointer to the jacobian function
          index_type const n,
          // n is the maximum number of iteration allowed
          value_type const & tol,
          // tolerance for the stopping critertia
          vec2 & x
          // x contains the computed solution
         ) ;
  
  /*
    A less difficult to read definition can be 
    obtained by splitting the prototype definition
    as follows
   */
  
  typedef void F_type  ( vec2 const &, vec2 &) ;
  typedef void JF_type ( vec2 const &, mat2 &) ;

  extern
  bool
  solver( vec2 const & x0,
          F_type     * pF,
          JF_type    * pJF,
          index_type const n,
          value_type const & tol,
          vec2 & x
         ) ;
}

#endif

17.1.cpp

File: 17.1.cpp

#include "17.h"

namespace newton_solver_2d {

  // F(x) = / x0 - x1           \
  //        |                   |
  //        \ x0**2 + x1**3 - 2 /
  vec2
  F( vec2 const & in ) {
    vec2 Fx ; 
    Fx.c1 = in.c1 - in.c2 ;
    Fx.c2 = in.c1*in.c1 - in.c2*in.c2*in.c2 - 2 ;
    return Fx ;
  }
  
  //
  mat2x2
  JF( vec2 const & in ) {
    mat2x2 JFx ;
    JFx.A11 = 1              ; // DF[0]/Dx[0]
    JFx.A12 = -1             ; // DF[0]/Dx[1]
    JFx.A21 = 2*in.c1        ; // DF[1]/Dx[0]
    JFx.A22 = -3*in.c2*in.c2 ; // DF[1]/Dx[1]
    return JFx ;
  }
  
  /*
   Small 2x2 linear solver
   */
  
  bool
  linear_solver( mat2x2 const & A, // coefficients matrix
                 vec2   const & b, // r.h.s
                 vec2         & x  // the solution of the linear system
               ) {
    // Using Cramer rule
    
    value_type det = A.A11 * A.A22 - A.A21 * A.A12 ;
    if ( det == 0 ) return false ;
    
    x.c1 = ( b.c1 * A.A22 - b.c2 * A.A12 ) / det ;
    x.c2 = ( A.A11 * b.c2 - A.A21 * b.c1 ) / det ;

    return true ;
  }

  /*
     Newton scheme
                         -1
     x    = x   - JF(x  )   F(x )
      n+1    n        n        n
   
     the implementation:
    
       SOLVE: JF(x ) h = F(x )
                  n         n

       UPDATE:  x   = x + h
                n+1    n
   */
  bool
  solver( vec2 const & x0,
          F_type     * pF,
          JF_type    * pJF,
          index_type const n,
          value_type const & tol,
          vec2 & x
          ) {
    
    // set x to the initial value
    x = x0 ; // not reccomended
    
    for ( index_type i = 0 ; i < n ; ++i ) {
      vec2 f, h ;
      mat2x2 jf ;
      // perform a Newton step 
      // substep 1: compute F(x)
      f = (*pF)(x) ; // shortcut borrowed for C syntax ==> pF(x) ;
      // substep 2: compute JF(x)
      jf = (*pJF)(x) ; // shortcut borrowed for C syntax ==> pJF(x) ;
      // substep 3: compute the correction h = JF(x)^(-1)F(x)
      bool ok = linear_solver(jf, f, h) ;
      if ( !ok ) return false ;
      // substep 4: update x 
      x.c1 -= h.c1 ;
      x.c2 -= h.c2 ;
      // check for the stopping criteria
      value_type res = sqrt( h.c1*h.c1+h.c2*h.c2 ) ;
      if ( res <= tol ) return true ; 
    }
    return false ;
  }
}

17.2.cpp

File: 17.2.cpp

#include "17.h"

int
main() {
  using newton_solver_2d::vec2 ;
  using newton_solver_2d::F ;
  using newton_solver_2d::JF ;
  using newton_solver_2d::solver ;
  using namespace std ;

  vec2 x0, x ;
  
  x0.c1 = 1 ;
  x0.c2 = 1 ;
  solver( x0, F, JF, 100, 1E-6, x ) ; 
  
  cout << "The solution is x = [ " 
       << x.c1 << ", " << x.c2 << "]\n" ;
}

17.h

File: 17.h

// standard trick for only one time inclusion
#ifndef SEVENTEEN_H
#define SEVENTEEN_H

#include <iostream>
#include <iomanip>
#include <cmath>

namespace newton_solver_2d {
  using namespace std ;
  
  // set the value_type to double,
  // value_type is the generic name for "number"
  // (non integer)type
  typedef double value_type ;

  // set the index_type to int,
  // index_type is the generic name for "integer" type
  typedef int index_type ;
  
  typedef struct {
    value_type c1 ; // component one
    value_type c2 ; // component two
  } vec2 ;

  typedef struct {
    value_type A11 ;
    value_type A12 ;
    value_type A21 ;
    value_type A22 ;
  } mat2x2 ;
 
  // an equivalent definition
  // typedef vec2 mat2[2] ;
  
  /*
   now we set the definition (prototipe)
   of the functions
   */
    
  // prototype of the non-linear system 
  extern vec2 F( vec2 const & in ) ;

  // prototype of the Jacobian of the non-linear system
  extern mat2x2 JF( vec2 const & in ) ;
  
  // prototype for the solver
  
  extern
  bool // return true = if the solution was found
  solver( vec2 const & x0,
          // the starting point of the Newton scheme
          vec2   (*pF)  ( vec2 const & x ),
          // pF is a pointer to the function
          mat2x2 (*pJF) ( vec2 const & x ),
          // pF is a pointer to the jacobian function
          index_type const n,
          // n is the maximum number of iteration allowed
          value_type const & tol,
          // tolerance for the stopping critertia
          vec2 & x
          // x contains the computed solution
         ) ;
  
  /*
    A less difficult to read definition can be 
    obtained by splitting the prototype definition
    as follows
   */
  
  typedef vec2   F_type  ( vec2 const & ) ;
  typedef mat2x2 JF_type ( vec2 const & ) ;

  extern
  bool
  solver( vec2 const & x0,
          F_type     * pF,
          JF_type    * pJF,
          index_type const n,
          value_type const & tol,
          vec2 & x
         ) ;
}

#endif

Lesson of 24/5/2005

18.1.cpp

File: 18.1.cpp

#include "18.h"

namespace newton_solver {
  
  /*
   Small 2x2 linear solver
   */
  
  bool
  linear_solver( mat const & A, // coefficients matrix
                 vec const & b, // r.h.s
                 vec       & x  // the solution of the linear system
               ) {
    // Using Cramer rule
    value_type det = A.values[0][0] * A.values[1][1]
                   - A.values[1][0] * A.values[0][1] ;
    if ( det == 0 ) return false ;
    
    x.values[0] = ( b.values[0] * A.values[1][1] 
                  - b.values[1] * A.values[0][1] ) / det ;
    x.values[1] = ( A.values[0][0] * b.values[1]
                  - A.values[1][0] * b.values[0] ) / det ;

    return true ;
  }

  /*
     Newton scheme
                         -1
     x    = x   - JF(x  )   F(x )
      n+1    n        n        n
   
     the implementation:
    
       SOLVE: JF(x ) h = F(x )
                  n         n

       UPDATE:  x   = x + h
                n+1    n
   */
  bool
  solver( vec        const & x0,
          F_type           * pF,
          JF_type          * pJF,
          index_type const   n,
          value_type const & tol,
          vec              & x
          ) {
    
    // set x to the initial value
    x = x0 ; // not reccomended
    
    for ( index_type i = 0 ; i < n ; ++i ) {
      vec f, h ;
      mat jf ;
      // perform a Newton step 
      // substep 1: compute F(x)
      f = pF(x) ;
      // substep 2: compute JF(x)
      jf = pJF(x) ;
      // substep 3: compute the correction h = JF(x)^(-1)F(x)
      bool ok = linear_solver(jf, f, h) ;
      if ( !ok ) return false ;
      // substep 4: update x 
      for ( index_type i = 0 ; i < dim ; ++i )
        x.values[i] -= h.values[i] ;
      // check for the stopping criteria
      value_type res = 0 ;
      for ( index_type i = 0 ; i < dim ; ++i )
        res += h.values[i]*h.values[i] ;
      res = sqrt(res) ;
      if ( res <= tol ) return true ; 
    }
    return false ;
  }
}

18.2.cpp

File: 18.2.cpp

#include "18.h"

namespace newton_solver {
  // F(x) = / x0 - x1           \
  //        |                   |
  //        \ x0**2 + x1**3 - 2 /
  vec
  F( vec const & in ) {
    vec Fx ;
    value_type const * x = in . values ;
    Fx.values[0] = x[0] - x[1] ; // in . values[0] - in . values[1]
    Fx.values[1] = x[0]*x[0] - x[1]*x[1]*x[1] - 2 ;
    return Fx ;
  }

  //
  mat
  JF( vec const & in ) {
    mat JFx ;
    value_type const * x = in . values ;
    JFx.values[0][0] = 1            ; // DF[0]/Dx[0]
    JFx.values[0][1] = -1           ; // DF[0]/Dx[1]
    JFx.values[1][0] = 2*x[0]       ; // DF[1]/Dx[0]
    JFx.values[1][1] = -3*x[1]*x[1] ; // DF[1]/Dx[1]
    return JFx ;
  }
}

int
main() {
  using newton_solver::vec ;
  using newton_solver::F ;
  using newton_solver::JF ;
  using newton_solver::solver ;
  using namespace std ;

  vec x0, x ;
  
  x0.values[0] = 1 ;
  x0.values[1] = 1 ;
  solver( x0, F, JF, 100, 1E-6, x ) ; 
  
  cout << "The solution is x = [ " 
       << x.values[0] << ", "
       << x.values[1] << "]\n" ;
}

18.h

File: 18.h

// standard trick for only one time inclusion
#ifndef EIGHTEEN_H
#define EIGHTEEN_H

#include <iostream>
#include <iomanip>
#include <cmath>

namespace newton_solver {
  using namespace std ;
  
  // set the value_type to double,
  // value_type is the generic name for "number"
  // (non integer)type
  typedef double value_type ;

  // set the index_type to int,
  // index_type is the generic name for "integer" type
  typedef int index_type ;
  
  index_type const dim = 2 ;
  
  typedef struct {
    value_type values[dim] ;
  } vec ;

  typedef struct {
    value_type values[dim][dim];
  } mat ;
 
  /*
   now we set the definition (prototipe)
   of the functions
   */
    
  // prototype of the non-linear system 
  extern vec F( vec const & in ) ;

  // prototype of the Jacobian of the non-linear system
  extern mat JF( vec const & in ) ;
  
  // prototype for the solver
    
  typedef vec F_type  ( vec const & ) ;
  typedef mat JF_type ( vec const & ) ;

  extern
  bool
  solver( vec        const & x0,
          F_type           * pF,
          JF_type          * pJF,
          index_type const   n,
          value_type const & tol,
          vec              & x
         ) ;
}

#endif

19.1.cpp

File: 19.1.cpp

#include "19.h"

namespace newton_solver {
  
  /*
   Small 2x2 linear solver
   */
  
  bool
  linear_solver( mat const & A, // coefficients matrix
                 vec const & b, // r.h.s
                 vec       & x  // the solution of the linear system
               ) {
    // Using Cramer rule
    value_type det = A(0,0) * A(1,1) - A(1,0) * A(0,1) ;
    if ( det == 0 ) return false ;
    
    x[0] = ( b[0] * A(1,1) - b[1] * A(0,1) ) / det ;
    x[1] = ( A(0,0) * b[1] - A(1,0) * b[0] ) / det ;

    return true ;
  }

  /*
     Newton scheme
                         -1
     x    = x   - JF(x  )   F(x )
      n+1    n        n        n
   
     the implementation:
    
       SOLVE: JF(x ) h = F(x )
                  n         n

       UPDATE:  x   = x + h
                n+1    n
   */
  bool
  solver( vec        const & x0,
          F_type           * pF,
          JF_type          * pJF,
          index_type const   n,
          value_type const & tol,
          vec              & x
          ) {
    
    // set x to the initial value
    x = x0 ; // not reccomended
    
    for ( index_type i = 0 ; i < n ; ++i ) {
      vec f, h ;
      mat jf ;
      // perform a Newton step 
      // substep 1: compute F(x)
      f = pF(x) ;
      // substep 2: compute JF(x)
      jf = pJF(x) ;
      // substep 3: compute the correction h = JF(x)^(-1)F(x)
      bool ok = linear_solver(jf, f, h) ;
      if ( !ok ) return false ;
      // substep 4: update x 
      for ( index_type i = 0 ; i < dim ; ++i )
        x[i] -= h[i] ;
      // check for the stopping criteria
      value_type res = 0 ;
      for ( index_type i = 0 ; i < dim ; ++i )
        res += h[i]*h[i] ;
      res = sqrt(res) ;
      if ( res <= tol ) return true ; 
    }
    return false ;
  }
}

19.2.cpp

File: 19.2.cpp

#include "19.h"

namespace newton_solver {
  // F(x) = / x0 - x1           \
  //        |                   |
  //        \ x0**2 + x1**3 - 2 /
  vec
  F( vec const & x ) {
    vec Fx ;
    Fx[0] = x[0] - x[1] ; // in . values[0] - in . values[1]
    Fx[1] = x[0]*x[0] - x[1]*x[1]*x[1] - 2 ;
    return Fx ;
  }

  //
  mat
  JF( vec const & in ) {
    mat JFx ;
    value_type const * x = in . values ;
    JFx(0,0) = 1            ; // DF[0]/Dx[0]
    JFx(0,1) = -1           ; // DF[0]/Dx[1]
    JFx(1,0) = 2*x[0]       ; // DF[1]/Dx[0]
    JFx(1,1) = -3*x[1]*x[1] ; // DF[1]/Dx[1]
    return JFx ;
  }
}

int
main() {
  using newton_solver::vec ;
  using newton_solver::F ;
  using newton_solver::JF ;
  using newton_solver::solver ;
  using namespace std ;

  vec x0, x ;
  
  x0.values[0] = 1 ;
  x0.values[1] = 1 ;
  solver( x0, F, JF, 100, 1E-6, x ) ; 
  
  cout << "The solution is x = [ " 
       << x.values[0] << ", "
       << x.values[1] << "]\n" ;
}

19.h

File: 19.h

// standard trick for only one time inclusion
#ifndef NINETEEN_H
#define NINETEEN_H

#include <iostream>
#include <iomanip>
#include <cmath>

namespace newton_solver {
  using namespace std ;
  
  // set the value_type to double,
  // value_type is the generic name for "number"
  // (non integer)type
  typedef double value_type ;

  // set the index_type to int,
  // index_type is the generic name for "integer" type
  typedef int index_type ;
  
  index_type const dim = 2 ;
  
  typedef struct {
    value_type values[dim] ;
    
    // declare the operator []
    value_type & operator [] ( index_type i ) 
    { return this -> values[i] ; }
    
    value_type const & operator [] ( index_type i )
    const // this const means that this method do not
          // modify the contents of the struct
    { return this -> values[i] ; }

  } vec ;

  typedef struct {
    value_type values[dim][dim];
    
    // declare the operator ()
    value_type & operator () ( index_type i, index_type j) 
    { return this -> values[i][j] ; }
    
    value_type const & operator () ( index_type i, index_type j)
      const // this const means that this method do not
            // modify the contents of the struct
    { return this -> values[i][j] ; }
  } mat ;
 
  /*
   now we set the definition (prototipe)
   of the functions
   */
    
  // prototype of the non-linear system 
  extern vec F( vec const & in ) ;

  // prototype of the Jacobian of the non-linear system
  extern mat JF( vec const & in ) ;
  
  // prototype for the solver
    
  typedef vec F_type  ( vec const & ) ;
  typedef mat JF_type ( vec const & ) ;

  extern
  bool
  solver( vec        const & x0,
          F_type           * pF,
          JF_type          * pJF,
          index_type const   n,
          value_type const & tol,
          vec              & x
         ) ;
}

#endif

20.1.cpp

File: 20.1.cpp

#include "20.h"

namespace newton_solver {
  
  /*
   Small 2x2 linear solver
   */
  
  bool
  linear_solver( mat const & A, // coefficients matrix
                 vec const & b, // r.h.s
                 vec       & x  // the solution of the linear system
               ) {
    // Using Cramer rule
    value_type det = A(0,0) * A(1,1) - A(1,0) * A(0,1) ;
    if ( det == 0 ) return false ;
    
    x[0] = ( b[0] * A(1,1) - b[1] * A(0,1) ) / det ;
    x[1] = ( A(0,0) * b[1] - A(1,0) * b[0] ) / det ;
    
    // b[0] is the same as:
    // b . operator [] (0)
    return true ;
  }

  /*
     Newton scheme
                         -1
     x    = x   - JF(x  )   F(x )
      n+1    n        n        n
   
     the implementation:
    
       SOLVE: JF(x ) h = F(x )
                  n         n

       UPDATE:  x   = x + h
                n+1    n
   */
  bool
  solver( vec        const & x0,
          F_type           * pF,
          JF_type          * pJF,
          index_type const   n,
          value_type const & tol,
          vec              & x
          ) {
    
    // set x to the initial value
    x = x0 ; // not reccomended
    
    for ( index_type i = 0 ; i < n ; ++i ) {
      vec f, h ;
      mat jf ;
      // perform a Newton step 
      // substep 1: compute F(x)
      f = pF(x) ;
      // substep 2: compute JF(x)
      jf = pJF(x) ;
      // substep 3: compute the correction h = JF(x)^(-1)F(x)
      bool ok = linear_solver(jf, f, h) ;
      if ( !ok ) return false ;
      // substep 4: update x 
      // x = x - h ;
      x -= h ;
      // check for the stopping criteria
      value_type res = h . norm2() ;
      if ( res <= tol ) return true ; 
    }
    return false ;
  }
}

20.2.cpp

File: 20.2.cpp

#include "20.h"

namespace newton_solver {
  // F(x) = / x0 - x1           \
  //        |                   |
  //        \ x0**2 + x1**3 - 2 /
  vec
  F( vec const & x ) {
    vec Fx ;
    Fx[0] = x[0] - x[1] ; // in . values[0] - in . values[1]
    Fx[1] = x[0]*x[0] - x[1]*x[1]*x[1] - 2 ;
    return Fx ;
  }

  //
  mat
  JF( vec const & in ) {
    mat JFx ;
    value_type const * x = in . values ;
    JFx(0,0) = 1            ; // DF[0]/Dx[0]
    JFx(0,1) = -1           ; // DF[0]/Dx[1]
    JFx(1,0) = 2*x[0]       ; // DF[1]/Dx[0]
    JFx(1,1) = -3*x[1]*x[1] ; // DF[1]/Dx[1]
    return JFx ;
  }
}

int
main() {
  using newton_solver::vec ;
  using newton_solver::F ;
  using newton_solver::JF ;
  using newton_solver::solver ;
  using namespace std ;

  vec x0, x ;
  
  x0.values[0] = 1 ;
  x0.values[1] = 1 ;
  solver( x0, F, JF, 100, 1E-6, x ) ; 
  
  cout << "The solution is x = [ " 
       << x.values[0] << ", "
       << x.values[1] << "]\n" ;
}

20.h

File: 20.h

// standard trick for only one time inclusion
#ifndef TWENTIETH_H
#define TWENTIETH_H

#include <iostream>
#include <iomanip>
#include <cmath>

namespace newton_solver {
  using namespace std ;
  
  // set the value_type to double,
  // value_type is the generic name for "number"
  // (non integer)type
  typedef double value_type ;

  // set the index_type to int,
  // index_type is the generic name for "integer" type
  typedef int index_type ;
  
  index_type const dim = 2 ;
  
  typedef struct vec {
    value_type values[dim] ;
    
    // declare the operator []
    value_type & operator [] ( index_type i ) 
    { return this -> values[i] ; }
    
    value_type const & operator [] ( index_type i )
    const // this const means that this method do not
          // modify the contents of the struct
    { return this -> values[i] ; }
    
    vec const &
    operator -= ( vec const & a ) {
      for ( index_type i = 0 ; i < dim ; ++i ) 
        values[i] -= a[i] ;
      return * this ;
    }

    vec const &
    operator += ( vec const & a ) {
      for ( index_type i = 0 ; i < dim ; ++i ) 
        values[i] += a[i] ;
      return * this ;
    }
    
    // compute the 2 norm of the vector
    value_type norm2() const {
      value_type res = 0 ;
      for ( index_type i = 0 ; i < dim ; ++i ) 
        res += values[i] * values[i] ;
      return sqrt(res) ;
    }

  } vec ;
  
  // external operator
  
  // add two vector
  inline
  vec
  operator + ( vec const & a, vec const & b ) {
    vec c ;
    for ( index_type i = 0 ; i < dim ; ++i )
      c[i] = a[i] + b[i] ;
  }

  // subtract two vector (definition only)
  vec
  operator - ( vec const & a, vec const & b ) ;
  
  typedef struct {
    value_type values[dim][dim];
    
    // declare the operator ()
    value_type & operator () ( index_type i, index_type j) 
    { return this -> values[i][j] ; }
    
    value_type const & operator () ( index_type i, index_type j)
      const // this const means that this method do not
            // modify the contents of the struct
    { return this -> values[i][j] ; }
  } mat ;
 
  /*
   now we set the definition (prototipe)
   of the functions
   */
    
  // prototype of the non-linear system 
  extern vec F( vec const & in ) ;

  // prototype of the Jacobian of the non-linear system
  extern mat JF( vec const & in ) ;
  
  // prototype for the solver
    
  typedef vec F_type  ( vec const & ) ;
  typedef mat JF_type ( vec const & ) ;

  extern
  bool
  solver( vec        const & x0,
          F_type           * pF,
          JF_type          * pJF,
          index_type const   n,
          value_type const & tol,
          vec              & x
         ) ;
}

#endif

21.1.cpp

File: 21.1.cpp

#include "21.1.h"

namespace vector_matrix {

  value_type &  // return value
  vec::operator []   // the operator [] is defined here
  ( index_type i ) // the argument of []
  { return this -> values[i] ; }
  
  value_type const &
  vec::operator [] ( index_type i ) const
  { return this -> values[i] ; }
    
  vec const &
  vec::operator -= ( vec const & a ) {
    for ( index_type i = 0 ; i < dim ; ++i ) 
      values[i] -= a[i] ;
    return * this ;
  }
  
  vec const &
  vec::operator += ( vec const & a ) {
    for ( index_type i = 0 ; i < dim ; ++i ) 
      values[i] += a[i] ;
    return * this ;
  }
    
  // compute the 2 norm of the vector
  value_type vec::norm2() const {
    value_type res = 0 ;
    for ( index_type i = 0 ; i < dim ; ++i ) 
      res += values[i] * values[i] ;
    return sqrt(res) ;
  }

  // external operator
  
  // add two vector
  vec
  operator + ( vec const & a, vec const & b ) {
    vec c ;
    for ( index_type i = 0 ; i < dim ; ++i )
      c[i] = a[i] + b[i] ;
    return c ;
  }

  // subtract two vector (definition only)
  vec
  operator - ( vec const & a, vec const & b ) {
    vec c ;
    for ( index_type i = 0 ; i < dim ; ++i )
      c[i] = a[i] - b[i] ;
    return c ;
  }
  
  // declare the operator ()
  value_type &
  mat::operator () ( index_type i, index_type j) 
  { return this -> values[i][j] ; }
    
  value_type const &
  mat::operator () ( index_type i, index_type j) const
  { return this -> values[i][j] ; }
  
  
  /*
   Small 2x2 linear solver
   */
  
  vec
  operator / ( vec const & b, mat const & A) {
    vec x ;
    // Using Cramer rule
    value_type det = A(0,0) * A(1,1) - A(1,0) * A(0,1) ;
    
    x[0] = ( b[0] * A(1,1) - b[1] * A(0,1) ) / det ;
    x[1] = ( A(0,0) * b[1] - A(1,0) * b[0] ) / det ;
    
    return x ;
  }

}

21.2.cpp

File: 21.2.cpp

#include "21.2.h"

namespace newton_solver {
  /*
     Newton scheme
                         -1
     x    = x   - JF(x  )   F(x )
      n+1    n        n        n
   
     the implementation:
    
       SOLVE: JF(x ) h = F(x )
                  n         n

       UPDATE:  x   = x + h
                n+1    n
   */
  bool
  solver( vec        const & x0,
          F_type           * pF,
          JF_type          * pJF,
          index_type const   n,
          value_type const & tol,
          vec              & x
          ) {
    
    // set x to the initial value
    x = x0 ; // not reccomended
    
    for ( index_type i = 0 ; i < n ; ++i ) {
      vec f, h ;
      mat jf ;
      // perform a Newton step 
      // substep 1: compute F(x)
      f = pF(x) ;
      // substep 2: compute JF(x)
      jf = pJF(x) ;
      // substep 3: compute the correction h = JF(x)^(-1)F(x)
      h = f / jf ;
      //bool ok = linear_solver(jf, f, h) ;
      //if ( !ok ) return false ;
      // substep 4: update x 
      // x = x - h ;
      x -= h ;
      // check for the stopping criteria
      value_type res = h . norm2() ;
      if ( res <= tol ) return true ; 
    }
    return false ;
  }
}

21.3.cpp

File: 21.3.cpp

#include "21.2.h"

namespace newton_solver {
  // F(x) = / x0 - x1           \
  //        |                   |
  //        \ x0**2 + x1**3 - 2 /
  vec
  F( vec const & x ) {
    vec Fx ;
    Fx[0] = x[0] - x[1] ; // in . values[0] - in . values[1]
    Fx[1] = x[0]*x[0] - x[1]*x[1]*x[1] - 2 ;
    return Fx ;
  }

  //
  mat
  JF( vec const & x ) {
    mat JFx ;
    JFx(0,0) = 1            ; // DF[0]/Dx[0]
    JFx(0,1) = -1           ; // DF[0]/Dx[1]
    JFx(1,0) = 2*x[0]       ; // DF[1]/Dx[0]
    JFx(1,1) = -3*x[1]*x[1] ; // DF[1]/Dx[1]
    return JFx ;
  }
}

int
main() {
  using newton_solver::vec ;
  using newton_solver::F ;
  using newton_solver::JF ;
  using newton_solver::solver ;
  using namespace std ;

  vec x0, x ;
  
  x0[0] = 1 ;
  x0[1] = 1 ;
  solver( x0, F, JF, 100, 1E-6, x ) ; 
  
  cout << "The solution is x = [ " 
       << x[0] << ", "
       << x[1] << "]\n" ;
}

21.1.h

File: 21.1.h

// standard trick for only one time inclusion
#ifndef TWENTIONE_ONE_H
#define TWENTIONE_ONE_H

#include <iostream>
#include <iomanip>
#include <cmath>

namespace vector_matrix {
  using namespace std ;

  // set the value_type to double,
  // value_type is the generic name for "number"
  // (non integer)type
  typedef double value_type ;

  // set the index_type to int,
  // index_type is the generic name for "integer" type
  typedef int index_type ;
  
  index_type const dim = 2 ;
  
  class mat ; // forward declaration
  
  class vec {
    // private member of the class
    value_type values[dim] ;
    
  public: // from now the public interface
    // declare the operator []
    value_type & operator [] ( index_type i ) ;    
    value_type const & operator [] ( index_type i ) const ;
    
    vec const & operator -= ( vec const & a );
    vec const & operator += ( vec const & a );    

    // compute the 2 norm of the vector
    value_type norm2() const ;

  } ;
  
  // external operator
  
  // add two vector
  vec operator + ( vec const & a, vec const & b ) ;

  // subtract two vector (definition only)
  vec operator - ( vec const & a, vec const & b ) ;
  
  vec operator / ( vec const & b, mat const & A ) ;
  
  class mat {
    value_type values[dim][dim];
    
  public:
    // declare the operator ()
    value_type & operator () ( index_type i, index_type j) ;
    value_type const & operator () ( index_type i, index_type j) const ;
  } ;


}

#endif

21.2.h

File: 21.2.h

// standard trick for only one time inclusion
#ifndef TWENTIONE_TWO_H
#define TWENTIONE_TWO_H

#include <iostream>
#include <iomanip>
#include <cmath>

#include "21.1.h"

namespace newton_solver {
  using namespace std ;
  using namespace vector_matrix ;
 
  /*
   now we set the definition (prototipe)
   of the functions
   */
    
  // prototype of the non-linear system 
  extern vec F( vec const & in ) ;

  // prototype of the Jacobian of the non-linear system
  extern mat JF( vec const & in ) ;
  
  // prototype for the solver
    
  typedef vec F_type  ( vec const & ) ;
  typedef mat JF_type ( vec const & ) ;

  extern
  bool
  solver( vec        const & x0,
          F_type           * pF,
          JF_type          * pJF,
          index_type const   n,
          value_type const & tol,
          vec              & x
         ) ;
}

#endif

22.1.cpp

File: 22.1.cpp

#include "22.1.h"

namespace vector_matrix {

  // the constructor
  vec::vec( index_type dim ) {
    this -> dim = dim ;
    // new allocate dynamicaly the memory
    // new TYPE [ SIZE ]
    values = new value_type [ dim ] ;
  }

  // the destructor
  vec::~vec() {
    delete [] values ;
  }
  
  vec const & vec::operator = ( vec const & a ) {
    for ( index_type i = 0 ; i < dim ; ++i )
      values[i] = a[i] ;
    return *this ;
  }

  value_type &  // return value
  vec::operator []   // the operator [] is defined here
  ( index_type i ) // the argument of []
  { return this -> values[i] ; }
  
  value_type const &
  vec::operator [] ( index_type i ) const
  { return this -> values[i] ; }
    
  vec const &
  vec::operator -= ( vec const & a ) {
    for ( index_type i = 0 ; i < dim ; ++i ) 
      values[i] -= a[i] ;
    return * this ;
  }
  
  vec const &
  vec::operator += ( vec const & a ) {
    for ( index_type i = 0 ; i < dim ; ++i ) 
      values[i] += a[i] ;
    return * this ;
  }
    
  // compute the 2 norm of the vector
  value_type vec::norm2() const {
    value_type res = 0 ;
    for ( index_type i = 0 ; i < dim ; ++i ) 
      res += values[i] * values[i] ;
    return sqrt(res) ;
  }

  // external operator
  
  // add two vector
  vec
  operator + ( vec const & a, vec const & b ) {
    index_type dim = a.get_dim() ;
    vec c(dim)  ;
    for ( index_type i = 0 ; i < dim ; ++i )
      c[i] = a[i] + b[i] ;
    return c ;
  }

  // subtract two vector (definition only)
  vec
  operator - ( vec const & a, vec const & b ) {
    index_type dim = a.get_dim() ;
    vec c(dim)  ;
    for ( index_type i = 0 ; i < dim ; ++i )
      c[i] = a[i] - b[i] ;
    return c ;
  }
  
  // the constructor
  mat::mat( index_type dim ) {
    this -> dim = dim ;
    // new allocate dynamicaly the memory
    // new TYPE [ SIZE ]
    values = new value_type [ dim * dim ] ;
  }
  
  // the destructor
  mat::~mat() {
    delete [] values ;
  }
  
  // the copy contructor
  mat const & mat::operator = ( mat const & a ) {
    for ( index_type i = 0 ; i < dim*dim ; ++i )
      values[i] = a . values[i] ;
    return *this ;
  }
  
  // declare the operator ()
  value_type &
  mat::operator () ( index_type i, index_type j) 
  { return this -> values[i+j*dim] ; }
    
  value_type const &
  mat::operator () ( index_type i, index_type j) const
  { return this -> values[i+j*dim] ; }
  
  
  /*
   Small 2x2 linear solver
   */
  
  vec
  operator / ( vec const & b, mat const & A) {
    vec x(2) ;
    // Using Cramer rule
    value_type det = A(0,0) * A(1,1) - A(1,0) * A(0,1) ;
    
    x[0] = ( b[0] * A(1,1) - b[1] * A(0,1) ) / det ;
    x[1] = ( A(0,0) * b[1] - A(1,0) * b[0] ) / det ;
    
    return x ;
  }

}

22.2.cpp

File: 22.2.cpp

#include "22.2.h"

namespace newton_solver {
  /*
     Newton scheme
                         -1
     x    = x   - JF(x  )   F(x )
      n+1    n        n        n
   
     the implementation:
    
       SOLVE: JF(x ) h = F(x )
                  n         n

       UPDATE:  x   = x + h
                n+1    n
   */
  bool
  solver( vec        const & x0,
          F_type           * pF,
          JF_type          * pJF,
          index_type const   n,
          value_type const & tol,
          vec              & x
          ) {
    
    index_type dim = x0 . get_dim();

    // set x to the initial value
    x = x0 ; // not reccomended
    
    for ( index_type i = 0 ; i < n ; ++i ) {
      vec f(dim), h(dim) ;
      mat jf(dim) ;
      // perform a Newton step 
      // substep 1: compute F(x)
      f = pF(x) ;
      // substep 2: compute JF(x)
      jf = pJF(x) ;
      // substep 3: compute the correction h = JF(x)^(-1)F(x)
      h = f / jf ;
      //bool ok = linear_solver(jf, f, h) ;
      //if ( !ok ) return false ;
      // substep 4: update x 
      // x = x - h ;
      x -= h ;
      // check for the stopping criteria
      value_type res = h . norm2() ;
      if ( res <= tol ) return true ; 
    }
    return false ;
  }
}

22.3.cpp

File: 22.3.cpp

#include "22.2.h"

namespace newton_solver {
  // F(x) = / x0 - x1           \
  //        |                   |
  //        \ x0**2 + x1**3 - 2 /
  vec
  F( vec const & x ) {
    vec Fx(2) ;
    Fx[0] = x[0] - x[1] ; // in . values[0] - in . values[1]
    Fx[1] = x[0]*x[0] - x[1]*x[1]*x[1] - 2 ;
    return Fx ;
  }

  //
  mat
  JF( vec const & x ) {
    mat JFx(2) ;
    JFx(0,0) = 1            ; // DF[0]/Dx[0]
    JFx(0,1) = -1           ; // DF[0]/Dx[1]
    JFx(1,0) = 2*x[0]       ; // DF[1]/Dx[0]
    JFx(1,1) = -3*x[1]*x[1] ; // DF[1]/Dx[1]
    return JFx ;
  }
}

int
main() {
  using newton_solver::vec ;
  using newton_solver::F ;
  using newton_solver::JF ;
  using newton_solver::solver ;
  using namespace std ;

  vec x0(2), x(2) ;
  
  x0[0] = 1 ;
  x0[1] = 1 ;
  solver( x0, F, JF, 100, 1E-6, x ) ; 
  
  cout << "The solution is x = [ " 
       << x[0] << ", "
       << x[1] << "]\n" ;
}

22.1.h

File: 22.1.h

// standard trick for only one time inclusion
#ifndef TWENTITWO_ONE_H
#define TWENTITWO_ONE_H

#include <iostream>
#include <iomanip>
#include <cmath>

namespace vector_matrix {
  using namespace std ;

  // set the value_type to double,
  // value_type is the generic name for "number"
  // (non integer)type
  typedef double value_type ;

  // set the index_type to int,
  // index_type is the generic name for "integer" type
  typedef int index_type ;
  
  class mat ; // forward declaration
  
  class vec {
    
    index_type dim ;
    
    // private member of the class
    value_type * values ;
    
    vec() ;
        
  public: // from now the public interface
      
    // the constructor
    vec( index_type dim ) ;

    // the destructor
    ~vec() ;

    index_type const & get_dim() const { return dim ; }

    // the copy constructor
    vec const & operator = (vec const & ) ;

    // declare the operator []
    value_type & operator [] ( index_type i ) ;    
    value_type const & operator [] ( index_type i ) const ;
    
    vec const & operator -= ( vec const & a );
    vec const & operator += ( vec const & a );    

    // compute the 2 norm of the vector
    value_type norm2() const ;

  } ;
  
  // external operator
  
  // add two vector
  vec operator + ( vec const & a, vec const & b ) ;

  // subtract two vector (definition only)
  vec operator - ( vec const & a, vec const & b ) ;
  
  vec operator / ( vec const & b, mat const & A ) ;
  
  class mat {
    
    index_type dim ;
    
    // private member of the class
    value_type * values ;
    
    mat() ;

  public: // from now the public interface
      
    // the constructor
    mat( index_type dim ) ;
    
    // the destructor
    ~mat() ;

    // the copy constructor
    mat const & operator = (mat const & ) ;

    // declare the operator ()
    value_type & operator () ( index_type i, index_type j) ;
    value_type const & operator () ( index_type i, index_type j) const ;
  } ;


}

#endif

22.2.h

File: 22.2.h

// standard trick for only one time inclusion
#ifndef TWENTITWO_TWO_H
#define TWENTITWO_TWO_H

#include <iostream>
#include <iomanip>
#include <cmath>

#include "22.1.h"

namespace newton_solver {
  using namespace std ;
  using namespace vector_matrix ;
 
  /*
   now we set the definition (prototipe)
   of the functions
   */
    
  // prototype of the non-linear system 
  extern vec F( vec const & in ) ;

  // prototype of the Jacobian of the non-linear system
  extern mat JF( vec const & in ) ;
  
  // prototype for the solver
    
  typedef vec F_type  ( vec const & ) ;
  typedef mat JF_type ( vec const & ) ;

  extern
  bool
  solver( vec        const & x0,
          F_type           * pF,
          JF_type          * pJF,
          index_type const   n,
          value_type const & tol,
          vec              & x
         ) ;
}

#endif

Lesson of 25/5/2005

24.cpp

File: 24.cpp

#include <iostream>

/*
  Overloading and templates 
 
 */

/*
  compute the maximum of two numbers 
 */

/* function max for double precision numbers */
// first definition

double max( double const & a, double const & b ) {
  if ( a > b ) return a ;
  else         return b ;
}

/* function max for single precision numbers */
// second definition

float max( float const & a, float const & b ) {
  if ( a > b ) return a ;
  else         return b ;
}

/* overloading means the same name for differents functions */

/*
  The template mechanism can produce skeleton
  of code that can be expanded when is needed
 */

template <typename T>
T max( T const & a, T const & b ) {
  if ( a > b ) return a ;
  else         return b ;
}


template <typename Ta, typename Tb>
Ta max( Ta const & a, Tb const & b ) {
  if ( a > b ) return a ;
  else         return b ;
}


int
main() {
  {
    double a = 1 ;
    double b = 2 ;
    double c = max(a,b) ;// use first definition
  }
  {
    float a = 1 ;
    float b = 2 ;
    float c = max(a,b) ;// use second definition
  }
  {
    int a = 1 ;
    int b = 2 ;
    int c = max(a,b) ; // use the template with T=int
  }
  {
    int  a = 1 ;
    long b = 2 ;
    int  c = max(a,b) ; // use the template with T=int
  }
}

Lesson of 30/5/2005

25.cpp

File: 25.cpp

#include <iostream>

using namespace std ;
/*
  Overloading and templates 
 */

/*
  The template mechanism can produce skeleton
  of code that can be expanded when is needed
 */

template <typename T, int N>
T sum( T const v[] ) {
  T res = 0 ;
  for ( int i = 0 ; i < N ; ++i ) 
    res += v[i] ;
  return res ;
}

extern double b[] ; // extern double b[] ;
                    // somewhere there is double b[12] ;

int
main() {
  double a[] = {1,2,3,4,5,6,7,8,9} ;
  int const n = sizeof(a)/sizeof(a[0]); //sizeof(a)/sizeof(double) ; 
  cout << n << endl ;
  double res = sum<double,5>(a) ;
  cout << res << endl ;
  res = sum<double,n>(a) ;
  cout << res << endl ;
  res = sum<double,n+10000>(a) ;
  cout << res << endl ;
}

26.cpp

File: 26.cpp

#include <iostream>
#include <string>

using namespace std ;

/*
 A class can be thinked as a structure
 with all the methods and data public.
 */

// this class
class A {
  public:
    int a, b ;
} ;

// is "nearly" equivalent to this struct

struct {
  int a, b ;
} A ;

/*
 One of the important feature of a class is
 that we can protect some information
 */

class B {
  string name ;
  double x, y ;
public:
  // constructor method
  B(string const & _n, 
    double const & _x,
    double const & _y) 
  : name(_n), x(_x), y(_y)
  { cout << "passing to the FIRST constructor for " << name << "\n" ; }
  // another constructor
  B(string const & _n) : name(_n), x(0), y(0)
  { cout << "passing to the SECOND constructor for " << name << "\n" ; }
  // another constructor
  B(string const & _n, int const _x, int const _y) 
  : name(_n), x(_x), y(_y)
  { cout << "passing to the THIRD constructor for " << name << "\n" ; }
  
  // the destructor
  ~B() {
    cout << "passing to the destructor for " << name << "\n" ; 
  }
} ;

B a("a",0,1) ; // call of the third constuctor

int
main() {
  cout << "STARTING THE PROGRAM\n" ;
  B b("b",0.0,1.0)  ; // call of the first constuctor

  B * pointer_to_d ;
  pointer_to_d = new B("pointer",0.0,1.0) ;
  cout << "INSIDE THE PROGRAM\n" ;
  delete pointer_to_d ; // free memory for pointer_to_d
  
  B c("c")          ; // call of the second constuctor
  
  cout << "ENDING THE PROGRAM\n" ;
  return 0 ; 
}

27.cpp

File: 27.cpp

#include <iostream>
#include <string>

using namespace std ;

/*
 Example of a class with many costructor
 */

template <typename T>
class vector {
public:
  typedef T value_type ; // alias of T to value_type ;
private:
  string       the_name ;
  int          the_size ;
  value_type * the_data ;
public:
  // constructor method
  vector(string const & n, int const size) : the_name(n), the_size(size) {
    cout << "Building vector " << the_name << " size = " << the_size << "\n" ;
    the_data = new value_type [size] ;
  }
  
  // the destructor
  ~vector() {
    cout << "passing to the destructor for " << the_name << "\n" ; 
    delete [] the_data ;
  }
  
  // define the [] operator
  value_type const & operator [] ( int i ) {
    if ( i < 0 || i >= the_size ) {
      cerr << "bad index for " << the_name << "[" << i << "]\n" ;
      exit(0) ; // exit from the program
    }
    return the_data[i] ;
  }

} ;

typedef vector<double> VECTOR ;

int
main() {
  cout << "STARTING THE PROGRAM\n" ;
  VECTOR v("v",100) ;
  
  for ( int i = 0 ; i < 200 ; ++i ) {
    VECTOR::value_type a ;
    a = v[i] ;
  }
  cout << "ENDING THE PROGRAM\n" ;
  return 0 ; 
}

28.cpp

File: 28.cpp

#include <iostream>
#include <string>
#include <sstream>

using namespace std ;

/*
 Example of a class with many costructor
 */

template <typename T>
class vector {
public:
  typedef T value_type ; // alias of T to value_type ;
private:
  string       the_name ;
  int          the_size ;
  value_type * the_data ;
public:
  // constructor method
  vector(string const & n, int const size) : the_name(n), the_size(size) {
    cout << "Building vector " << the_name << " size = " << the_size << "\n" ;
    the_data = new value_type [size] ;
  }
  
  // the destructor
  ~vector() {
    cout << "passing to the destructor for " << the_name << "\n" ; 
    delete [] the_data ;
  }
  
  // define the [] operator
  value_type const & operator [] ( int i ) {
    if ( i < 0 || i >= the_size ) {
      // throw 1.234 ;
      ostringstream the_error ;
      the_error << "bad index for " << the_name << "[" << i << "]\n" ;
      throw the_error.str(); // throw a string object
    }
    return the_data[i] ;
  }

} ;

typedef vector<double> VECTOR ;

int
main() {
  cout << "STARTING THE PROGRAM\n" ;
  VECTOR v("v",100) ;
  
  try {
    for ( int i = 0 ; i < 200 ; ++i ) {
      VECTOR::value_type a ;
      a = v[i] ;
    }
  }
  catch (int err_code) {
    cerr << "found error number " << err_code << "\n" ; 
  }
  catch (string const & err ) {
    cerr << "found error : " << err << "\n" ;     
  }
  catch (...) {
    cerr << " unknown error found\n" ;
  }
  cout << "ENDING THE PROGRAM\n" ;
  return 0 ; 
}

29.cpp

File: 29.cpp

#include <iostream>
#include <string>
#include <sstream>

using namespace std ;

/*
  Example of an error class
 */

class error {
  string the_error_message ;  
  int    the_error_number ;
public:
  error(string const & s, int n) :
    the_error_message(s),
    the_error_number(n) {}

  void
  print( ostream & stream ) const {
    stream << "ERROR: " << the_error_message 
           << "\nERROR number: " << the_error_number
           << "\n" ;
  }
};

/*
 Example of a class with many costructor
 */

template <typename T>
class vector {
public:
  typedef T value_type ; // alias of T to value_type ;
private:
  string       the_name ;
  int          the_size ;
  value_type * the_data ;
public:
  // constructor method
  vector(string const & n, int const size) : the_name(n), the_size(size) {
    cout << "Building vector " << the_name << " size = " << the_size << "\n" ;
    the_data = new value_type [size] ;
  }
  
  // the destructor
  ~vector() {
    cout << "passing to the destructor for " << the_name << "\n" ; 
    delete [] the_data ;
  }
  
  // define the [] operator
  value_type const & operator [] ( int i ) {
    if ( i < 0 || i >= the_size ) {
      ostringstream the_error ;
      the_error << "bad index for " << the_name << "[" << i << "]\n" ;
      throw error(the_error.str(),1) ; // throw a string object
    }
    return the_data[i] ;
  }

} ;

typedef vector<double> VECTOR ;

int
main() {
  cout << "STARTING THE PROGRAM\n" ;
  VECTOR v("v",100) ;
  
  try {
    for ( int i = 0 ; i < 200 ; ++i ) {
      VECTOR::value_type a ;
      a = v[i] ;
    }
  }
  catch (int err_code) {
    cerr << "found error number " << err_code << "\n" ; 
  }
  catch (string const & err ) {
    cerr << "found error : " << err << "\n" ;     
  }
  catch (error const & err ) {
    err . print(cerr) ;     
  }
  catch (...) {
    cerr << " unknown error found\n" ;
  }
  cout << "ENDING THE PROGRAM\n" ;
  return 0 ; 
}

30.cpp

File: 30.cpp

#include <iostream>
#include <string>

using namespace std ;

/*
  Example of inheritance
 */

// define the base class

class Base {
  string the_name ;
  double x, y ;
public:
  Base( string const & n ) : the_name(n), x(0), y(0) {}
  Base( string const & n,
        double const & _x,
        double const & _y) : the_name(n), x(_x), y(_y) {}
  
  string const & name() const { return the_name ; }
  
  void print( ostream & stream ) const {
     stream << "class: " << the_name << " " << x << " " << y << "\n" ; 
  }
} ;

class Derived : public Base {
  public:
  Derived(string const & n) : Base(n,1,1) {}
} ;

int
main() {
  cout << "STARTING THE PROGRAM\n" ;
  Base    a("a") ;
  Derived b("b") ;
  a . print(cout) ;
  b . print(cout) ;
  cout << "ENDING THE PROGRAM\n" ;
  return 0 ; 
}

31.cpp

File: 31.cpp

#include <iostream>
#include <string>

using namespace std ;

/*
  A classical Example of inheritance
 */

// define the base class

class Object {
  string the_name ;
protected:
  double cx, cy ; // center of the object
  double area ;   // area of the object
public:
  Object( string const & n ) : the_name(n) {}

  string const & name() const { return the_name ; }
  
  void print( ostream & stream ) const {
     stream << "object: " << the_name 
            << "\ncenter " << cx << " " << cy
            << "\narea   " << area << "\n" ; 
  }
} ;

class Circle : public Object {
  double ray ;
  // here cx, cy, area are in the private part of Circle
  double const PI ; 
public:
  Circle(double const & _cx,
         double const & _cy,
         double const & _ray) :
    Object("circle"),
    ray(_ray),
    PI(3.1415)
  {
    // You can use the scope :: operator to access
    // inherited part of the class
    Object::cx = _cx ;
    Object::cy = _cy ;
    Object::area = 0 ;
      
    // If there is no possibility of confusion
    // you can avoid the scope operator ::
    cx = _cx ;
    cy = _cy ;
    area = ray*ray*PI ;
  }

  void print( ostream & stream ) const {
    Object::print(stream) ;
    stream << "perimeter: " << 2*ray*PI << "\n" ; 
  }
  
} ;

class Square : public Object {
  double length ;
  // here cx, cy, area are in the private part of Circle
public:
    Square(double const & _cx,
           double const & _cy,
           double const & _length) :
    Object("square"),
    length(_length)
  {
    // If there is no possibility of confusion
    // you can avoid the scope operator ::
    cx   = _cx ;
    cy   = _cy ;
    area = length*length ;
  }
} ;

void
print( Object const & o ) {
  o . print ( cout ) ; 
}

int
main() {
  cout << "STARTING THE PROGRAM\n" ;
  Circle c(2,3,3) ;
  Square q(2,3,3) ;
  c . print(cout) ;
  print(c) ; // loose the print of the perimeter
  print(q) ; 
  cout << "ENDING THE PROGRAM\n" ;
  return 0 ; 
}

32.cpp

File: 32.cpp

#include <iostream>
#include <string>

using namespace std ;

/*
  A classical Example of inheritance
 */

// define the base class

class Object {
  string the_name ;
protected:
  double cx, cy ; // center of the object
  double area ;   // area of the object
public:
  Object( string const & n ) : the_name(n) {}

  string const & name() const { return the_name ; }
  
  virtual // declare the routine as virtual 
  void print( ostream & stream ) const {
     stream << "object: " << the_name 
            << "\ncenter " << cx << " " << cy
            << "\narea   " << area << "\n" ; 
  }
} ;

class Circle : public Object {
  double ray ;
  // here cx, cy, area are in the private part of Circle
  double const PI ; 
public:
  Circle(double const & _cx,
         double const & _cy,
         double const & _ray) :
    Object("circle"),
    ray(_ray),
    PI(3.1415)
  {
    // You can use the scope :: operator to access
    // inherited part of the class
    Object::cx = _cx ;
    Object::cy = _cy ;
    Object::area = 0 ;
      
    // If there is no possibility of confusion
    // you can avoid the scope operator ::
    cx = _cx ;
    cy = _cy ;
    area = ray*ray*PI ;
  }

  void print( ostream & stream ) const {
    Object::print(stream) ;
    stream << "perimeter: " << 2*ray*PI << "\n" ; 
  }
  
} ;

class Square : public Object {
  double length ;
  // here cx, cy, area are in the private part of Circle
public:
    Square(double const & _cx,
           double const & _cy,
           double const & _length) :
    Object("square"),
    length(_length)
  {
    // If there is no possibility of confusion
    // you can avoid the scope operator ::
    cx   = _cx ;
    cy   = _cy ;
    area = length*length ;
  }
} ;

void
print( Object const & o ) {
  o . print ( cout ) ; 
}

int
main() {
  cout << "STARTING THE PROGRAM\n" ;
  Circle c(2,3,3) ;
  Square q(2,3,3) ;
  c . print(cout) ;
  print(c) ; // NOW do not loose the print of the perimeter
  print(q) ; 
  cout << "ENDING THE PROGRAM\n" ;
  return 0 ; 
}

Lesson of 30/5/2005 (pomeriggio)

33.cpp

File: 33.cpp

#include <iostream>

using namespace std ;

/*
  A classical Example of generic programming
 */

// a simple routine that sort a vector of number
void
sort( int v[], int const dim ) {
  // sorting by bubble sort
  for ( int i = 0 ; i < dim-1 ; ++i ) 
    for ( int j = i+1 ; j < dim ; ++j ) 
      if ( v[i] > v[j] ) {
        // swap the values
        int k = v[i] ;
        v[i] = v[j] ;
        v[j] = k ;
      }
}

/*
  Using template we can generalize the algorithm
  for any data type, other than integer
 */

template <typename T>
void
sort( T v[], int const dim ) {
  // sorting by bubble sort
  for ( int i = 0 ; i < dim-1 ; ++i ) 
    for ( int j = i+1 ; j < dim ; ++j ) 
      if ( v[i] > v[j] ) {
        // swap the values
        T k = v[i] ;
        v[i] = v[j] ;
        v[j] = k ;
      }
}


/*
  Sorting algorithm in the spirit of 
  the Standard Template Library
 
 */

template <typename T>
void
sort( T * begin, // pointer to the first element of the vector
      T * end) { // pointer to PAST to the last element of the vector
  // sorting by bubble sort
  for ( T * pi = begin ; pi < end-1 ; ++pi ) 
    for ( T * pj = pi+1 ; pj < end ; ++pj ) 
      if ( *pi > *pj ) {
        // swap the values
        T buff = *pi ;
        *pi = *pj ;
        *pj = buff ;
      }
}

int
main() {
  //int v[] = {1,3,-3,2,23,123,-23,232,-3,-45,-2} ;
  //int dim = sizeof(v)/sizeof(v[0]) ;
  long v[] = {1,3,-3,2,23,123,-23,232,-3,-45,-2} ;
  long dim = sizeof(v)/sizeof(v[0]) ;

  sort(v,v+dim) ; // v       = pointer to the first element
                  // v + dim = pointer to past to the last element
  cout << v[0] ;
  for ( int i = 1 ; i < dim ; ++i ) 
    cout << ", " << v[i] ;
  cout << "\n" ;
  return 0 ; 
}

34.cpp

File: 34.cpp

#include <iostream>

// Standard template library for algorithm
#include <algorithm>

using namespace std ;

int
main() {
  //int v[] = {1,3,-3,2,23,123,-23,232,-3,-45,-2} ;
  //int dim = sizeof(v)/sizeof(v[0]) ;
  
  long a[] = {1,3,-3,2,23,123,-23,232,-3,-45,-2} ;
  long adim = sizeof(a)/sizeof(a[0]) ;

  long b[] = {22,-34,12} ;
  long bdim = sizeof(b)/sizeof(b[0]) ;

  sort(a,a+adim) ; // a       = pointer to the first element
                   // a + dim = pointer to past to the last element

  sort(b,b+bdim) ; // b       = pointer to the first element
                   // b + dim = pointer to past to the last element

  for ( int i = 0 ; i < adim ; ++i ) 
    if ( a[i] == 23 ) 
      cout << "FOUND at rank " << i << endl ;
  
  /*
     lower_bound implements the binary search
     in a sorted vector ( o vector like ).
     If found return the pointer (iterator) of the found
     thing. If not found return the pointer (iterator)
     of past to the last element.
  */
  long * res = lower_bound( a, a + adim, 1111123) ;
  if ( res == a + adim || *res != 1111123 ) 
    cout << "NOT FOUND\n" ;
  else
    cout << " FOUND = " << *res
         << " RANK = " << res - a
         << endl ;

  cout << a[0] ;
  for ( int i = 1 ; i < adim ; ++i ) 
    cout << ", " << a[i] ;
  cout << "\n" ;
  
  /*
    merge: merge two sorted vector into a 
           big one sorted vector
   */
  long c[1000] ;
  merge( a, a+adim,
         b, b+bdim,
         c ) ;
  
  cout << c[0] ;
  for ( int i = 1 ; i < adim+bdim ; ++i ) 
    cout << ", " << c[i] ;
  cout << "\n" ;
  
  cout << "MAX a = " << * max_element( a, a + adim) << endl ;
  
  c[0] = 1 ;
  c[1] = 2 ;
  c[2] = 3 ;
  for ( int i = 0 ; i < 7 ; ++i ) {
    cout << c[0] << " " << c[1] << " " << c[2] << endl ;
    next_permutation( c, c+3 ) ;
  }
  
  return 0 ; 
}

35.cpp

File: 35.cpp

#include <iostream>

// Standard template library for algorithm
#include <algorithm>

// Standard template library for vector container
#include <vector>

using namespace std ;

int
main() {
  long avec[] = {1,3,-3,2,23,123,-23,232,-3,-45,-2} ;
  long adim = sizeof(avec)/sizeof(avec[0]) ;

  vector<long> a ;

  a . reserve( 10 ) ; // reserve memory
  
  for ( int i = 0 ; i < adim ; ++i ) {
    a . push_back(avec[i]) ;
    cout << "insert " << avec[i]
         << " size = " << a . size()
         << " capacity = " << a . capacity() << endl ;
  }

  sort( a . begin(), a . end() ) ;
  
  vector<long> b(3) ; // 3 is the size of the vector
  b[0] = 22 ;
  b[1] = -34 ;
  b[2] = 12 ;

  sort( b . begin(), b . end() ) ;

  /*
     lower_bound implements the binary search
     in a sorted vector ( o vector like ).
     If found return the pointer (iterator) of the found
     thing. If not found return the pointer (iterator)
     of past to the last element.
  */
  
  vector<long>::iterator
    res = lower_bound( a . begin(),
                       a . end(),
                       23 ) ;
  
  if ( res == a . end() || *res != 23 ) 
    cout << "NOT FOUND\n" ;
  else
    cout << " FOUND = " << *res
         << " RANK = " << res - a . begin()
         << endl ;

  cout << a . front() ; // equivalent *(a.begin())
  for ( vector<long>::const_iterator i = a.begin()+1 ;
        i < a . end() ; ++i ) 
    cout << ", " << *i ;
  cout << "\n" ;

  /*
    merge: merge two sorted vector into a 
           big one sorted vector
   */
  
  vector<long> c(a.size()+b.size()) ;
  merge( a . begin(), a . end(),
         b . begin(), b . end(),
         c.begin() ) ;
  
  cout << c . front() ;
  for ( vector<long>::const_iterator i = c.begin()+1 ;
        i < c . end() ; ++i ) 
    cout << ", " << *i ;
  cout << "\n" ;
  
  cout << "MAX a = " << * max_element( a.begin(), a.end()) << endl ;

  c . resize(3) ;
  c[0] = 1 ;
  c[1] = 2 ;
  c[2] = 3 ;
  for ( int i = 0 ; i < 7 ; ++i ) {
    cout << c[0] << " " << c[1] << " " << c[2] << endl ;
    next_permutation( c . begin(), c . end() ) ;
  }

  return 0 ; 
}

---