We will use Newton's method to practice programming with the C11 thread library. Give a function f, say f(x) = x^3 - 1, where x^3 denotes the third power of x, Newton's method iteratively strives to find a root of f. Starting with a value x_0, one sets x_{k+1} = x_k - f(x_k) / f'(x_k). Convergence of this is unclear, and even if it does converge, it is unclear what root of f is approached.

Here is an example run of Newton's method for three different (complex) values of x_0.

real valuecomplex valuecomplex conjucate of previous one
x_04.00000000000000-3.00000000000000 + 1.00000000000000i-3.00000000000000 - 1.00000000000000i
x_12.68750000000000-1.97333333333333 + 0.686666666666667i-1.97333333333333 - 0.686666666666667i
x_21.83781773931855-1.25569408339907 + 0.505177535283118i-1.25569408339907 - 0.505177535283118i
x_31.32390198935733-0.705870398114610 + 0.462793270024375i-0.705870398114610 - 0.462793270024375i
x_41.07278230395299-0.284016534261040 + 0.737606186335893i-0.284016534261040 - 0.737606186335893i
x_51.00482620484690-0.585120941952759 + 0.849582155034797i-0.585120941952759 - 0.849582155034797i
x_61.00002314326798-0.501764864834321 + 0.859038339628915i-0.501764864834321 - 0.859038339628915i
x_71.00000000053559-0.499955299939949 + 0.866052541624328i-0.499955299939949 - 0.866052541624328i

One can experimentally find which points on the complex plane converge to which root. Below there are two pictures: In the left one every pixel is colored in such a way that its color corresponds to the root which it converges to. In the right hand side, every pixel is colored in such a way that its color corresponds to the number of iterations which are necessary to get close to a root of f. Note that this computation can be carried out for each pixel in parallel. There is no interdependence among the iterations.

Implementation

Implement in C using the C11 thread library a program called newton that computes similar pictures for a given functions f(x) = x^d - 1 for complex numbers with real and imaginary part between -2 and +2. The program should write the results as PPM files (see the description of PPM below) to the current working directory. The files should be named newton_attractors_xd.ppm and newton_convergence_xd.ppm, where d is replaced by the exponent of x that is used.

Your program should accept command line arguments

./newton -t5 -l1000 7
./newton -l1000 -t5 7

The last argument is d, the exponent of x in x^d - 1. The argument to -l gives the number of rows and columns for the output picture. The argument to -t gives the number of threads used to compute the pictures. For example, the above command computes images of size 1000x1000 for the polynomial x^7 - 1. It uses 5 threads for the Newton iteration.

Recall from Assignment 0 that command line arguments can, in general, be passed in arbitrary order. The degree is a positional argument, i.e. it will be the last argument whenever the program newton is invoked.

For this assignment we do not list recommended videos. It is also part of the assignment to origient yourself in the material that you have learned.

Implement details: Iteration

You have to run the iteration for each point until either of these criteria is met. In particular, the cut off for the number of iterations mentioned later may only be applied after completing the iteration.

Implement details: Writing to file

PPM files

PPM (Portable PixMap) files are an ASCII based image format. You can write and read them as a text file. Wikipedia contains a description with examples. In this assignment the header must look like

P3
L L
M

where L is the number of rows and columns and M is the maximal color value used in the body of the file.

Implementation details: Simplifications

Implementation strategy

The following is a possible implementation strategy. While it is not mandatory to follow it, I highly recommend any student without extensive programming experience to adhere to it.

Splitting the assignment into smaller tasks

This assignment comprises six largely independent tasks:

Parsing

You can adjust the code from Assignment 0 to this end. Use global variables to store the number of threads and lines to make them available to all threads.

Synchronization

A framework for synchronizing compute and write threads is presented in the video material.

When finishing these steps you have functional synchronization facilities as required for the assignment. It is helpful to

Data transfer

Recall that the result of a single computation in the video was represented by an array

int * result;

Since items correspond to rows in this implementation, each result has to include the attractor index and the convergence index for each pixel in a row.

Implement, allocate, and free two global arrays 

TYPE_ATTR ** attractors;
TYPE_CONV ** convergences;

and corresponding local arrays in the computation and writing threads:

TYPE_ATTR * attractor;
TYPE_CONV * convergence;

Writing

It is advantageous to implement the writing step prior to the computation step, since you will have accessible feedback once the writing is in place. To guarantee meaningful input values in the writing thread add the following initialization to your computing threads:

for ( size_t cx = 0; cx < nmb_lines; ++cx ) {
  attractor[cx] = 0;
  convergence[cx] = 0;
}

Make sure that the attractor and convergence value 0 indeed is assigned meaning in your encoding. Most likely this is the case, but it is important to check these assumptions to avoid unnecessary searching for errors in your code.

One might be tempted to implement the writing as

fprintf(attr_file, "%d %d %d ", COLOR_VALUE0, COLOR_VALUE1, COLOR_VALUE2);
fprintf(conv_file, "%d %d %d ", GRAY_VALUE, GRAY_VALUE, GRAY_VALUE);

but recall that fprintf is slow; fwrite is a better solution.

Checking

Before providing code for the iteration, the checking should be in place to guarantee that each step of the computation is carried out under the valid assumption that the norm of x be not too small.

There are three kinds of checks:

When implementing the arithmetic part of this step, it is helpful to first consider two questions:

The checking for the d+2 states (d the degree of the given polynomial) inside of a repeated computation like the iteration in this assignment requires an adapted programming pattern. Consider one similar to the following:

for ( conv = 0, attr = DEFAULT_VALUE; ; ++conv ) {
  if ( CHECK CONDITION ) {
    attr = VALUE;
    break;
  }
  if ( CHECK CONDITION ) {
    attr = VALUE;
    break;
  }
  for ( EXPRESSION )
    if ( CHECK CONDITION ) {
      attr = VALUE_NOT_EQUAL_TO_THE_DEFAULT_ONE;
      break;
    }
  if ( attr != DEFAULT_VALUE )
    break;

  // computation
}

Computation

It remains to implement the computation in order to complete the program. Since you can use functionality from complex.h and you can hardcode your formula, this step is largely about finding an expression for the iteration step that is as efficient as possible. Recall that using cpow is not an option.

When hardcoding the expression that you derive from this, the following syntax is convenient:

switch ( degree ) {
  case 1:
    STATEMENTS FOR DEGREE 1;
    break;
  case 2:
    STATEMENTS FOR DEGREE 2;
    break;
  case 3:
    STATEMENTS FOR DEGREE 3;
    break;
  // insert further cases

  default:
    fprintf(stderr, "unexpected degree\n");
    exit(1);
}

To understand this code, read up on

Going into depth: Further questions to consider

This assignment presents plenty of opportunities for improvement, both from the algorithmic point of view and the technical one. Here are some inspirations for further exploration.

Timing Goals

If you take this course as a student at Chalmers/GU, you can find the timing goals in the Assignment Timing Goals.