C++ For Mathematicians (2006) [eng]

416		C++ for Mathematicians
	10
	12
	13
	10
	13.3333

	Only cout << (4./3)*10 gives the proper output.
2.3 The output of the program looks like this:
	-2
	2
	-2
	2

The mathematician’s mod operation a mod n is usually only defined for n > 0 and a mod n is an element of {0,1,2,...,n − 1}. Thus, −5 mod 3 = 1, but C++ returns −2 for (-5)%3. In addition, C++ allows the modulus, n, to be negative.

2.4Both ’3’ and 3 are integer types (the first is a char which is considered to be an integer type). The result of this expression is false because the character ’3’ does not have the same value as the integer 3. Indeed, on many computers ’3’ has the value 51 because the character ’3’ is the 51st in the ASCII character set.

2.5The first occurrence of the expression a==b evaluates to false because, when this expression is reached, a and b hold different values. When a bool value of false is written to the screen, the number 0 is used.

Next, the expression a=b has the effect of copying b’s value into a and returns the value now held in common in a and b. Hence, 10 is printed.

Finally, when the second occurrence of a==b is evaluated, both a and b hold the value 10 and so the expression evaluates to true which is printed as 1 on the screen.

2.6The result of your experiment might depend on the computer on which it is run. The following is the typical result.

If the numerator and denominator are both int types, both 00 and 10 evaluate to zero.

However, if the numerator or the denominator (or both) are float then the

special values nan and inf are returned for 00 and 10 , respectively. Clearly inf stands for infinity whereas nan stands for not a number.

2.7400.

2.8(a) A variable name may not begin with a digit.

(b)The minus sign is an operator and may not be part of a variable’s name.

Answers

417

(c)The word double is a C++ keyword and may not be used as a variable name.

(d)A variable name may not include a colon.

(e)The ampersand is an operator and may not be part of a variable’s name.

(f)This begins with a digit, contains a minus sign, and is a number (equal to 0.01).

3.1d = 17, x = 6, and y = −1. Other answers for x and y are possible.

3.2If this procedure is invoked with n equal to −1 we fall into an infinite recursion.

There is also the issue that if a large value of n is passed to this procedure, the result would be too large to fit in a long, overflow would result, and the answer returned would be incorrect.

3.3In this implementation, we take the input argument to be a real number (type double) and the return value to be an int. Your choice may vary.

int signum(double x) { if (x < 0.) return -1; if (x > 0.) return 1; return 0;

}

3.4 Here is an iterative version.

long fibonacci(long n) { if (n < 0) return -1; if (n==0) return 1; if (n==1) return 1; long a = 1;

long b = 1; long c;

for(int k=2; k<=n; k++) { c = a+b;

a = b; b = c;

}

return c;

}

This is a recursive version.

long fibonacci(long n) { if (n < 0) return -1; if (n==0) return 1; if (n==1) return 1;

return fibonacci(n-1)+fibonacci(n-2);

}

418	C++ for Mathematicians

3.5The iterative version is significantly faster than the recursive version. When the iterative version computes Fn, it does n − 1 additions. However, when the recursive version computes Fn, it requests the calculation of Fn−1 and Fn−2. The calculation of Fn−1 also requests the computation of Fn−2, and this is wasted effort. The amount of work needed to calculate Fn grows exponentially with n.

When asked to calculate Fn, the recursive procedure presented here is called

2Fn − 1 times (your program may have different behavior). This is easy to show by induction. Let c(n) denote the number of times fibonacci is called

when asked to calculate Fn. For n = 0,1 we have c(n) = 1 and for n > 1 we note that c(n) = 1 + c(n −1) + c(n −2). Now one checks that 2Fn −1 satisfies this recurrence with the given initial conditions.

The values you were requested to find are

F20 = 10,946

F30 = 1,346,269

F40 = 165,580,141.

3.6 Here is a complete program.

#include <iostream> #include <cmath> using namespace std;

const float zeta2 = M_PI*M_PI/6.;

float zeta2up(long n) { float ans = 0.;

for (long k=1; k<=n; k++) { ans += 1./float(k)/float(k);

}

return ans;

}

float zeta2down(long n) { float ans = 0.;

for (long k=n; k>=1; k--) { ans += 1./float(k)/float(k);

}

return ans;

}
int main() {
long n = 1000000; // 10ˆ6
cout << "zeta2up	= " <<	zeta2up(n) << endl;
cout << "zeta2down = " <<		zeta2down(n) << endl;
cout << "truth	= " <<	zeta2 << endl;
return 0;
}

The output looks like this.

zeta2up	= 1.64473
zeta2down =		1.64493
truth	=	1.64493

Answers

419

Observe that zeta2down, which begins the sum with 1/N2, gives the more accurate result. The discrepancy is due to roundoff issues. A float holds only so many digits and (on my computer) the epsilon value for a float is around 10−7. This means that by the time the sum reaches k = 106, the value held in ans is too large for the addition of 10−12 to have any effect. Many terms at the end of the series have little or no effect on the sum; their contribution is lost. However, those tail terms do have a cumulative contribution on the actual sum and are needed to give an accurate result. By summing in the reverse order these small terms do not have their contributions lost.

3.7The important issue here is to use call by reference for the third and fourth arguments. Here is a complete program and its output.

#include <iostream> #include <cmath> using namespace std;

void xy2polar(double x, double y, double& r, double& t) { r = hypot(x,y);

t = atan2(x,y);

}

void polar2xy(double r, double t, double& x, double &y) {

x= r*cos(t);

y= r*sin(t);

}

int main() { double x,y,r,t;

x = -4.; y = -1.; xy2polar(x,y,r,t);

cout << "The point (" << x << "," << y

<<") in polar coordinates is ("

<<r << "," << t << ")" << endl;

r = 2.; t = 2.5; polar2xy(r,t,x,y);

cout << "The point in polar coordinates (" << r

<< "," << t << ") is (" << x << "," << y << ")" << endl;

return 0;

}

	The point (-4,-1) in polar coordinates is (4.12311,-1.81577)
	The point in polar coordinates (2,2.5) is (-1.60229,1.19694)
3.8 Here is a program. Change long to long long (or		int64) if that type is
	available to you.

#include <iostream> using namespace std;

420	C++ for Mathematicians

bool is_zero_one(long n) { long ones_digit = n % 10;

if (ones_digit > 1) return false; if (n == 0) return true;

return is_zero_one(n/10);

}

long find_zero_one_mult(long n) { for (long k=1; ; k++) {

long m = k*n;

if (m < 0) return -1; // overflow without success if (is_zero_one(m)) return m;

}

return 0;

}

int main() {

for(long k=1; k<=100; k++) {

cout << k << "\t" << find_zero_one_mult(k) << endl;

}

return 0;

}

The first procedure, is_zero_one, checks if its input argument consists entirely of 0s and 1s (in base ten) by a recursive algorithm. The second procedure, find_zero_one_mult, does its work by trying all positive multiples of the input argument until it succeeds. However, if the multiplies overflow (detected by testing if(m<0)) then we return −1 to signal failure.

Changing long to long long avoids the overflow problem, but then it takes much longer to find the proper multiple.

This program is not efficient. You may notice the computer pausing while it works to find the least multiple of 9 of the required form: 111111111. The program is not likely to find the least such multiple of 99, because it is

111111111111111111.

	18
		digits

Instead of testing successive multiples of n until we find the result we want, we can step through decimal numbers composed entirely of 0s and 1s and check if n is a divisor. However, this is trickier to program.

To step through the decimal values 1, 10, 11, 100, 101, 110, and so on, we realize that if we interpret these numbers as binary, then we are simply counting 1, 2, 3, 4, 5, and so on. If we had a procedure to convert binary numbers to the corresponding decimal number with the same digits, the task would be easy.

To this end, we create a procedure named bin2dec that performs the mapping

n = ∑ b j2 j → f (n) = ∑ b j10 j

j≥0

j≥0

Answers

421

where the b j are in {0,1}. There’s a nice recursive way to define this transformation:

f (n) = 10 f		2n	1		if n is even, and
10 f		n−2		+ 1	if n is odd

with base cases f (0) = 0 and f (1) = 1. [Note: The recursion can be written with a single algebraic expression like this: f (n) = 10 f ( n/2 ) + (n mod 2).]

With this idea in place, we present a second solution to this problem. This program runs much faster than the first.

#include <iostream> using namespace std;

long long bin2dec(long n) { if (n==0) return 0;

if (n==1) return 1;

return 10*bin2dec(n/2) + (n%2);

}

long long find_zero_one_mult(long long n) { for (long long k=1; ; k++) {

long long m = bin2dec(k);

if (m < 0) return -1; // overflow without success if (m%n == 0) return m;

}

return 0;

}

int main() {

for(long long k=1; k<=100; k++) {

cout << k << "\t" << find_zero_one_mult(k) << endl;

}

return 0;

}

4.1 Here is a program that does the job.

#include "uniform.h" #include <iostream> #include <cmath> using namespace std;

int main() { double x,y,u,v;

const long nreps = 100000000; double sum;

seed(); sum = 0;

for (long k=0; k<nreps; k++) { x = unif();

y = unif(); u = unif(); v = unif();

sum += sqrt( (x-u)*(x-u) + (y-v)*(y-v) );

league James Fill informs me that it is 151 proximately 0.521405.

422 C++ for Mathematicians

}

cout << "The average length of the segment is " << sum/nreps << endl;

return 0;

}

	This gives the following output.
	The average length of the segment is 0.521435

This naturally leads to the question: What is the analytic answer? My col-

√

2 + 2 + sinh−1 1 which is ap-

4.2Here is a program. Note the inclusion of the uniform.h header; to compile this program we need this file (call it buffon.cc) and the files uniform.h and uniform.cc.

#include <iostream> #include <cmath> #include "uniform.h" using namespace std;

/**

*This procedure simulates one drop of Buffon’s needle.

*It returns true if the needle crosses a line and false if not.

*/

bool drop() {
double x1;		// x-coord of one end point of the needle
double x2;		// x-coord of the other end point
double theta;		// orientation of the needle
x1 = unif();
theta =	2*M_PI*unif();
x2 = x1	+ cos(theta);

if (floor(x1) == floor(x2)) return false;

return true;

}

int main() {

long nreps; // number of times to drop the needle cout << "Enter number of needle drops -> ";

cin >> nreps;

seed();

long count = 0; // number of successful drops for (int k=0; k<nreps; k++) {

if (drop()) ++count;

}

	Answers	423
cout << count <<	" successes out of " << nreps
<< " drops"	<< endl;

cout << "frequency = " << double(count)/double(nreps) << endl;

}

	Here is a sample run of this program.
	Enter number of needle drops -> 100000000
	63653758 successes out of 100000000 drops
	frequency = 0.636538
	This is reasonably close to the theoretical value, π2 ≈ 0.63662.

4.3Here is some advice on how to approach this exercise. The procedures to generate points are random in a circle and in a triangle should be of the following form,

void point_in_circle(double& x, double& y); void point_in_triangle(double& x, double& y);

Using reference arguments enables these procedures effectively to return two values.

For point_in_circle, use a rejection algorithm. That is, generate a point (x,y) uniformly at random in the square [−1,1] × [−1,1]. If x2 + y2 ≤ 1, then return (x,y); if not, try again. This is a good opportunity to use the do-while control structure.

For point_in_triangle, it does not matter in which triangle the points are generated. (Reason: Whether four points lie at the vertices of a convex quadrilateral is invariant under invertible affine transformations.) For simplicity, use the triangle with vertices located at (0,0), (1,0), and (0,1). A rejection algorithm may be used (generate points in a square until finding one in the triangle). However, it is more efficient to find a point uniformly at random in [0,1] × [0,1] and if it lies in the wrong half of the square, reflect it across the line through (0,1) and (1,0).

To test if four points determine the corners of a convex quadrilateral, first write a procedure to see if the points (x1,y1) and (x2,y2) lie on opposite sides of the line through (x3,y3) and (x4,y4). They are if and only if the two triples

of points (x1,y1),(x2,y2),(x3,y3) and (x1,y1),(x2,y2),(x4,y4) have opposite orientation (see the figure).

(x4 , y4 )(x2 , y2 )

(x1, y1 ) (x3, y3 )

424	C++ for Mathematicians

Thus, the points (x1,y1) and (x2,y2) lie on opposite sides of the line through (x3,y3) and (x4,y4) if and only if

	1 x1 y1 1 x1 y1
det	1 x2 y2	1 x2 y2		< 0.
	1 x3 y3	· 1 x4 y4

Call this procedure something like opposite_sides and use it to build a procedure that tests if four points determine a convex quadrilateral.

Your experiments should yield the following results. For a triangle, P(K) = 23 . For a circle, P(K) = 1 − 1235π2 ≈ 0.70448.

4.4Despite the presence of trig functions, the first version is somewhat faster than the rejection method on my computer.

4.5The probability that a point (x,y,z) chosen uniformly at random in [−1,1]3

lies within distance one of the origin equals the volume of the ball of radius 1 divided by the volume of the cube: 43 π ÷ 8 = π/6 ≈ 0.5236.

In high-dimensional space, the unit ball occupies only a tiny fraction of the cube [−1,1]n. Therefore, the probability that a given iteration of the rejection algorithm succeeds is small and many iterations are required to generate a point within the unit ball.

An efficient way to generate a point x = (x1,x2,... ,xn) uniformly at random on the unit sphere in Rn is to assign to each xi a Gaussian N(0,1) random value and then to scale by 1/ x .

4.6In order for the procedure to remember the position of the particle from one invocation to the next, use a static variable. Here is the code.

#include "uniform.h"

int random_walk() { static int position = 0; if (unif() > 0.5) {

++position;

}

else { --position;

}

return position;

}

4.7Naturally, we need to use a static variable to remember the value between calls. However, a static variable declared in up is not available for use by down, or vice versa. (It is possible to use a global variable, defined outside the scope of any procedure, but this practice is frowned upon.) The solution is to put the static variable in another procedure! Here’s a solution.

Answers

425

int updown(int change) { static int value = 0; value += change; return value;

}

int up() {

return updown(1);

}

int down() {

return updown(-1);

}

5.1ϕ(100) = 4, ϕ(29) = 256, and ϕ(6!) = 192.

5.2The problem is with the declaration int vals[n];. Because the value of n cannot be determined when the program is compiled, this statement is illegal. Instead, declare vals with the statement int* vals; and then allocate space with vals = new int[n];. When vals is no longer needed, release the allocated memory with delete[] vals;.

5.3x ≡ 23 (mod 180).

5.5First, there is nothing wrong with using bool types for the array theSieve. The decision to use char is based on an oddity. Most compilers use one byte to house a char, but some use four bytes to house a bool. Use the expressions sizeof(char) and sizeof(bool) to learn the behavior on your computer.

If we are using an array that has only a million entries, then the difference between one and four bytes per cell is not important. But if the array has hundreds of millions (or a billion) entries, then the difference is important as the program may exhaust available memory.

In Chapter 8 we introduce the vector<bool> type that provides greater memory efficiency; this uses only one bit for each entry.

5.6Did you remember to declare your array to hold 21 values?

#include <iostream> using namespace std;

int main() { long fib[21]; fib[0] = 1; fib[1] = 1;

for(int k=2; k<=20; k++) { fib[k] = fib[k-1] + fib[k-2];

}

for (int k=0; k<=20; k++) {

cout << k << "\t" << fib[k] << endl;

}

return 0;

}