/*******************************************************************************
Treibergs                                                           Mar. 1, 2006   

Solution to HW 5a.   the sigma function sigma(n)

sigma(n)  is the sum of the proper divisors of n.

We print a table of the first 200 sigmas. It is evident that
6  and  28  are the only perfect numbers  < 100.  496  is also perfect.

Then it asks one number at a time and determines sigma, which is coded as a 
function. This part of the program is all that was required for your HW.


sigma.c
*******************************************************************************/

# include <stdlib.h>
# include <stdio.h>

int
sigma( int n );

int
main( void )
{
    int i, j, m=0, n;

    printf("sigma(n) = The sum of the proper divisors of  n.\n\n             ");

    for (i=0; i<= 9; i++)
            printf ( "%7d", i);
    for( i = 0; i <= 190 ; i=i+10 )
    {
    printf("\n            +----------------------------------------------------------------------\n%10d  |",i);

    
          for( j = 0; j<= 9; j++)
          {     
                n = i+j;
               if ( n > 0 )
                {
                    m = sigma(n);
                    printf("%7d", m);
                    if (m == n)
                      {
                       printf("\a <-------- HEY!  %d  IS PERFECT!\a\n            |", n);
                       for(m=0; m<= j; m++)
                                printf("       ");
                       }
                }
                else
                    printf("       ");

           }
    }
        printf("\n            +----------------------------------------------------------------------\n");

     m = 1;
     do
    {
          printf( "\n Enter positive   n  (Nonpositive stops the program.) : ");
          scanf("%d", &i);

          if (i <= 0)
                  exit(0);
          else
                  { j = sigma(i);
                   printf("  The sum of proper divisors  sigma(%d) = %d", i, j );
                   if (i == j)
                       printf("\a <-------- HEY!  %d  IS PERFECT!\a", i);
                   printf("\n\n");
                   }
    }
     while( (m++ <= 10) && (i > 0) );


    return EXIT_SUCCESS;
}

int
sigma( int n )
{
    int i,k=0;
         for(i = 1; 2*i <= n; i++)
             if( n % i == 0 )
                   k = k + i;
     
     return k;
}

/***********************************************************************************
Treibergs                                                              Mar. 14, 2006 

HW 5c. Use Simpson's Rule to integrate (a + x^3)^(1/2) between 1 and 3
with error <= 0.5e-6.. 

Input a. Integrate given function using Simpson's rule. Double 
the number of intervals each step. Note that at each step, we only need
to evaluate the function at the points midway between the points from the 
previous step. Use the fact that Simpson's rule is (trap + 2 midpt)/3

For Simpson's rule with  n  intervals, the error is  E = (b-a)^5 M / ( 2880 n^4 ) 
where  M  = sup{ |f''''(c)| : c real }.  

Now I'll estimate  M. You can do this by plotting f'''' to see what values it takes. 
I'll compute f'''' and group the terms so that the estimate is obvious.

In our case  f' = (3/2) x^2 (a  +  x^3)^(-1/2)
            f'' = 3x (a  +  x^3/4) (a  +  x^3)^(-3/2)
           f''' = 3 (a^2 - (5/2) a x^3 - (1/8) x^6) (a + x^3)^(-5/2)
          f'''' = 9 x^2 (-5 a^2  +  (7/2) a x^3  +  (1/16) x^6) (a + x^3)^(-7/2).
In the range  1 <= x <= 3,  how big does |f''''(x)| get? The key observation is
Q = -5 a^2  +  (7/2) a x^3  +  (1/16) x^6 <= (7/4) (a  +  x^3)^2  (Multiply it out!)
The other thing is that  Q's  cubes of zeros are
       (7a +/- sqrt(49 + 5a^2) )  /  20    <=    ( 7 + sqrt(54) ) / 20 <= .8 < 1
which says that 0 < Q  so  f''''>0  for 1 <= x.
Hence, we can partition the 
(a + x^3)^(-7/2) = (a + x^3)^(-5/6) (a + x^3)^(-2/3)(a + x^3)^(-2)
so
|f''''|<= 9(a + x^3)^(-5/6) x^2(a + x^3)^(-2/3) *
                *(-5 a^2  +  (7/2) a x^3  +  (1/16) x^6) (a + x^3)^(-2)
       <= 9  *  1  *  (7/4)  <  8.
 Thus  M = 8 suffices for  1 <= a <= b.
(If you were to graph  f'''', you'll see thar |f''''| <= 1.2 over  1 <= x <= 3.)
      
Finally, we just use our Simpson's program with the correct function,
and stop when the error < 0.5e-6.

M2160s5c.c
*************************************************************************************/
      
# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# define  ERRORTOL  (0.5e-6)

double aaa;

double
f ( double x );

int 
main ( void ) 
{
    double a, ahalf, b, e, error=1.0e9, h, s, simps, t, M= 8.0;
    int i = 1, m = 1, n = 1;
    
    
    printf("Enter the value for  a  : ");
    scanf ( "%lf", &aaa );
    
    printf ( "\n\n\n Simpson's Rule\n\n" );
    
    a = 1.0;
    b = 3.0;
   
    printf ( " The left endpoint a = %f.  The right endpoint  b = %f.\n", a, b );
    
   
    printf("\n\t   n \t  Simpson's Rule\t\t\t  Error in Simpson's Rule\n");
    
    h = b  -  a;
    e = h * h * h * h * h * M / 2880.0;
    t = ( f(a)  +  f(b) ) / 2.0;
    while ( (error > ERRORTOL ) && (  m++ <= 15) )
    {
       ahalf = a  +  h / 2.0 ;
       s = 0.0;
       
       for( i = 0; i < n; i++ )
                s = s + f ( ahalf + i*h );
       
       error =  e/((double)n*n*n*n);
       simps = (t  +  2.0 * s) * h / 3.0;
       printf ("%11d %22.15f %30.15e\n", n, simps, error);
       t = t  +  s;
       h = h / 2.0;
       n = 2.0 * n;
     }
     printf ( "\n\nSimpson's Rule Integral = %f\n", simps );
    
     return EXIT_SUCCESS;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

double
f ( double x )
{
    return sqrt( aaa  +  x*x*x );
}

/*****************************************************************************
Treibergs                                                        Feb. 24, 2006 

Homework Problem  5b: (To see what a solution may look like in C++, see
extra program below!) 


Find all zeros of

f(x) = x^5  +  2x^4  -  3x^3  -  4x^2  +  4x   +  aaa
f'(x) = 5x^4  +  8x^8  -  9 x^3  -  8x  +  4;

where  0 < aaa < 1.

A fifth degree polynomial has at most five real roots. The roots fall
between the two tails and are spliced between the two relative minima and maxima.
If one of the minima is positive, as is is for this polynomial, then there
are two fewer zeros. 

We tabulate  f  to see that there seem to be three zeros and a possible
positive relative minimum. We tabulate  f'  to see where the critical points
are. By doing bisection on  f'  we see that there is a critical point at
x = 1,  f  is decreasing for  x<1  and increasing for  x>1  and that  f(1)>0
so that  x=1  is a positive minimum. Hence there are at most three real zeros
of the polynomial. We use the tabulation of  f  to bracket the zeros.

In doing this project I expect you to run the bisections as separate runs, 
and just tell me what the conclusions are. There may be other arguments
to show that you are not missing any roots. I have programmed the several 
steps into one program. Run the following program. It does all of these things
in sequence to finds the zeros of the polynomial. To find the zeros, it calls 
the bisection subroutine, which runs until it finds an approximation within 
0.5e-6  of the zero.

There is one techicality in the program about using functions as arguments.
We wish to run bisection on both the derivatve  df  and the function  f.
The idea is to call  bisect(a,b,f)  but C doesn't understand this. Instead
you have to pass a pointer to the function. The prototype says bisect is a
the double valued function

double
bisect (double a, double b, double (*)(double) );

and that one of the arguments is a pointer to a double valued function
of a double argument.Inside the function subprogram, the variable function 
is called by
                 else if ( (*g)(a) * (*g)(m) <  0.0 )
where  g  is a pointer,  *g  is the value at the pointed address.
Note that extra parentheses are required for correct parsing by  C. This
code basically says to use the value returned by the function pointed to 
by  g.  Finally, inside the calling program, we use
    m  =  bisect(.86, 1.12, (double(*)(double))(df) );
which tells  C  that   df  is a pointer to a double valued function of a 
double. (This is tricky. See section 5.11 of Kernighan & Ritchie.)


M2160s5b.c
*****************************************************************************/
      
# include <stdio.h>
# include <stdlib.h>
# include <math.h>

double aaa;     /* global variable accessible by all program units. */

double
f ( double x );

double
df ( double x );


double
bisect (double , double , double (*)(double) ); /* function passed as arg */

int 
main ( void ) 
{
    double   m, x1, x2, x3 ;
    int i;
    
    printf("Zeros of the given polynomial using the bisection blgorithm\n\n");
    printf(" Enter the constant  a : ");
    scanf( "%lf", &aaa );
    
    printf("\n Table of Values of the function  f(x)\n   x\t\t\t  f(x)\n");
    for(i=-10; i<= 10; i++)
    printf("%12.6f\t%21.15f\n",-.18+0.26*i,f(-.18+0.26*i));
    
    printf("\nThere seems to be a positive local minimum at about  x=1"
           " and three sign changes.\n\n");
    printf("Use bisection to find the valus of zero near x=1 of  f'(x).\n");
           
    printf("\n Table of Values of its derivative  f'(x)\n   x\t\t\t  f'(x)\n");
    for(i=-10; i<= 10; i++)
    printf("%12.6f\t%21.15f\n",-.18+0.26*i,df(-.18+0.26*i));
    
    printf("\nThere seems to be four zeros of f' at about  x=-2,-1,-0.4,1.\n\n");
    printf("Using the  f'  table, we run bisection on  f'  to find the\n"
           "critical point between %f < x < %f\n", .86, 1.12 );

    m  =  bisect(.86, 1.12, (double(*)(double))(df) );
    
    printf("\nf' is negative to the left and positive to the right, "
    "\nhence there is a minimum of  f  at  x = %f\n",m);
    printf("But  f(%f) = %f  so this critical point is a positive\n",m,f(m));
    printf("local minimum of  f  so there is no zero of  f  near  x = %f\n",m);
    
    printf("\nNow use the  f  table to find brackets for the three remaining\n"
           "zeros. Then run bisection on the function  f  near those zeros.\n");
    
    x1 = bisect(-2.26,-2.0, (double(*)(double))(f) );
    x2 = bisect(-2.0, -1.48, (double(*)(double))(f) );
    x3 = bisect(-0.44, 0.08, (double(*)(double))(f) );
      
    printf("\n\nBecause one of the relative minima of  f  is positive, there\n"
      "are only three possible remaining roots of the 5th degree polynomial.\n");
   
    printf("  x  =  %f  is the first root of the polynomial.\n",x1);
    printf("  x  =  %f  is the second root of the polynomial.\n",x2);
    printf("  x  =  %f  is the third root of the polynomial.\n",x3);

    return EXIT_SUCCESS;
}

double
bisect( double a, double b, double(*g)(double))
{
     double error, m;
     int i;
               error = fabs(b - a);
     printf("\n n\t\tLeft\t\tMiddle\t\tRight\t\t\tError\n");
              for(i=1; i <= 25;  i++)
            {
                 m = (a+b)/2.0;
                 error = error / 2.0;
                 printf("%3d%20.15f%20.15f%20.15f%15.6e\n",i,a,m,b,error);
                 if( error < 0.5e-6 )
                 {
                    printf(" **===** The function is zero at %20.15g \n", m);
                       break;
                 }
                 else if ( (*g)(a) * (*g)(m) <  0.0 )
                       b = m;
                 else
                       a = m;       
             }
    return m;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

double
f ( double x )
{
   return x*(x*(x*(x*(x + 2.0) -3.0)-4.0)  + 4.0 )  +  aaa;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/* the derivative of the given function */

double
df ( double x )
{
   return x*(x*(x*(5.0* x + 8.0) -  9.0)  -  8.0)  + 4.0   ;
}




/*********************************************************************************
**********************************************************************************
**********************************************************************************
**********************************************************************************
The following code, written in the C++ language, has been generously provided
by professor Kai-Uwe Bux of the University of Virginia. Despite my protestations,
he volunteered to solve an assigned problem to demonstrate the power of C++.
It may give you some idea of what extensions C++ offers beyond ANSI C.  C++ 
supports a wider variety of data types, which enable one to think about object in
more of a mathematical, rather than programming way. A.T.
**********************************************************************************
**********************************************************************************
**********************************************************************************
*********************************************************************************/
Hi Andrejs,


here is my C++ solution to problem 5/B. Feel free to use
it in any way you like. You will probably find that the
program is of little to no use to your students. However, I
think, it illustrates nicely some of the points I made about
how one can write C++ in a way that parallels the naive way
of thinking about the problem.

A few points to note:

a) I define (ad hoc) a polynomial type that supports evaluation:
   If P is a polynomial, I can write P(5.0) and it will just
   evaluate.

b) This is not just syntactic sugar: in one place I evaluate
   a polynomial on all points in a vector. Since polynomials
   support evaluation, I can do that using std::transform().

c) I freely used iterator magic. You may want to ignore those
   lines in a first reading.

-----------
Kai-Uwe Bux
e-mail:    bux_2002@kubux.net
home page: http://www.kubux.net



// polynomial_solver.cc (C) Kai-Uwe Bux [2006]
// ===========================================

#include <vector>
#include <utility>
#include <cmath>
#include <iterator>
#include <algorithm>


// generalities on real numbers

bool have_same_sign ( double a, double b ) {
  if ( a >= 0 && b >= 0 ) {
    return ( true );
  }
  if ( a <= 0 && b <= 0 ) {
    return ( true );
  }
  return ( false );
}


// polynomials

struct Polynomial : public std::vector< double > {

  // FIXME: [bad hack]
  /*
    | This is a quick hack. A good implementation of polynomials
    | should support many more operations:
    | - adding, subtracting, and multiplying polynomials
    | - long division of polynomials
    | - composition of polynomial functions
    | - maybe polynomials in several variables
    | Also, this should be a template; and it probably should not
    | inherit from std::vector!
    |
    | In real life, I would just use the polynomial template
    | from my library :-)
  */

  int degree ( void ) const {
    int result = this->size() - 1;
    while ( ( result >= 0 ) && (*this)[result] == 0 ) {
      -- result;
    }
    return ( result );
  }

  double operator() ( double arg ) const {
    double result = 0;
    for ( const_reverse_iterator iter = this->rbegin();
	  iter != this->rend(); ++iter ) {
      result *= arg;
      result += *iter;
    }
    return ( result );
  }

  double leading_coefficient ( void ) const {
    int index = this->degree();
    if ( index >= 0 ) {
      return ( (*this)[index] );
    } else {
      return ( 0.0 );
    }
  }

}; 


// taking derivatives:

Polynomial diff ( Polynomial const & P ) {
  Polynomial result;
  for ( Polynomial::size_type i = 1; i < P.size(); ++i ) {
    result.push_back( P[i] * i );
  }
  return ( result );
}


// asymptotic behavior of polynomials

bool right_limit_positive ( Polynomial const & P ) {
  return ( P.leading_coefficient() > 0 );
}

bool left_limit_positive ( Polynomial const & P ) {
  return ( ( P.leading_coefficient() < 0 ) == ( P.degree() % 2 != 0 ) );
}


// finding points left and right that have the sign of the limit:

double far_right_point ( Polynomial const & P, double from ) {
  while ( ( right_limit_positive(P) &&  P( from ) <= 0 )
	  ||
	  ( ( ! right_limit_positive(P) ) && P( from ) >= 0 ) ) {
    from += 1;
  }
  return ( from );
}

double far_left_point ( Polynomial const & P, double from ) {
  while ( ( left_limit_positive(P) && P( from ) <= 0 )
	  ||
	  ( ( ! left_limit_positive(P) ) &&P( from ) >= 0 ) ) {
    from -= 1;
  }
  return ( from );
}


// the bisection method in one interval

double bisect ( Polynomial const & P, double left, double right ) {
  double left_value = P(left);
  double right_value = P(right);
  double middle = ( left + right ) / 2;
  // WARNING: [tricky loop condition]
  /*
    | Don't try:
    |
    |   while ( left == right ) {
    |     ...
    |   }
    |
    | as that could loop forever.
  */
  while ( ( left < middle ) && ( middle < right ) ) {
    assert( have_same_sign( left_value, -right_value ) );
    double middle_value = P( middle );
    if ( have_same_sign( left_value, -middle_value ) ) {
      right_value = middle_value;
      right = middle;
    } else {
      left_value = middle_value;
      left = middle;
    }
    middle = ( left + right ) / 2;
  }
  if ( std::abs( left_value ) < std::abs( right_value ) ) {
    return ( left );
  } else {
    return ( right );
  }
}


// find all zeros of a real polynomial

typedef std::vector< double > DoubleVector;

DoubleVector zeros ( Polynomial const & P ) {
  DoubleVector result;
  if ( P.degree() <= 0 ) {
    return ( result );
  }
  if ( P.degree() == 1 ) {
    result.push_back( - P[0]/P[1] );
  } else {
    // Key idea: Between two consecutive points with horizontal tangent,
    // there is at most one zero; thus, we can run the bisection method
    // in those intervals.
    DoubleVector horizontal = zeros( diff( P ) );
    // WARNING: [rounding errors]
    /*
      | Because of rounding errors, the call above might not
      | exactly yield the points with horizontal tangents. Thus, we
      | cannot really guarantee that we will find all zeros. A careful
      | analysis of the numercial stabilty seems difficult.
    */
    DoubleVector breaks;
    // We need to take care of possibly one more zero to the left
    // an one more zero to the right.
    if ( horizontal.empty() ) {
      breaks.push_back( far_left_point( P, -1.0 ) );
      breaks.push_back( far_right_point( P, 1.0 ) );
    } else {
      breaks.push_back( far_left_point( P, horizontal.front() ) );
      std::copy( horizontal.begin(), horizontal.end(),
		 std::back_inserter( breaks ) );
      breaks.push_back( far_right_point( P, horizontal.back() ) );
    }
    DoubleVector values;
    std::transform( breaks.begin(), breaks.end(),
		    std::back_inserter( values ), P );
    for ( DoubleVector::size_type n = 0; n + 1 < breaks.size(); ++n ) {
      if ( have_same_sign( values[n], -values[n+1] ) ) {
	result.push_back( bisect( P, breaks[n], breaks[n+1] ) );
      }
    }
  }
  return ( result );
}


#include <iostream>

int main ( void ) {
  Polynomial P;
  P.push_back( 0.123 );
  P.push_back( 4 );
  P.push_back( -4 );
  P.push_back( -3 );
  P.push_back( 2 );
  P.push_back( 1 );
  DoubleVector output = zeros( P );
  std::cout << "I found the following zeros:\n";
  std::copy( output.begin(), output.end(),
	     std::ostream_iterator< double >( std::cout, "\n" ) );
}