Lesson 1.16: Looping statements (Part 1: while loops)
Up to now, we have seen a while loop. There may be certain cases where it is necessary to execute the body of the loop at least once; for example, where the execution of the continued loop depends on an action taken during the execution of some condition set up during the exeuction.
The easiest example is the execution of the body of a loop based on a user response:
Suppose we are querying the user as to what action should be performed:
unsigned int n; // Uninitialized
std::cout << "Enter a number to test if it is prime (enter '0' to quit): ";
std::cin >> n;
while ( n != 0 ) {
// Test if 'n' is prime...
std::cout << "Enter a number to test if it is prime (enter '0' to quit): ";
std::cin >> n;
}
It seems unnecessary to have to have code that asks twice. A do-while loop guarantees that the loop will execute at least once, and the test is performed at the end of the exeuction of the body (or, if you prefer, at the start of the execution of the next loop of the body):
unsigned int n; // Uninitialized
do {
std::cout << "Enter a number to test if it is prime (enter '0' to quit): ";
std::cin >> n;
if ( n != 0 ) {
// Test if 'n' is prime...
}
} while ( n != 0 );
Now, the body of the loop must be executed at least once.
High-low
In playing the game high-low, we ahve a similar situation: If we were to use a while loop, we would have to initially ask the user for a guess, and then if the user was wrong, we would loop:
int target{ (std::rand() % 100) + 1 };
unsigned int guess; // Uninitialized
std::cout << "Enter a guess (1-100): ";
std::cin >> guess;
if ( (guess == 0) || (guess > 100) ) {
std::cout << "Your guess must be between 1 and 100: ";
std::cin >> guess;
}
while ( guess != target ) {
if ( guess < target ) {
std::cout << "Low" << std::endl;
} else if ( guess > target ) {
std::cout << "High" << std::endl;
}
std::cout << "Enter a guess (1-100): ";
std::cin >> guess;
if ( (guess == 0) || (guess > 100) ) {
std::cout << "Your guess must be between 1 and 100: ";
std::cin >> guess;
}
}
It seems unnecessary to have to have code that asks twice. A do-while loop guarantees that the loop will execute at least once, and the test is performed at the end of the exeuction of the body (or, if you prefer, at the start of the execution of the next loop of the body):
will be executed once if the condition is true, and it will never execute if the condition is false.
As a more serious example, in calculating $n!$, we start with $2$, multiply that by $3$ to get $6$, multiply that by $4$ to get $24$, and so on until we finally multiply the last result by $n$.
To implement such an operation, we need some way of performing an operation multiple times, and for that we introduce a while statement:
while ( some-condition ) {
// Body of the while loop
}
If the condition is never true, the body of the loop will never be executed. If the condition is always true, the body of the loop will be executed repeatedly forever. The flow chart for such a statement is:
We continue doing something until the condition is false.
Now, don't run this next example as we have a potential issue: suppose that the condition is always true. In this case, the loop will continue executing forever:
#include <iostream>
int main();
int main() {
int n{0};
while ( n == 0 ) {
std::cout << "Hello world!" << std::endl;
}
return 0;
}
This loop will never exit—it will continue printing Hello world! forever or until you end the execution (Ctrl-c), the computer malfunctions or the next blackout occurs.
A loop that never terminates is called an infinite loop. An infinite loop is likely a problem if you're trying to solve a problem: there are two possible sources of an infinite loop:
- an error with your program, or
- a problem with no solution.
Almost all of your infinite loops for the next two years will be of the first kind; however, when you get to second year, you will be finding approximations to mathematical problems, and the likelihood of the second will increase.
There is, however, one kind of infinite loop that is very common, and actually required: an infinite loop in an embedded system. Consider your Tickle-Me-ElmoTM where a microcontroller controls the behavior. A microcontroller is a processor usually less powerful than a general-purpose microprocessor found in your laptop or desktop computers, but where some peripherals of laptops and desktops are built into the microcontroller. When you turn on your Tickle-Me-Elmo, the microcontroller runs in a loop waiting for the child to do something to which it should respond.
However, for the purpose of this course, an infinite loop will always be considered a semantic programming error—if you get into an infinite loop, you did something wrong.
Collatz conjecutre (or the $3n + 1$ problem)
The $3n + 1$ problem is defined as follows: Start with any positive integer $a_0 = n$. We will now define a sequence of integers $a_1, a_2, ...$ where the value of $a_n$ depends on $a_{n - 1}$ where
- if $a_{n - 1}$ is odd, $a_n = 3a_{n - 1} + 1$, and
- if $a_{n - 1}$ is even, $a_n = a_{n - 1}/2$.
We keep going until $a_n = 1$, for after this, it repeats: $1, 4, 2, 1, 4, 2, 1, ...$.
Some examples are:
- $3, 10, 5, 16, 8, 4, 2, 1$ and
- $7, 22, 11, 34, 17, 52, 26, 13, 40, 20, 10, 5, 16, 8, 4, 2, 1$
It has been conjectured (the Collatz conjecture) that all such sequences ultimately end at $1$, but this has never been proven: there may be some sequence that, over time, goes to $\infty$.
Write a function collatz( int n ) that returns the value $n$ such that $a_n = 1$. We will say that $n$ is the Collatz number of $a_0$.
Here is algorithm we will follow:
- Set $n \leftarrow 0$.
- If $a \ne 1$,
- Increment $n$ (set $n \leftarrow n + 1$).
- If $n$ is even, divide $n$ by two,
otherwise, set $n$ to the value $3n + 1$. - Go back to Step 2.
Otherwise, return $n$.
This can be shown in the flow chart
Writing this as a function, we have:
unsigned int collatz( unsigned int a0 );
unsigned int collatz( unsigned int a0 ) {
unsgined int counter{0};
while ( a0 != 1 ) {
++counter;
if ( (a0 % 2) == 0 ) {
a0 /= 2; // identical to 'a0 = a0/2;'
} else {
a0 = 3*a0 + 1;
}
}
return counter;
}
For the counter variable, we use counter as opposed to n so that the reader gets a better idea as to what the local variable is for. The parameter is labelled as a0 just to indicate that what the user is passing is the initial value in the sequence.
Now, let us print the first 100 Collatz numbers:
- Set an integer $n \leftarrow 1$.
- If $n \le 100$,
- Call collatz($n$) and print the statement "collatz(n) = that value.
- Increment $n$ (set $n \leftarrow n + 1$).
- Go back to Step 2.
A flow chart for this program would be:
Writing this as a function main(), we have:
#include <iostream>
int main();
unsigned int collatz( unsigned int n );
int main() {
unsigned int n{1};
while ( n <= 100 ) {
std::cout << "collatz(" << n << ") = " << collatz( n ) << std::endl;
++n;
}
return 0;
}
// Definition of collatz()...
Now, we already know that if $a_0 = 2^n$, then $a_n = 1$. However, if $a_0 = 27$, then we must wait until $a_{111} = 1$. Let us only print $n$ if the Collatz number is larger than any previous Collatz number. If you were asked to do this, you would start calculating ${\rm collatz}(n)$ for ever increasing values of $n$, but you would track the largest value you had found so far. Then, you would only announce that a new Collatz number if it is greater than the largest value you've already found; after which, you would update that largest value:
- Set $n \leftarrow 1$.
- Set $max \leftarrow 0$.
- If $n \le 10000$,
- If ${\rm collatz}(n) > max$,
- Update $max$ with that new value, and
- print the statement "collatz(n) = that value.
- Increment $n$ (set $n \leftarrow n + 1$).
- Go back to Step 3.$
- If ${\rm collatz}(n) > max$,
- Return 0.
A flow chart for this program would be:
Writing this as a function main(), we have:
#include <iostream>
int main();
unsigned int collatz( unsigned int n );
int main() {
unsigned int n{1};
unsigned int largest_found{0};
while ( n <= 10000 ) {
if ( collatz( n ) > largest_found ) {
largest_found = collatz( n );
std::cout << "collatz(" << n << ") = "
<< largest_found << std::endl;
}
++n;
}
return 0;
}
// Definition of collatz()...
Now, you may have noticed a small issue with the previous: each time we update largest_found, we calculate collatz( n ) three times. We can calculate it once, store its value in a local variable, and then use this local variable if necessary:
#include <iostream>
int main();
unsigned int collatz( unsigned int n );
int main() {
unsigned int n{1};
unsigned int largest_found{0};
while ( n <= 10000 ) {
unsigned int cn{collatz( n )};
if ( cn > largest_found ) {
largest_found = cn;
std::cout << "collatz(" << n << ") = "
<< largest_found << std::endl;
}
++n;
}
return 0;
}
// Definition of collatz()...
Determining if a number is prime
Our first serious example of an algorithm that you learned in secondary school is determining if an integer $n$ is prime. This requires you to try to divide the integer by all integers up to $\sqrt{n}$. We could visualize this as a flow chart:
It is important to understand what the problem you are trying to solve is first before you code it. This is why flow charts are useful: you can give this flow chart to a non-programmer, and that person could reasonably be expected to execute this algorithm.
Note, $\sqrt{n}$ may not be an integer, so instead of checking if $k \le \sqrt{n}$, we will check if $k^2 \le n$.
Long division
Another algorithm you saw in secondary school (and likely in primary school) is long division. While you may have learned the procedure, it's quite unlikely that you could easily write down how someone else could implement this algorithm. Here is one means of performing the operation of long division: Suppose we are calculating $m \div n$. First, we must assume that each number has a finite number of decimal places.
- First, if both $m$ and $n$ have the same sign, the result will be positive; otherwise, the result will be negative. Thus, our solution begins with either "+" or "-". Now discard the signs of $m$ and $n$.
- Second, multiply both $m$ and $n$ by the smallest power of $10$ that allows $n$ to be an integer. For example, instead of calculating $39.527 \div 1.73$, we are now calculating $3952.7 \div 173$. Similarly, instead of calculating $35.25 \div 1300$, we are now calculating $0.3525 \div 13$. Update both $m$ and $n$ to be these two new numbers.
- Third,
- if $|m| < 1$, append a "0." to the solution and let $r$ equal the first digit after the decimal place.
- Otherwise, let $r$ equal the first non-zero digit.
- Now, iterate:
- Find the largest multiple $d$ of $n$ such that $r - dn \ge 0$, append $d$ to the solution, and let $r - dn$ be the new $r$.
- Add the next digit of $m$ to the end of the new $r$ and if the next digit is the first digit after the decimal place, append a decimal point to the end of the solution.
Factorial
The factorial $n!$ is defined as the product of all numbers between $2$ and $n$.
Binomial coefficients
The binomial coefficients $n \choose k$ can be found to be equal to $\frac{n!}{k! (n - k)!}$, which requires the calculation of three separate factorials.
A faster algorithm in finding $n \choose k$ is to first let $k$ be the minimum of $k$ and $n - k$, and to then calculate $\frac{n(n - 1)\cdots(n - k - 1)}{k!}$, which requires a total of only $2k - 2$ multiplications and not $2n - 6$ multiplications.
Greatest common denominator
In calculating the greatest common demoninator (gcd) of two positive integers, we have the following flow chart:
Factorial function
How would you calculate $n!$? First, you would start with a result $r = 1$, and then you would start multiplying the result by $2$, then $3$, and every integer up to $n$. You could describe your algorithm as follows:
- Set a result to $1$,
- Set a multiplier to $2$,
- If the multiplier is less than or equal to $n$,
- Multiply the result by the multiplier, and
- add one to the multiplier.
- Go back to Step 3.
When this algorithm finishes, the result will be the value of $n!$. In a flow chart, this algorithm may be shown graphically as follows:
To program this into a computer, we can start with:
unsigned int result{1};
unsigned int multiplier{2};
This initializes the result and the multiplier. Next, we have a loop
while ( multiplier <= n ) {
// Multiply the result by 'multiplier' and increment 'multiplier'
}
Thus, we have the loop
while ( multiplier <= n ) {
// Multiply the result by 'multiplier' and increment 'multiplier'
result *= multiplier; // identical to 'result = result*multiplier;'
++multiplier; // identical to 'multiplier = multiplier + 1;'
}
Thus, our factorial function will be:
unsigned int factorial( unsigned int n );
unsigned int factorial( unsigned int n ) {
unsigned int result{1};
unsigned int multiplier{2};
while ( multiplier <= n ) {
result *= multiplier;
++multiplier;
}
return result;
}
By observation, every time the body of the loop executes, the value of multiplier must be incremented by 1, so this loop should execute a finite number of times, so this will not lead to an infinite loop.
Greatest common divisor
You will recall that the greatest common divisor of two positive integers $m$ and $n$ is that positive integer $d$ such that $d$ divides both $m$ and $n$ and no number larger than $d$ divides both $m$ and $n$.
How can we find the greatest common divisor of $m$ and $n$; that is, how do we find ${\rm gcd}(m,n)$?
If $m \ge n$ and $n$ divides $m$ (that is, the remainder of $m \div n$ is zero or $m {\rm mod} n = 0$), then the gcd of these two is $n$. If, however, one does not divide the other, we have a problem: how do we find the greatest common factor?
Suppose that $d$ is the greatest common divisor of both $m$ and $n$. In that case, if $m > n$, we can both write $m = d\tilde{m}$, $n = d\tilde{n}$ but not only that but $m = qn + r$ where $n < r < 0$. Combining these, we have $d\tilde{m} = qd\tilde{n} + r$, so $d(\tilde{m} - \tilde{n}) = r$, so the remainder, too, must be divisible by the greatest common divisor $d$.
Thus, ${\rm gcd}(m, n) = {\rm gcd}(n, r)$. Now, we repeat the algorithm.
To give an example of this, the ${\rm gcd}(27720, 5265) = 45$. We note that $27720 = 616 \times 45$ and $5265 = 117 \times 45$ where ${\rm gcd}(616, 117) = 1$.
- Now, $27720 = 5\times 5265 + 1395$, and we see that $1395 = 31 \times 45$.
- Next, $5265 = 3\times 1395 + 1080$, and we see that $1080 = 24 \times 45$.
- Next, $1395 = 1\times 1080 + 315$, and we see that $315 = 7 \times 45$.
- Next, $1080 = 3\times 315 + 135$, and we see that $135 = 3 \times 45$.
- Next, $315 = 2\times 135 + 45$, and we see that $45 = 1 \times 45$.
- Finally, $135 = 3\times 45 + 0$, so we are done.
Thus, our program (or algorithm) is expressible by:
An outline of the function is:
unsigned int gcd( unsigned int m, unsigned int n {
// While 'n' does not divide 'm'
// - Determine the remainder of 'm' divided by 'n'
// - Let 'm' be assigned the value of 'n'
// - Let 'n' be assigned the remainder that was found
// Return the second
}
Now, the phrase if $n$ does not divide $n$ is the same as saying that $m {\rm mod} n \ne 0$, or in C++, (m % n) != 0. To find the remainder, we can assign the result of this operation to a local variable: r = m % n. Thus, our algorithm is:
unsigned int gcd( unsigned int m, unsigned int n );
unsigned int gcd( unsigned int m, unsigned int n ) {
while ( (m % n) != 0 ) {
unsigned int r{m % n};
m = n;
n = r;
}
return n;
}
Finding the integer square root of an integer
The integer square root of an integer $n$ is that integer $s$ that is less than or equal to $\sqrt{n}$. Thus, the following two conditions must be satisfied:
- The integer square root $s$ must be such that $s^2 \le n$, and
- because it is the largest integer that satisfies this condition, it must be true that $(s + 1)^2 > n$.
Thus, this condition is true only when we have found the integer square root. While this condition is not true, we must get a better approximation. As you can read on Wikipedia, the integer square root can be found by starting with any positive integer $s_0$ and from this calculating $s_1 = \left \lfloor \frac{1}{2}\left( s_0 + \frac{n}{s_0} \right ) \right \rfloor$, and continuing this, in general
$s_{k + 1} = \left \lfloor \frac{1}{2}\left( s_k + \frac{n}{s_k} \right ) \right \rfloor$
where $\lfloor x \rfloor$ is the floor function, also described as the grin function or greatest integer less than or equal to $x$. For example, $\lfloor -3.2 \rfloor = -4$, $\lfloor -3.0 \rfloor = -3$, $\lfloor 3.2 \rfloor = 3$ and $\lfloor 4.0 \rfloor = 4$.
Now, why does this formula give us a better approximation of the square root? First, we see that if $s_0 > \sqrt{n}$ then $\frac{1}{s_0} < \frac{1}{\sqrt{n}}$, so multiplying both sides by $n > 0$ gives us that $\frac{n}{s_0} < \frac{n}{\sqrt{n}} = \sqrt{n}$. Thus, $\frac{n}{s_0} < \sqrt{n} < s_0$, and therefore the average of these two ($s_1$) must always be closer to $\sqrt{n}$ than the original value of $s_0$. Additionally, can see that $\sqrt{n} = \frac{1}{2}\left( \sqrt{n} + \frac{n}{\sqrt{n}} \right ) = \frac{1}{2}\left( \sqrt{n} + \sqrt{n} \right )$ is true; thus, each successive approximation will always get closer and closer to $\sqrt{n}$.
Incidentally, the ancient Babylonians were aware of this formula for approximating the square root of two.
In searching for the integer square root, we take the floor of this approximation.
For example, suppose we are trying to find the integer square root of $298$. As our initial guess, let us start with $s_0 = 298$ (yes, a horrible guess); however, as we progress:
$s_1 = 149$, $s_2 = 75$, $s_3 = 39$, $s_4 = 23$ and $s_5 = 17$,
and this last value $s_5^2 = 289 \le 298$ while $(s_5 + 1)^2 =324 > 298$, so this is our integer square root of $298$.
Now, notice that once we calculated the value of $s_1$, we never again needed $s_0$, and once we calculated $s_2$, we never needed the value of $s_1$. We will use this in our implementation.
Additionally, integer arithmetic in C++ already performs the floor function. Thus, to find the integer square root, we need to know:
- an initial guess, and
- when to stop.
One guess of the integer square root is $n$, but this is wrong for all values of $n$ other than $0$ and $1$. Consequently, we will always start with this value, so why not just apply one application of our formula to an initial guess n:
If $s_0 = n$, $s_1 = \left \lfloor \frac{1}{2}\left( n + \frac{n}{n} \right ) \right \rfloor = \left \lfloor \frac{n + 1}{2} \right \rfloor$.
We are finished if (s*s <= n) && ((s + 1)*(s + 1) > n). Now, consider the statement:
"I am not both taller than five feet and shorter than six feet."
This is equivalent to saying:
"I am either not taller than five feet or not shorter than six feet."
However, being not taller than five feet is the same as being five feet or shorter, and being not shorter than six feet is the same as being six feet or taller. Thus, we may rephrase this as:
"I am either five feet or shorter or six feet or taller."
The first step is called the De Morgan's law: $\lnot (p \land q) \leftrightarrow (\lnot p) \lor (\lnot q)$. Thus,
while ( !((s*s <= n) && ((s + 1)*(s + 1) > n)) ) {
is equivalent to:
while ( !(s*s <= n) || !((s + 1)*(s + 1) > n) ) {
which is, in turn, equivalent to:
while ( (s*s > n) || ((s + 1)*(s + 1) <= n) ) {
Thus, our flow chart would look like
Finish implementing this function, and author a main() function that tests your isqrt(...) function. For example, isqrt(8) == 2 while isqrt(9) == 3.
unsigned int isqrt( unsigned int n ) {
unsigned int isqrt_n{(n + 1)/2}; // Our initial guess
while ( our condition is not satisifed ) {
// update isqrt_n
}
return isqrt_n;
}
As an aside, this iteration is guaranteed to converge. The number of iterations required is no more than $2.68 + \log_{3.44}(n)$ steps. For example, in calculating the integer square root of $2106534608$ (roughly two billion), the intermediate approximations are
$1053267304$, $526633653$, $263316828$, $131658417$, $65829216$, $32914623$, $16457343$, $8228735$, $4114495$, $2057503$, $1029263$, $515654$, $259869$, $133987$, $74854$, $51497$, $46201$, $45898$, $45897$
until we finally converge on the correct answer of $45896$ (where $45896^2 = 2106442816$ and $45897^2 = 2106534609$. Knowing how long an algorithm takes to find a solution is as important for an engineer as is the solution. An algorithm that takes a very long time to run is often not one that is useful or practical.
Next, we will look at a different type of algorithm for solving an algebraic equation.
Finding the $x$ such that $x = cos( x )$
Get out your calculator, first ensure that the angular units are radians and not degrees. Next, enter any number you want; for example, $2.5$, $-73.2$ or $53245322.323523$. Now, start pressing the cos key repeatedly. Suppose we start with 1 and press cos 88 times:
1
0.5403023058681397
0.8575532158463934
0.6542897904977792
0.7934803587425656
0.7013687736227565
0.7639596829006542
0.7221024250267077
0.7504177617637605
0.7314040424225099
0.7442373549005569
0.7356047404363473
0.7414250866101092
0.7375068905132428
0.7401473355678758
0.7383692041223232
0.7395672022122561
0.7387603198742113
0.7393038923969058
0.7389377567153445
0.7391843997714936
0.7390182624274122
0.7391301765296710
0.7390547907469174
0.7391055719265362
0.7390713652989450
0.7390944073790912
0.7390788859949921
0.7390893414033927
0.7390822985224023
0.7390870426953322
0.7390838469650002
0.7390859996481299
0.7390845495752127
0.7390855263619244
0.7390848683867142
0.7390853116067619
0.7390850130484204
0.7390852141609171
0.7390850786891230
0.7390851699445544
0.7390851084737987
0.7390851498812395
0.7390851219886894
0.7390851407774467
0.7390851281211138
0.7390851366465719
0.7390851309037208
0.7390851347721744
0.7390851321663375
0.7390851339216606
0.7390851327392538
0.7390851335357372
0.7390851329992164
0.7390851333606234
0.7390851331171754
0.7390851332811649
0.7390851331706996
0.7390851332451102
0.7390851331949862
0.7390851332287504
0.7390851332060064
0.7390851332213270
0.7390851332110069
0.7390851332179586
0.7390851332132759
0.7390851332164302
0.7390851332143054
0.7390851332157367
0.7390851332147726
0.7390851332154220
0.7390851332149846
0.7390851332152792
0.7390851332150808
0.7390851332152145
0.7390851332151244
0.7390851332151851
0.7390851332151442
0.7390851332151717
0.7390851332151532
0.7390851332151656
0.7390851332151572
0.7390851332151629
0.7390851332151591
0.7390851332151617
0.7390851332151599
0.7390851332151611
0.7390851332151603
0.7390851332151603
As you keep pressing cos, more and more digits stop changing, and it seems to converge to a fixed value. If you change the angular unit to degrees and try this again, you get a different result:
1
0.99984769515639124
0.99984774154521179
0.99984774153108381
0.99984774153108811
0.99984774153108811
It converges again, but to a different value, and it converges really quickly: after only pressing cos five times.
#include <cmath> // This is for std::cos
double cosine_convergence( double x );
double cosine_convergence( double x ) {
while ( x != std::cos( x ) ) {
x = std::cos( x );
}
return x;
}
We can now use this to see that it doesn't matter which value we start with, we will get the same answer:
#include <iostream>
#include <cmath> // This is for std::cos
int main();
double cosine_convergence( double x );
int main() {
std::cout << "Starting with 3.2: "
<< cosine_convergence( 3.2 ) << std::endl;
std::cout << "Starting with -4.5: "
<< cosine_convergence( -4.5 ) << std::endl;
std::cout << "Starting with 552323.23: "
<< cosine_convergence( 552323.23 ) << std::endl;
return 0;
}
// Definition of 'cosine_convergence'
If you run this, however, the output is...shall we say, disappointing:
Starting with 3.2: 0.739085 Starting with -4.5: 0.739085 Starting with 552323.23: 0.739085
Only six significant digits are being printed. The computer, however, uses many more than six significant digits. We can tell std::cout to print the full 16 digits that are stored (it's not really 16, as the number is stored in binary but it's close enough):
std::cout << "Starting with 3.2: " << std::setprecision( 16 ) << cosine_convergence( 3.2 ) << std::endl; std::cout << "Starting with -4.5: " << std::setprecision( 16 ) << cosine_convergence( -4.5 ) << std::endl; std::cout << "Starting with 552323.23: " << std::setprecision( 16 ) << cosine_convergence( 552323.23 ) << std::endl;
This, however, is in yet a different library:
#include <iomanip> // This is for std::setprecision
As you may guess, this library is an abbreviation of input-output manipulation. Thus, you must include this file, as well. Now the output is:
Starting with 3.2: 0.7390851332151607 Starting with -4.5: 0.7390851332151607 Starting with 552323.23: 0.7390851332151607
Now, here is something for you to try:
- What happens if you use std::cos( 0.5*x )?
- What happens if you use std::cos( 0.01745329251994330*x )?
- What happens if you use std::cos( 0.7*x )?
The first two appear to be okay—it converges to some value. However, at this point, you're probably a touch annoyed: each time you update a value, you must go back and change two places in the function cosine_convergence(...). If you forget to change it in one place, you run into an infinite loop. Thus, why not make that coefficient a parameter? For example,
#include <cmath> // This is for std::cos
double cosine_convergence( double x, coefficient c );
double cosine_convergence( double x, coefficient c ) {
while ( x != std::cos( c*x ) ) {
x = std::cos( c*x );
}
return x;
}
Now we can run:
std::cout << "Iterating cos( 0.5*x ): " << std::setprecision( 16 ) << cosine_convergence( 3.2, 0.5 ) << std::endl; std::cout << "Iterating cos( 0.01745329251994330*x ): " << std::setprecision( 16 ) << cosine_convergence( 3.2, 0.01745329251994330 ) << std::endl; std::cout << "Iterating cos( 0.7*x ): " << std::setprecision( 16 ) << cosine_convergence( 3.2, 0.7 ) << std::endl;
The first two produce an output, but the third seems to be a problem: nothing is being printed to the screen, but the program is still executing. Let us modify our code to log the values that are being calculated:
#include <cmath> // This is for std::cos
double cosine_convergence( double x, double coefficient );
double cosine_convergence( double x, double coefficient ) {
while ( x != std::cos( coefficient*x ) ) {
x = std::cos( coefficient*x );
std::cout << "x: " << std::setprecision( 16 ) << x << std::endl;
}
return x;
}
If you run this again with cosine_convergence( 3.2, 0.7 ), you get:
x: -0.6203616120126795 x: 0.9071845287932304 x: 0.805054393050757 x: 0.8453704353013128 x: 0.8299605850925448 x: 0.8359293508100089 x: 0.8336289452706264 x: 0.8345172639425491 x: 0.8341744894355629 x: 0.834306793666074 x: 0.8342557325219017 x: 0.8342754397725121 x: 0.8342678338074653 x: 0.8342707693298367 x: 0.8342696363677173 x: 0.8342700736338073 x: 0.8342699048712691 x: 0.8342699700050606 x: 0.8342699448667177 x: 0.8342699545688445 x: 0.834269950824315 x: 0.8342699522695137 x: 0.8342699517117402 x: 0.8342699519270125 x: 0.8342699518439283 x: 0.8342699518759946 x: 0.8342699518636186 x: 0.8342699518683951 x: 0.8342699518665516 x: 0.8342699518672632 x: 0.8342699518669885 x: 0.8342699518670945 x: 0.8342699518670537 x: 0.8342699518670694 x: 0.8342699518670633 x: 0.8342699518670657 x: 0.8342699518670648 x: 0.8342699518670651 x: 0.834269951867065 x: 0.8342699518670651 x: 0.834269951867065
At the end, it is bouncing between two values, two values that are insignificantly different—they only different in the 16th significant digit!
Important concept:
Engineers are not interested in exact answers: they are only interested in solutions
that are good enough for the problem they are working on. Good enough will change
depending on the problem: if you are designing a mirror to be used on a car rear-view mirror,
you only have to make sure that a human's eye will not notice a distortion; however, if
you are designing a mirror for the Hubble Space
Telescope, your requirements will be significantly more strict. No one will pay
you for an answer that is more precise than it has to be.
16 digits of precision would be enough to measure the circumference of the orbit of the center of the earth to within a tenth of a millimeter—the width of a thick human hair. If you wanted to find that value of $x$ such that $x = \cos(0.7x)$, you are welcome to look at this two-thousand digit approximation of the solution.
Instead of stopping when the two are equal, let's stop when they're close enough: suppose we decide to stop when $|x - \cos( x )| < 10^{-10}$. In this case, we must continue iterating while $|x - \cos( x )| \ge 10^{-10}$:
#include <cmath> // This is for std::cos
#include <iomanip> // This is for std::cos
double cosine_convergence( double x, double coefficient );
double cosine_convergence( double x, double coefficient ) {
while ( std::abs( x - std::cos( x ) ) >= 1e-10 ) {
x = std::cos( x );
std::cout << "x: " << std::setprecision( 16 ) << x << std::endl;
}
return x;
}
At this point, you are probably getting bored of modifying the code, so let's instead get information from the user:
#include <iostream>
#include <cmath> // This is for std::cos
#include <iomanip> // This is for std::setprecision
int main();
double cosine_convergence( double x, double coefficient );
int main() {
double c;
while ( true ) {
std::cout << "Enter a coefficient: ";
std::cin >> c;
std::cout << "Iterating x = cos(" << c << "*x) converges to "
<< std::setprecision( 10 )
<< cosine_convergence( 1.0, c ) << std::endl;
}
return 0;
}
// Definition of 'cosine_convergence'
If you run this again with std::cos( 1.5*x ), you get (starting with x == 3.2), you will see that it doesn't stop again, but this time for a different reason:
x: 0.08749898343944727 x: 0.9913992759811771 x: 0.08359972673191708 x: 0.9921477692285247 x: 0.08248086465080451 x: 0.9923562779825784 x: 0.0821691631845706 x: 0.992413868216025 x: 0.08208306964890315 x: 0.9924297366966927 x: 0.08205934722714282 x: 0.9924341062278271 x: 0.08205281503340675 x: 0.9924353091990008 x: 0.08205101666117245 x: 0.9924356403711588 x: 0.08205052157792903 x: 0.9924357315399964 x: 0.08205038528578053 x: 0.9924357566378954 x: 0.08205034776587024 x: 0.9924357635470969 x: 0.0820503374370128
It no longer appears to be converging to one number, but rather it appears to be bouncing between two different numbers: 0.9924357 and 0820503. If you keep watching, it ends up jumping between two values:
x: 0.08205033351347035 x: 0.9924357661716413 x: 0.08205033351347035 x: 0.9924357661716413
Now, we are most certainly in an infinite loop, and for different values again, it will start bouncing between four values, etc.
In reality, what would you do? In all likelihood, you would try this perhaps 100 times, and then give up and say that the sequence doesn't seem to converge. We could start with an integer that is initially zero, and then each time we run the body of the loop, we will increase this counter. If the counter ever reaches 100, we will throw an exception:
#include <cmath> // This is for std::cos
#include <iomanip> // This is for std::cos
#include <cassert> // This is for assert
double cosine_convergence( double x );
double cosine_convergence( double x, double coefficient ) {
int counter{0};
while ( std::abs( x - std::cos(coefficient*x) ) >= 1e-10 ) {
++counter;
assert( counter < 100 );
x = std::cos(coefficient*x);
std::cout << "x: " << std::setprecision( 16 ) << x << std::endl;
}
return x;
}
Now, if counter ever equals 100, the program will terminate. Later, we will see how to throw exceptions and then we will see how to properly deal with such exceptions.
Questions and practice:
1. In our implementation of cosine_convergence, we have two function calls of the form std::cos( coefficient*x ) each time the body of the loop is run: once in the condition, and another time in the body of the loop. Can you reduce this to one function call by using local variables?
2. The ancient Babylonians were aware that if you started with $x_0 = 1$, if you iterate $x_{k + 1} = \frac{n + x_k^2}{2x_k}$, this sequence would converge to $\sqrt{n}$ for any $n \ge 0$. For example, if $n = 2$, then $x_1 = {2 + 1^2}{2 \cdot 1} = 1.5$. Then, $x_2 = \frac{2 + 1.5^2}{2\times 1.5} = \frac{4.25}{3}= 1.41\overbar{6}$, although the Babylonians used base 60.
Implement a function that calculates this by starting with $x_0 = 1$. Stop when $|x_k^2 - n| < 10^{-10}$. If the user enters a value less than 0, return std::nan( "" ).
double iterative_sqrt( double n );
3. Write a function that returns the arithmetic series:
$$a + 2a + 3a + 4a + \cdots + na$$by explicitly adding the terms.
double arithmetic_series( double scalar_multiple, unsigned int n );
Question: do you want to calculate the series $1 + 2 + 3 + 4 + \cdots + n$ first, and then multiply the result by $a$, or do you want to calculate $a + 2a + 3a + 4a + \cdots + na$ directly? Which do you think is more efficient and why?
4. You should have learned the closed formula for an arithmetic series where
$1 + 2 + 3 + 4 + \cdots + n = \frac{n(n + 1)}{2}$.
Use this to write a function:
double arithmetic_fast( double scalar_multiple, unsigned int n );
5. Suppose you let n = 1000000 and you calculated arithmetic_series( 0.3, n ) and you saw that the time it took for this calculate was $0.308 {\rm ms}$. How long do you think it would take to calculate arithmetic_series( 0.3, 2*n )? How long would it take to calculate arithmetic_series( 0.3, 10*n )?
6. In contrast with Question 6, does the time it takes to calculate arithmetic_series( 0.3, n ) change if you change the value of n?
7. In our function gcd(...), we find the remainder twice with each iteration of the loop. Can you come up with a way of using a local variable so that you need only calculate the remainder once per loop?
8. In our function gcd(...), we find the remainder using the modulus operator. Implicit in this is the fact that we must be performing a long-division algorithm. In some cases, however, an engineer may be working on a system so small that the processor does not even implement integer long division. It is possible to calculate the gcd without using long division by observing that if $m > n$ and $d$ divides both, then $d$ must also divide $m - n$. This is Dijkstra's algorithm for calculating the gcd without using division:
- If $m = n$, then ${\rm gcd}(m,n) = m$,
- if $m > n$, ${\rm gcd}(m,n) = {\rm gcd}(m - n,n)$,
- otherwise $m < n$, so ${\rm gcd}(m,n) = {\rm gcd}(m,n - n)$.
Try to implement this program:
double gcd_no_div( unsigned int m, unsigned int n );
Write a program main() that tests your implementation by comparing the output of both gcd(...) and gcd_no_div(...).
If you need help, try to draw a flow chart first. If you really need help, look here. Note: if you jump to this link without even spending ten minutes on this problem, you are welcome to do so; however, please choose institutions other than the University of Waterloo for your undergraduate application. However, if you spend a solid ten minutes attempting to find this solution, but just don't get it, don't let this discourage you; try sleeping on it first before following the link—persistence and perseverance are valuable traits and will help you.
9. Write a function that returns the geometric series:
$$ar^0 + ar^1 + ar^2 + ar^3 + \cdots + ar^n = a + ar + ar^2 + ar^3 + \cdots + ar^n$$given a scalar multiplier $a$, a ratio $r$ and an upper limit $n$ by explicitly adding the terms.
double geometric_series( double scalar_multiplier, double ratio, unsigned int n );
You will still need a counter for this program, but you will not ever use its value. Instead it will only be used to ensure that you have added the right number of terms to result.
Note: If you tried using r^k to calculate the ratio $r^k$ for $k = 0, 1, \ldots, n$, you will will find you are not getting the answer you expect. Instead, you will have to start with a term $t$ that is initialized to $1$, and with each iteration of the loop, you will have to add this term to the result and then multiply it by $r$.
If you are having trouble with this program, you can always look at the flow chart here.
10. You should have learned the closed formula for an geometric series where
$1 + r + r^2 + r^3 + \cdots + r^n = \frac{1 - r^{n + 1}}{1 - r}$.
Use this to write a function:
double geometric_series_fast( double scalar_multiplier,
double ratio,
unsigned int n );
Note: Arithmetic and geometric series form a overall expectation of Section C Discrete functions of the Grade 11 curriculum in Ontario; therefore if you graduated from an Ontario secondary school, you must have seen these concepts as no professional secondary school teacher would choose to deliberately deviate from the required curriculum. Now is a good time to refresh your memory, as you will use these concepts throughout your engineering studies at university. You can read more on this on Wikipedia: arithmetic progression and geometric series.
11. If $m > n$, argue that the gcd algorithm will stop after a finite number of loops by recalling that if $r = m {\rm mod} n$ then $0 \le r < n$.
12. What happens with the gcd algorithm if $m < n$? How many more times will the loop in ${\rm gcd}(m, n)$ execute versus ${\rm gcd}(n, m)$?
13. Some programmers who have never seen the modulus operator before may try to begin using % in their documentation. What is the problem with this approach if those programmers are trying to communicate with the mathematical community in general? For example, in Matlab, % is the comment operator (like // in C++), and the function rem( 315, 23 ) returns the remainder; and in Maple, % is used to refer to the result of the most recent calculation, while the remainder may be calculated as either irem( 315, 23 ); or 315 mod 23; Conversely, would you ever begin writing $3 * 4$ instead of $3 \times 4$ or $3*x^2$ instead of $3x^2$ just because you are now a programmer? Also, would would you write $3*(x - 3)*(2*x - 1)$ instead $3(x - 3)(2x - 1)$? Which is easier to parse?
14. Let us test if our isqrt(...) function is correct: write a program that calculates the integer square root of every integer from 1 to 100000000 (100 million) and test whether or not it is true that isqrt(n)*isqrt(n) <= n and also that (isqrt(n) + 1)*(isqrt(n) + 1) > n. If both of these conditions are not true, print the warning:
Warning: isqrt(some-value) != the-answer-returned-by-your-implementation
15. In your previous program include the ctime library:
#include <ctime>
Next, add the following line as the first statement in main():
time_t t0 = time( nullptr );
This will store the number of seconds since midnight, January 1st, 1970, in a variable t0, our initial time. The type time_t is a signed integer type that is large enough to store the value returned by time(...).
Finally, just before the return 0; statement, call this function again and subtract off t0:
std::cout << "This took " << (time( nullptr ) - t0) << " s." << std::endl;
Department of Electrical and Computer Engineering
University of Waterloo
200 University Avenue West
Waterloo, Ontario, Canada N2L 3G1
Phone: +1 519 888 4567 extension 37023
Fax: +1 519 746 3077
Website maintained by Douglas Wilhelm Harder