/*******************************************************************************
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
# include
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.5e6..
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 = (ba)^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.5e6.
M2160s5c.c
*************************************************************************************/
# include
# include
# include
# define ERRORTOL (0.5e6)
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.5e6 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
# include
# include
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.5e6 )
{
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 KaiUwe 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.

KaiUwe Bux
email: bux_2002@kubux.net
home page: http://www.kubux.net
// polynomial_solver.cc (C) KaiUwe Bux [2006]
// ===========================================
#include
#include
#include
#include
#include
// 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
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" ) );
}