Primality Testing — Part 3. A Further Optimized Method

The previous two posts cover the naive approaches for primality testing and a method with precomputation. This post introduces a new method which could perform better than previous approaches and it can be used together with the precomputation.

The optimization is based on the observation that all prime numbers are in the form of 6k+1 or 6k+5, with 2, 3 and 5 being the only exceptions. This is because all integers can be expressed in 6k+i for some integer k, and i  is from the set {0, 1, 2, 3, 4, 5}. And 2 divides 6k+0, 6k+2 and 6k+4, 3 divides 6k+3.

Therefore, the primality test for a given integer n can be done by checking if n is divisible by 2 and 3 first, and then checking if n is divisible for all numbers of form 6k+1 or 6k+5 and less than sqrt(n).

The idea of this optimization is actually similar to the Sieve of Eratosthenes algorithm described in previous post. We eliminated testing divisibility for numbers that is divisible by the prime factors.

In general, all prime numbers are in the form gk + i, where i < k and it represents the numbers that are coprime to g. In the case above, g = 6 and the coprime of 6 are 1 and 5. If i and g are not coprime, then i must be divisible by some prime factors of g.

We go through another example, suppose g = 2*3*5 = 30, then all integers can be expressed in the form 30*k + i where i is in the set {0, 1, 2….29}. Now we mark the numbers that divisible by 2, 3 and 5 respectively. 30*k + {0, 2, 4, 6, 8, 10, 12, …28} are divisible by 2; 30*k + {3, 6, 9, 12, 15, 18, 21, 24, 27} are divisible by 3; and 30*k + {5, 10, 15, 20, 25} are divisible by 5. The rest of the numbers {1, 7, 11, 13, 17, 19, 23, 29} are not divisible by 2, 3 and 5 and are therefore coprimes of 30.

So to test if a number n is prime number, we first test if the number is divisible by 2, 3, and 5 first, and then check if the number is divisible by numbers in the form 30*k + {1, 7, 11, 13, 17, 19, 23, 29}  up to sqrt(n).

Below is an implementation of the optimization for primality testing when g is set to 6,

/**

generate prime numbers

**/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

 

int rdiv(int n, int m) {

    if (n%m ==0) {

        printf("%d is divisible by %dn", n, m);

        return 0;

    }

    return 1;

}

 

int prime(int n) {

    int i, k, ed, m;

    int cpr[2] = {1, 5};

    //check if we can divide the prime factors first

    if (n == 2 || n == 3) {return 1;}

    if (rdiv(n, 2) == 0) {

        return 0;

    } else if (rdiv(n, 3) == 0) {

        return 0;

    }

    //printf("ifprime: %d", ifprime);

    //check up to sqrt(n) of the form 6k + {coprime of 6} 

    ed = sqrt(n) + 1;

    for (k = 0; ;++k) {

        for (i = 0; i < 2; ++i) {

            m = 6*k + cpr[i];

            if (m < 3) {continue;}

            if (m > ed) {

                return 1;

            } 

            if (n%m == 0) {

                printf("%d is divisible by %dn", n, m);

                return 0;

            }

        }

    }

    return 1;

}

 

int main(int argc, char **argv) {

    int n = atoi(argv[1]);

    int i;

    struct timeval stTime, edTime;

    gettimeofday(&stTime, NULL);

    if (prime(n)) {

        printf("%d is prime numbern", n);

    } else {

        printf("%d is not prime numbern", n);

    }

    gettimeofday(&edTime, NULL);

    printf("time: %u:%un", (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec)); 

    return 0;

}

Save the code to prim4.c and compile the code with command below,

gcc -o prim4 prim4.c –lm

Below are some simple tests:

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim4 2

2 is prime number

time: 0:49

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim4 3

3 is prime number

time: 0:62

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim4 5

5 is prime number

time: 0:52

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim4 2147483647

2147483647 is prime number

time: 0:267

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim4 2147483641

2147483641 is divisible by 2699

2147483641 is not prime number

time: 0:72

Note that the code can be further improved by adding the pre-computation we mentioned in part2.

Primality Testing — Part 2. Pre-computation with Sieve of Eratosthenes

Previous post covers two naive methods for primality testing. This post introduces an algorithm to generate prime numbers called Sieve of Eratosthenes. With this algorithm, an optimization can be applied to speed up the primality testing.

Sieve of Eratosthenes

Suppose we want to generate prime numbers up to N. The algorithm consists of two loops, and it can be described as below,

  • Create a mask for all integers up to N.
  • In the first loop, let x be 2, mark 2x, 3x, 4x for all n such that nx < N.
  • finding an integer bigger than x and not marked, let it be the new x, and repeat the step above
  • the first loop terminates when we cannot find an integer bigger than x and not marked
  • In the second loop, starting from 2, collects the numbers that are not marked, which are the prime numbers up to N.

Below is an implementation of the Sieve of Eratosthenes algorithm,

/**

generate prime numbers

**/

#include <stdio.h>

#include <stdlib.h>

 

#define MAXV 10000000

#define ININ 1000

 

unsigned char bits[MAXV/8];

 

void setbit(int pos) {

    bits[pos/8] |= (0x01 << (pos%8));

}

 

int testbit(int pos) {

    return (bits[pos/8] & (0x01 << (pos%8)));

}

 

void genprime(int n, int **pr, int *npr) {

   int i, j, k;

   for (i = 2; i < n; ++i) {

       if (!testbit(i)) {

           for (j = 2; i*j < n; ++j) {

               setbit(i*j);

           }

       }

   } 

   k = ININ;

   *pr = malloc(sizeof(int) * k);

   *npr = 0;

   for (i = 2, j = 0; i < n; ++i) {

       if (!testbit(i)) {

           (*pr)[j++] = i;

           (*npr)++;

           if (j == k) {

               k += ININ;

               *pr = realloc(*pr, sizeof(int)*k);

           }

       }

   }

}

 

int main(int argc, char **argv) {

    int n = atoi(argv[1]);

    int *pr;

    int npr;

    int i;

    struct timeval stTime, edTime;

    if (n > MAXV) {

        printf("the max n is %dn", MAXV);

        return 0;

    }

    gettimeofday(&stTime, NULL);

    genprime(n, &pr, &npr);

    printf("No. of prime numbers: %dn", npr);

    for (i = 0; i < npr; ++i) {

        printf("%d  ", pr[i]);

    }

    printf("n");

    gettimeofday(&edTime, NULL);

    printf("time: %u:%un", (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec)); 

    if (pr != NULL) {

        free(pr);

    }

    return 0;

}

Note that the code above use a single bit to indicate if an integer is marked or not.

Save the code to prim2.c, and then compile the code with the command below,

gcc -o prim2 prim2.c

Below is a screenshot of some simple tests,


Figure 1. Simple Tests of Implementation for Sieve of Eratosthenes

Primality Testing with Sieve of Eratosthenes for Precomputation

One way to speed up primality testing is to precompute a list of prime numbers up to N. Then we first check if an input integer n is divisible by the prime numbers generated. If the number is not divisible by any of those prime numbers, we then proceed to the test mentioned in the naive approach post.

Below is an implementation of this precomputation approach,

/**

a primality test with precomputation

**/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

 

#define MAXV 10000000

#define ININ 1000

 

unsigned char bits[MAXV/8];

 

void setbit(int pos) {

    bits[pos/8] |= (0x01 << (pos%8));

}

 

int testbit(int pos) {

    return (bits[pos/8] & (0x01 << (pos%8)));

}

 

void genprime(int n, int **pr, int *npr) {

   int i, j, k;

   for (i = 2; i < n; ++i) {

       if (!testbit(i)) {

           for (j = 2; i*j < n; ++j) {

               setbit(i*j);

           }

       }

   } 

   k = ININ;

   *pr = malloc(sizeof(int) * k);

   *npr = 0;

   for (i = 2, j = 0; i < n; ++i) {

       if (!testbit(i)) {

           (*pr)[j++] = i;

           (*npr)++;

           if (j == k) {

               k += ININ;

               *pr = realloc(*pr, sizeof(int)*k);

           }

       }

   }

}

 

int prime(int n, int *pr, int npr, int x) {

    int i;

    int j;

    

    for (i = 0; i < npr; ++i) {

        if (n == pr[i]) {

            return 1;

        }

        if (n%pr[i] == 0) {

            printf("%d is divisble by %dn", n, pr[i]);

            return 0;

        }

    }

    if (n == 2) {

        return 1;

    }

    if (!(n&1)) {

        printf("%d is divisible by 2n", n);

        return 0;

    }

    j = sqrt(n);

    //to be safe, run to j+1

    i = (x+1) > 3 ? (x+1):3;

    if (!(i&1)) {

        ++i;   //make sure i starts as odd number

    }

    for (; i <= j + 1; i+=2) {

        if (n%i == 0) {

            printf("%d is divisible by %dn", n, i);

            return 0;

        }

    }

    return 1;

}

 

int main(int argc, char **argv) {

    int m, n;

    int *pr;

    int npr;

    int i;

    struct timeval stTime, stTime2, edTime;

    if (argc < 3) {

        printf("./prim3 <primality test number> <prime number generation limit>n");

        return 0;

    }

    m = atoi(argv[1]);

    n = atoi(argv[2]);

    if (n > MAXV) {

        printf("the max n is %dn", MAXV);

        return 0;

    }

    gettimeofday(&stTime, NULL);

    genprime(n, &pr, &npr);

    gettimeofday(&stTime2, NULL);

    //printf("No. of prime numbers: %dn", npr);

/*    for (i = 0; i < npr; ++i) {

        printf("%dn", pr[i]);

    }*/

    if (prime(m, pr, npr, n)) {

        printf("%d is prime numbern", m);

    } else {

        printf("%d is not prime numbern", m);

    }

    gettimeofday(&edTime, NULL);

    printf("total time: %u:%un", (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec)); 

    printf("test prime time: %u:%un", (unsigned int)(edTime.tv_sec - stTime2.tv_sec), (unsigned int)(edTime.tv_usec - stTime2.tv_usec)); 

    if (pr != NULL) {

        free(pr);

    }

    return 0;

}

Save the source code to prim3.c, and then compile the code with the command below,

gcc -o prim3 prim3.c –lm

Below are some simple tests for the compiled program,

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 2 10

2 is prime number

total time: 0:144

test prime time: 0:48

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 3 10

3 is prime number

total time: 0:60

test prime time: 0:22

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 2147483647 1000

2147483647 is prime number

total time: 0:469

test prime time: 0:302

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 2147483643 1000

2147483643 is divisble by 3

2147483643 is not prime number

total time: 0:223

test prime time: 0:57

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 2147483641 1000

2147483641 is divisible by 2699

2147483641 is not prime number

total time: 0:250

test prime time: 0:80

References:
1. Wikipedia: http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes

Primality Testing–Part 1. Naïve Approaches

Primality refers to determining if a given input number is prime or not. This post covers two naive approaches for primality testing.

Naive Approach 1

This most straightforward approach is to check if a given input integer n is divisible by integers from 2 to n-1. If a number is prime, then it should not divisible by any numbers in the range [2, n-1]. This is implemented as below,

/**

a naive method of determining if n is prime or not

divide from 2 to n-1

**/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

 

int prime(int n) {

    int i;

    for (i = 2; i <= n-1; i++) {

        if (n%i == 0) {

            printf("%d is divisble by %dn", n, i);

            return 0;

        }

    }

    return 1;

}

 

int main(int argc, char **argv) {

    int n = atoi(argv[1]);

    struct timeval stTime, edTime;

    gettimeofday(&stTime, NULL);

    if (prime(n)) {

        printf("%d is primen", n);

    } else {

        printf("%d is not primen", n);

    }

    gettimeofday(&edTime, NULL);

    printf("time: %u:%un", (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec));

}

Save the code to prim0.c. Compile the code with the command below,

gcc -o prim0 prim0.c –lm

Below are some sample testing results,

./prim0 2147483647

2147483647 is prime

time: 7:145569

 

./prim0 2147483640

2147483640 is divisble by 2

2147483640 is not prime

time: 0:46

 

./prim0 2

2 is prime

time: 0:51

 

./prim0 3

3 is prime

time: 0:49

Naive Approach 2: Improved Version of NA 1

The above algorithm can be improved. Firstly, if a number is divisible by an even integer, it must be divisible by 2. So we only need to check if the input number n is divisible by 2 and skip the rest of the even integers.

Secondly, if a number n is divisible by a number x, where x > sqrt(n). Then n = x*y, where y < sqrt(n), and y is also divisible by n. Therefore, we only need to check up to integer sqrt(n).

Below is the implementation with the two optimizations described,

/**

a naive method of determining if n is prime or not

divide from 2 to sqrt(n), skip the even numbers

**/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

 

int prime(int n) {

    int i, j;

    if (n == 2) {

        return 1;

    }

    if (!(n&1)) {

        //n is even

        printf("%d is divisible by 2n", n);

        return 0;

    }

    j = sqrt(n);

    //to be safe, run to j+1

    for (i = 3; i <= j + 1; i+=2) {

        if (n%i == 0) {

            printf("%d is divisible by %dn", n, i);

            return 0;

        }

    }

    return 1;

}

 

int main(int argc, char **argv) {

    int n = atoi(argv[1]);

    struct timeval stTime, edTime;

    gettimeofday(&stTime, NULL);

    if (prime(n)) {

        printf("%d is primen", n);

    } else {

        printf("%d is not primen", n);

    }

    gettimeofday(&edTime, NULL);

    printf("time: %u:%un", (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec));

}

Note that sqrt is a floating point operation. We check until sqrt(n) + 1, just to be safe.

Save the code to prim1.c, and compile the code using the command below,

gcc -o prim1 prim1.c –lm

Below are some testing results,

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim1 2

2 is prime

time: 0:50

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim1 3

3 is prime

time: 0:60

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim1 2147483647

2147483647 is prime

time: 0:308

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim1 2147483641

2147483641 is divisible by 2699

2147483641 is not prime

time: 0:73

References:
1. Wikipedia Primality Test: http://en.wikipedia.org/wiki/Primality_test
2. Algorithms by Dasgupta, CH Papadimitriou and UV Vazirani

gnuplot — Basic 3D Plot

Previous post covers gnuplot basic 2D plot. This post covers the basics of gnuplot 3D plot.
Most of the 2D plot options apply to 3D plot, one should read the 3D plot post before reading this post.

gnuplot supports splot command for drawing 3D figures and graphs.

There are 3 options specific to 3D plot: isosamples, hidden3d, and surface.

isosamples: controls the number of grid points at which a function will be evaluated when using splot. It is therefore only relevant when plotting functions.

set isosamples {int:xlines} [,{int:ylines}]

hidden3d: use opaque surface rendering. It has the following format:

set hidden3d [offset {int:offset}] [trianglepattern {int:maks}]

  • offset: used to control the color and line type used for the bottom side of the surface.
  • trianglepattern: control which lines will be drawn to connect neighboring grid points.

surface: by default, the surface is shown. It can be switched off by using the command,

unset surface

or switch on with the command,

set surface

An Example

Below is a sample 3D plot script,

unset title

#unset key

#unset xtics

#unset ytics

#set format y "%.1t";

#set format x "%.1t";

set title "Test 3D" 0,1

set xlabel "X" -4,-2;

set ylabel "Y" 0,-2;

set zlabel "Z" 4,6;

#set border 0

set hidden3d offset 0

set surface

#set the view angle: one can adjust manually to get the proper value first

#and then modify the value here

set view 80,49 

#set xrange []

#set yrange []

#set zrange []

#set ticslevel 0

#splot '45050.data' u 1:2:3 with lines

splot '$0' using 1:2:3 title "data1" with linespoints, '$1' using 1:2:3 title "data2" with linespoints

#

set term push 

set term png size 800, 600

set output "$2"

replot

set output

set terminal pop

Suppose the file is saved as 3d.pg. And the two data files “3data.dat” and “3data2.dat” are listed below. Note that the data file must follow certain formats, otherwise the plot may not be expected. gnuplot supports two 3D data file formats, we only covers the grid format here. One can refer to reference 1 for more details.

The grid data format supported by gnuplot requires the follows,

  • each line must contain x and y coordinates and the z value.
  • data are organized into data blocks. And blocks are separated by a single blank line.
  • each block must contain all data points for a single row of data points parallel to x axis. Therefore, the x coordinate stays constant at each data block.
  • all data blocks must contain same number of data points. If the number is not the same, then gnuplot cannot draw a proper surface.

3data.dat

0 0 29.21

0 160 28.94

0 320 28.87

0 480 28.28

0 640 28.83

0 800 27.94

0 960 27.02

 

160 0 28.45

160 160 28.48

160 320 27.73

160 480 27.72

160 640 26.98

160 800 26.72

160 960 26.38

 

320 0 27.30

320 160 26.95

320 320 26.80

320 480 26.85

320 640 26.52

320 800 26.11

320 960 25.95

 

480 0 25.94

480 160 25.86

480 320 25.84

480 480 26.35

480 640 25.61

480 800 25.13

480 960 24.37

 

544 0 25.67

544 160 25.49

544 320 25.50

544 480 25.60

544 640 25.27

544 800 25.63

544 960 24.01

And 3data2.dat

0 0 17.67

0 160 17.67

0 320 17.67

0 480 17.67

0 640 17.67

0 800 17.67

0 960 17.67

 

160 0 17.67

160 160 17.67

160 320 17.67

160 480 17.67

160 640 17.67

160 800 17.67

160 960 17.67

 

320 0 17.67

320 160 17.67

320 320 17.67

320 480 17.67

320 640 17.67

320 800 17.67

320 960 17.67

 

480 0 17.67

480 160 17.67

480 320 17.67

480 480 17.67

480 640 17.67

480 800 17.67

480 960 17.67

 

544 0 17.67

544 160 17.67

544 320 17.67

544 480 17.67

544 640 17.67

544 800 17.67

544 960 17.67

One can execute the script using the command below,

gnuplot> call “./3d.pg” “./3data.dat” “./3data2.dat” “3test.png”

And the plotted figure will be as below,


Figure 1. a 3-D Plot using gnuplot

References:
1. gnuplot In Action

gnuplot–Basic 2D Plot

This is the third project that I needed to use gnuplot, and I’m learning it for the third time. I know I’ll need to write a note for it. So here it is.

gnuplot is easy to start with. You can plot a nice graph after reading a tutorial for 5 minutes. This post summarizes some basic stuff for 2D plotting.

gnuplot supports both interactive mode and scripting mode. We’ll cover both.

A Simple Plot

The plot command is used to plot 2D figures and graphs. It’s very simple to use. For example, you want to plot a sine function and a y=x function. First, you enter the “gnuplot” command in the Linux command shell to go the gnuplot command interface.

And then you enter the command below,

gnuplot> plot sin(x), x

This gives you the figure below,

Figure 1. a sine function and a y=x function drawn by gunplot

 

Plotting for a Data File

gnuplot is easy to integrate with other programs. Another program can output a data file and gnuplot plot the graph based on the data file.

Suppose we have a data file named “data.dat” with content as below,

0.1 31.15

0.2 30.42

0.3 29.94

0.4 27.65

0.5 25.76

0.6 23.61

0.7 21.16

0.8 20.13

0.9 18.30

1.0 16.81

To plot the graph based on the data file, we can use the command below,

gnuplot>  plot “data.dat” using 1:2 title “1-3” with lines

And we get the figure below,

Figure 2. gnuplot plots figure based on a data file

Note that “using 1:2” maps the first column of the data to x axis and the second column of the data to y axis. title “1-3” specifies the legend title for the line. “with line” specifies the drawing styles. Other common options of styles include “points” and “linespoints”.

Also note that it is easy to draw two lines with data from the same file. Suppose the data.dat file contains the third column, then the following command can give us two lines.

gnuplot> plot “data.dat” using 1:2 title “1-3” with lines, “data.dat” using 1:3 title “1-3” with points

 

Plotting for Multiple Files

It is easy to plot from multiple data files too. Suppose we have another file “data2.dat” with the content as below,

0.1 17.67

0.2 17.67

0.3 17.67

0.4 17.67

0.5 17.67

0.6 17.67

0.7 17.67

0.8 17.67

0.9 17.67

1.0 17.67

gnuplot> plot “data.dat” using 1:2 title “1-3” with lines, “data2.dat” using 1:2 with linespoints

This gives us the figure below,

Figure 3. gnuplot plot from two data files

Export Current Configuration as Script

Once you finished plotting a figure, you may want to save the command and configuration for later usage. You can do it using the command below,

gnuplot> save “test.gp”

Suppose you just finished the example above, and the save command will give you a file as below,

#!/usr/bin/gnuplot -persist

#

#    

#        G N U P L O T

#        Version 4.2 patchlevel 6 

#        last modified Sep 2009

#        System: Linux 2.6.32-40-generic

#    

#        Copyright (C) 1986 - 1993, 1998, 2004, 2007 - 2009

#        Thomas Williams, Colin Kelley and many others

#    

#        Type `help` to access the on-line reference manual.

#        The gnuplot FAQ is available from http://www.gnuplot.info/faq/

#    

#        Send bug reports and suggestions to <http://sourceforge.net/projects/gnuplot>

#    

# set terminal wxt 0

# set output

unset clip points

set clip one

unset clip two

set bar 1.000000

set border 31 front linetype -1 linewidth 1.000

set xdata

set ydata

set zdata

set x2data

set y2data

set timefmt x "%d/%m/%y,%H:%M"

set timefmt y "%d/%m/%y,%H:%M"

set timefmt z "%d/%m/%y,%H:%M"

set timefmt x2 "%d/%m/%y,%H:%M"

set timefmt y2 "%d/%m/%y,%H:%M"

set timefmt cb "%d/%m/%y,%H:%M"

set boxwidth

set style fill  empty border

set style rectangle back fc lt -3 fillstyle  solid 1.00 border -1

set dummy x,y

set format x "% g"

set format y "% g"

set format x2 "% g"

set format y2 "% g"

set format z "% g"

set format cb "% g"

set angles radians

unset grid

set key title ""

set key inside right top vertical Right noreverse enhanced autotitles nobox

set key noinvert samplen 4 spacing 1 width 0 height 0 

unset label

unset arrow

set style increment default

unset style line

unset style arrow

set style histogram clustered gap 2 title  offset character 0, 0, 0

unset logscale

set offsets 0, 0, 0, 0

set pointsize 1

set encoding default

unset polar

unset parametric

unset decimalsign

set view 60, 30, 1, 1  

set samples 100, 100

set isosamples 10, 10

set surface

unset contour

set clabel '%8.3g'

set mapping cartesian

set datafile separator whitespace

unset hidden3d

set cntrparam order 4

set cntrparam linear

set cntrparam levels auto 5

set cntrparam points 5

set size ratio 0 1,1

set origin 0,0

set style data points

set style function lines

set xzeroaxis linetype -2 linewidth 1.000

set yzeroaxis linetype -2 linewidth 1.000

set zzeroaxis linetype -2 linewidth 1.000

set x2zeroaxis linetype -2 linewidth 1.000

set y2zeroaxis linetype -2 linewidth 1.000

set ticslevel 0.5

set mxtics default

set mytics default

set mztics default

set mx2tics default

set my2tics default

set mcbtics default

set xtics border in scale 1,0.5 mirror norotate  offset character 0, 0, 0

set xtics autofreq  norangelimit

set ytics border in scale 1,0.5 mirror norotate  offset character 0, 0, 0

set ytics autofreq  norangelimit

set ztics border in scale 1,0.5 nomirror norotate  offset character 0, 0, 0

set ztics autofreq  norangelimit

set nox2tics

set noy2tics

set cbtics border in scale 1,0.5 mirror norotate  offset character 0, 0, 0

set cbtics autofreq  norangelimit

set title "" 

set title  offset character 0, 0, 0 font "" norotate

set timestamp bottom 

set timestamp "" 

set timestamp  offset character 0, 0, 0 font "" norotate

set rrange [ * : * ] noreverse nowriteback  # (currently [0.00000:10.0000] )

set trange [ * : * ] noreverse nowriteback  # (currently [-5.00000:5.00000] )

set urange [ * : * ] noreverse nowriteback  # (currently [-5.00000:5.00000] )

set vrange [ * : * ] noreverse nowriteback  # (currently [-5.00000:5.00000] )

set xlabel "" 

set xlabel  offset character 0, 0, 0 font "" textcolor lt -1 norotate

set x2label "" 

set x2label  offset character 0, 0, 0 font "" textcolor lt -1 norotate

set xrange [ * : * ] noreverse nowriteback  # (currently [-10.0000:10.0000] )

set x2range [ * : * ] noreverse nowriteback  # (currently [-10.0000:10.0000] )

set ylabel "" 

set ylabel  offset character 0, 0, 0 font "" textcolor lt -1 rotate by 90

set y2label "" 

set y2label  offset character 0, 0, 0 font "" textcolor lt -1 rotate by 90

set yrange [ * : * ] noreverse nowriteback  # (currently [-10.0000:10.0000] )

set y2range [ * : * ] noreverse nowriteback  # (currently [-10.0000:10.0000] )

set zlabel "" 

set zlabel  offset character 0, 0, 0 font "" textcolor lt -1 norotate

set zrange [ * : * ] noreverse nowriteback  # (currently [-10.0000:10.0000] )

set cblabel "" 

set cblabel  offset character 0, 0, 0 font "" textcolor lt -1 rotate by 90

set cbrange [ * : * ] noreverse nowriteback  # (currently [-10.0000:10.0000] )

set zero 1e-08

set lmargin  -1

set bmargin  -1

set rmargin  -1

set tmargin  -1

set locale "C"

set pm3d explicit at s

set pm3d scansautomatic

set pm3d interpolate 1,1 flush begin noftriangles nohidden3d corners2color mean

set palette positive nops_allcF maxcolors 0 gamma 1.5 color model RGB 

set palette rgbformulae 7, 5, 15

set colorbox default

set colorbox vertical origin screen 0.9, 0.2, 0 size screen 0.05, 0.6, 0 front bdefault

set loadpath 

set fontpath 

set fit noerrorvariables

GNUTERM = "wxt"

plot "data.dat" using 1:2 title "1-3" with lines, "data2.dat" using 1:2 with linespoints

#    EOF

Next time you want to draw the figure again, simple execute “load” or “call” command,

gnuplot> load “test.pg”

or

gnuplot> call “test.pg”

The differences between load and call is that call allows one to pass up to 10 parameters, which corresponds to $0 – $9 in the script file.

Saving the Drawing as Pictures

gnuplot allows one to export the plot in many formats, including “png”, “svg”, “postscript” etc. One can check out all the formats supported by “set term” command.

Below is the commands that exports the plot as png file.

gnuplot> set term png
gnuplot> set output “test.png”
gnuplot> replot

The first command set the output device type, the second command sets the output file and the third command replot the plot on the output file.

A more general script can be written as below to export plot,

set terminal push   # save the current terminal settings

set terminal png    # change terminal to PNG

set output "$0"     # set the output filename to the first option

replot              # repeat the most recent plot command

set output          # restore output to interactive mode

set terminal pop    # restore the terminal

Assume the script file is “test.gp”, then we can execute the command below in gnuplot command line interface to generate the graph.

call  “test.gp” “test.png”

A Simple Script

Now we provide a simple script,

unset title

#unset key

#unset xtics

#unset ytics

#set format y "%.1t";

#set format x "%.1t";

set title "Testndata.dat VS. data2.dat" 0,1

set xlabel "X";

set ylabel "Y";

#set border 0

#set xrange []

#set yrange []

plot '$0' using 1:2 title "data.dat" with linespoints, '$1' using 1:2 title "data2.dat" with linespoints

#

set term push

set term png size 800, 600

set output "$2"

replot

set output

set terminal pop

This script accepts three input parameters, $0 and $1 are two data files, and $2 is the output png filename.

Suppose the script is named “test.pg”, and we still use the two data files “data.dat” and “data2.dat”, and we output to “test.png” file. Then we can enter the command below,

gnuplog> call “test.pg” “data.dat” “data2.dat” “test.png”

Then we get the figure below,

Figure 4. gnuplot plot with a script

Once we’re done with gnuplot, we can type “quit” or “exit” to exit from the gnuplot command line interface.

Linux Signal–Part 2. Signal Handler

Previous post covers the basics of Linux signals. This post illustrates how to install a Linux signal handler with examples.

Linux provides two methods to install signal handler, signal and sigaction. sigaction is the preferred method over signal because signal behaves differently on different UNIX and UNIX-like OSes and therefore less compatible than sigaction. This post covers sigaction only.

Before going to sigaction system call, we go through several functions which we’ll need later. You can also go the end of the post to read the source code first and refer to the explaination later.

sigprocmask System Call

A signal can be blocked. In this case, it is not delivered to the process until it is unblocked. After the signal is generated and before it is delivered, the signal is in pending state. sigprocmask system call is used to query blocked signals. It has the following prototype,

#include <signal.h>

int sigprocmask(int how, const sigset_t *set, sigset_t *oldset);

how: specify how to change the block signal set. It can be one of the three values.  SIG_BLOCK: the set of blocked signal is the union of current set and the input parameter set SIG_UNBLOCK: the signals specified in set is removed from the set of blocked signals SIG_SETMASK: the set of blocked signals is set to same as the input parameter set

set: if not null, the set value is used to according to the description for how. If null, the blocked signal set is unchanged and how has no meeting.

oldset: if not null, the previous value of the signal mask is returned through oldset.

Signal Sets System Calls

The sigprocmask function returns a data structure of sigset_t type, several functions are defined to manipulate the signal set.

#include <signal.h>

int sigemptyset(sigset_t *set);

int sigfillset(sigset_t *set);

int sigaddset(sigset_t *set, int signum);

int sigdelset(sigset_t *set, int signum);

int sigismember(const sigset_t *set, int signum);

One can tell what the functions do from their names.

pause and sigsuspend System Calls

Two system calls are defined to suspend the execution of a process to wait for a signal to be caught, pause() and sigsuspend(). Their prototypes are as below,

#include <unistd.h>

int pause(void);

 

#include <signal.h>

int sigsuspend(const sigset_t *mask);

pause: Suspends execution until any signal is caught. It only returns when a signal was caught and the signal catching function returned. In this case, it returns -1 and errno is set to EINTR.

sigsuspend: Temporarily changes the signal mask and suspends execution until one of the unmasked signals is caught. The calling process’s signal maks is temporarily replaced by value given in mask input parameter until the sigsuspend returns or the process is terminated.

Note that sigsuspend always returns -1, normally with the error EINTR.

Now we explain the sigaction system call.

sigaction System Call

sigaction is a system call to query and change the actions taken when a process receives a particular signal. The sigaction system call has the following prototype,

#include <signal.h>

int sigaction(int signum, const struct sigaction *act, struct sigaction *oldact);

signum: specify the signal
act: if not NULL, the action specified by act is installed.
oldact: if not NULL, the previous action is returned.

And the data structure sigaction is defined as below,

struct sigaction {

      void     (*sa_handler)(int);

      void     (*sa_sigaction)(int, siginfo_t *, void *);

      sigset_t   sa_mask;

      int        sa_flags;

      void     (*sa_restorer)(void);

};

The parameters have the following meanings,

sa_handler: specifies the action taken to the signal, SIG_DFL (default action), SIG_IGN (ignore the signal) or a function pointer to a user defined signal handler. The defined signal handler should accepts the signal number as the only input argument.

sa_sigaction: if SA_SIGINFO flag is set in sa_flags, the handler pointed by sa_sigaction should be used.

The handler function should accept three arguments, the signal number, a pointer to an instance of siginfo_t structure (contains information about the signal) and a pointer to an object of type ucontext_t (contains the receiving process context which has been interrupted by the delivered signal. It is a void pointer can be casted to ucontext_t pointer).

The siginfo_t data structure has the following elements,

siginfo_t {

   int      si_signo;    /* Signal number */

   int      si_errno;    /* An errno value */

   int      si_code;     /* Signal code */

   int      si_trapno;   /* Trap number that caused

                            hardware-generated signal

                            (unused on most architectures) */

   pid_t    si_pid;      /* Sending process ID */

   uid_t    si_uid;      /* Real user ID of sending process */

   int      si_status;   /* Exit value or signal */

   clock_t  si_utime;    /* User time consumed */

   clock_t  si_stime;    /* System time consumed */

   sigval_t si_value;    /* Signal value */

   int      si_int;      /* POSIX.1b signal */

   void    *si_ptr;      /* POSIX.1b signal */

   int      si_overrun;  /* Timer overrun count; POSIX.1b timers */

   int      si_timerid;  /* Timer ID; POSIX.1b timers */

   void    *si_addr;     /* Memory location which caused fault */

   int      si_band;     /* Band event */

   int      si_fd;       /* File descriptor */

}

sa_mask: specify the signals should be blocked when the handler is in execution. In other words, the sa_mask adds signals to the signal mask of the process before signal handler is called. And when the signal handler function returns, the signal mask of the process is reset to its previous value. In this way, we can block certain signals when the signal handler is running. In addition, the signal that caused the signal handler to execute is also blocked, unless the SA_NODEFER flag is set.

sa_flags: specifies various options of handling the signal. The actions can be aggregated by the bitwise OR operation. The details of the flags can be referred from Linux man page.

sa_restorer: This is obsolete and should not be used.

An Example

Below is an example of installing a naive signal handler for SIGINT.

#include <stdio.h>

#include <signal.h>

 

void check_block(int signum) {

    sigset_t set;

    if (sigprocmask(0, NULL, &set)==-1) {

        perror("error sigprocmask:");

    } else {

        if (sigismember(&set, signum)) {

            printf("%d is blockedn", signum);

        } else {

            printf("%d is not blockedn", signum);

        }

    }

}

 

void my_handler(int signum) {

    check_block(signum);

    printf("inside signal handler: %dn", signum);

}

 

void my_handler2(int signum, siginfo_t *info, void *context) {

    check_block(signum);

    printf("inside signal handler: %dn", info->si_signo);

}

 

int main() {

    struct sigaction action;

//    action.sa_handler = my_handler;

//    action.sa_flags = SA_RESTART;

    action.sa_sigaction = my_handler2;

    action.sa_flags = SA_RESTART | SA_SIGINFO;

    sigaction(SIGINT, &action, NULL);

    printf("press cntl + cn");

    pause();

    printf("after signal handler:n");

    check_block(SIGINT);

    return 0;

}

The code first installs the signal handler for SIGINT with sigaction system call. Inside the signal handler, it checks if the SIGINT signal is blocked. As described in sa_mask of sigaction system call, the SIGINT  should be blocked inside the signal handler.

After installing the signal handler, the program calls pause() system call to wait for signal to occur. And it checks to see if the SIGINT is blocked again after the execution returns from the signal handler. This time, SIGINT should not be blocked.

Compile the code with command below,

gcc -o sigaction sigaction.c

And a sample execution is as below,

Untitled

Figure 1. Sample Execution of Sigaction Program

SA_RESTART and Interrupted System Call

You may notice that the we set the sa_flags with SA_RESTART in the example above. This has something to do with system calls.

When a signal is caught while a system call is blocked (e.g. waiting for IO), then two results can happen after the signal handler:

  • the call is restarted automatically after the signal handler returns
  • the call fails with errno set to EINTR

With SA_RESTART flag set, some system calls are restarted automatically after the signal handler, including open, wait, waitpid etc. Those system calls return failure if SA_RESTART is not set.

There’re some system calls return failure regardless SA_RESTART is set or not. A special case is the sleep function, which is never restarted but the number of seconds remaining is returned instead of failure.

The detailed list of system calls that are restarted can found with “man 7 signal” command.

The Async-signal-safe Functions

When a signal is caught by a signal handler, the normal execution is temporarily interrupted by the signal handler. If  the signal handler returns instead of terminating the process, the execution continues at where it is interrupted. Special attention needs to be paid for functions inside signal handler because we don’t want signal handler code to interfere with other code.

A bad example of this is the signal handler changes the global static variable which is used by the code after the handler.

POSIX standard defines a list of safe functions which can be called inside the signal handler, the detailed list can be obtained with “man 7 signal”.

Note that in the example above, we called printf() function inside the signal handler, which is actually not safe. If the signal is caught when we are calling printf inside the main function, the results may be unexpected.

References:

1. The Linux Signals Handling Model: http://www.linuxjournal.com/article/3985

2. Linux Signals for the Application Programmer: http://www.linuxjournal.com/article/6483

3. Linux sigaction man page.

4. Linux signal man page (man 7 signal).

5. Advaned Programming in the UNIX environment

Linux Signal–Part 1. The Basics

Linux Signals Overview

Linux supports both POSIX reliable/standard signals and real-time signals. The first 31 signals are standard signals. Real time signals ranges from 34 to 64. The Linux command “kill -l” lists all signals numbers and names.

This post discusses reliable signals only.

A signal can be synchronous or asynchronous to the process, depending what caused the signal. Synchronous signals are also referred as traps, because they cause a trap into a kernel trap handler. An example is signal caused due to illegal instruction. Asynchronous signals are also referred as interrupts, they’re external to the current execution context and often used to send asynchronous events to a process.

If a process is in interruptible sleep state, the kernel can deliver the signal to the process and wake it up to handle it. For example, the process is waiting for terminal IO. If a process is in uninterruptible sleep, the kernel will hold the signal until the process wakes up. For instance, the process is waiting for disk IO.

Linux kernel maintains signal information for each process. The information includes an array of signal handlers and related information. Once a signal is generated, the kernel sets a bit in corresponding to the signal. Since it’s a single bit, multiple occurrences and a single occurrence are equivalent.

Signals names always start with SIG, they are defined by positive integer constants called signal numbers in the header file signal.h. signal.h is normally found in /usr/include/ directory of Linux. The actual signal numbers are usually defined in /usr/include/bits/signum.h.

Signals can happens at random times and the process can tell kernel to do one of the three things when a signal occurs,

1. ignore the signal. All except two signals (SIGSTOP and SIGKILL) can be ignored. SIGKILL always terminate a process, while SIGSTOP always clears any pending/queued SIGCONT signals and stop the process (a stopped is process may be resumed later, which is different from a terminated process). In addition, some hardware generated signals may cause the process behavior undefined if ignored (e.g. SIGSEGV).

2. catch the signal. All except two signals (SIGSTOP and SIGKILL) can be caught. We tell the kernel to execute a handler function whenever the signal occurs.

3. apply default action. Every signal has a default action. There’re four possible actions if the signal is not ignored or caught.

  • ignore: nothing happens
  • terminate: the process is terminated
  • terminate + coredump: create a core dump file and then terminate the process
  • stop: stop all threads in the process. The process goes to TASK_STOPPED state.

Linux Signals One by One

SIGHUP(1): Hangup. The signal is sent to the controlling process when a disconnect is detected by its controlling terminal interface. By default, the process is terminated.

SIGINT(2): Interrupt. The signal is sent to all foreground processes when the interrupt key (often DELETE or Control – C) is pressed. By default, the foreground processes are terminated.

SIGQUIT(3): Quit. The signal is sent to all foreground processes when the quit key (often Control – ) is pressed. In addition to terminates the foreground processes, it creates a core dump file.

SIGILL(4): Illegal instruction. It is generated when a process has executed an illegal (malformed, privileged, or unknown) hardware instruction. By default, the process is terminated and a core dump is created.

Below is a program that generates SIGILL,

typedef void(*FUNC)(void);

int main(void) {

   const static unsigned char insn[4] = {0xff, 0xff, 0xff, 0xff};

   FUNC function = (FUNC) insn;

   function();

   return 0;

}

SIGTRAP(5): Trace trap. It is used as a mechanism to notify a debugger when the process execution hits a breakpoint. By default, the process is killed and a core dump is created.

SIGABRT(6): Abort. It is generated by calling the abort function. The process terminates abnormally and a core dump file will be created by default.

SIGBUS(7): Bus error, BUS refers to the address bus in the context of a bus error. It indicates implementation defined hardware fault. It is usually caused by improper memory handling. By default, the process is terminated and a core dump is created.

SIGFPE(8): Floating point exception. It is sent to a process when it performs an illegal arithmetic operation. Note that although it is named so mainly for backward compatibility. Some common cases are dividing by 0 or floating point overflow.By default, the process is terminated and a core dump is created.

int main(void)

{

      /* "volatile" needed to eliminate compile-time optimizations */

      volatile int x = 42; 

      volatile int y = 0;

      x=x/y;

      return 0; /* Never reached */

}

SIGKILL(9): Kill. This signal always cause the process to terminate. It cannot be ignored or caught. It provides a sure way to kill a process for users with proper privilege.

SIGUSR1(10): User defined signal 1. It is sent to a process to indicate user-defined conditions. By default, the process is terminated.

SIGSEGV(11): Segmentation fault. Invalid memory segment access. By default, the process is terminated and a core dump is created. Some common cases of SIGSEGV are buffer overflow, using uninitialized pointers, dereferencing NULL pointers, exceeding the allowable stack size, attempting to access the memory the program doesn’t own, etc.

#include <stdio.h>

#include <stdlib.h>

int main(void) {

   int* a;

   a[1] = 0;

   return 0;

}

SIGUSR2(12): User defined signal 2. It is sent to a process to indicate user-defined conditions. By default, the process is terminated.

SIGPIPE(13): Broken pipe. By default, the process is terminated. It is sent to a process when it tries writing to a pipe without a process connecting to the other end.

SIGALRM(14): Alarm clock. By default, the process is terminated. It is sent to a process when a time limit has exceeded. Programs usually use SIGALRM to make a long running operation time out, or to perform an action at regular intervals.

SIGTERM(15): Termination signal. It is the default signal sent by kill and killall command. By default, the process is terminated.

SIGSTKFLT(16): Stack fault. By default, the process is terminated.

SIGCHLD/SIGCLD(17): Child process has stopped or exited, changed. By default, the signal is ignored by the process. A process can create a child process using fork. The signal is sent to the parent process when the child process terminates.

SIGCONT(18): Continue executing, if stopped. It resumes the process from TASK_STOPPED state and also clears any pending/queued stop signals. This happens no matter the process blocks, catches, or ignores the SIGCONT signal. By default, the signal is ignored.

SIGSTOP(19): Stop executing. The signal cannot be caught or ignored. It clears any pending/queued SIGCONT signals and stops the process.

SIGTSTP(20): Terminal stop signal. It clears any pending/queued SIGCONT signals no matter the signal is blocked, ignored or caught. Though it may or may not stop the process later on. By default, it stops the process. It is the signal sent to the process by its controlling terminal when user presses the SUSP key combination (Ctrl + Z normally).

SIGTTIN(21): Background process trying to read, from TTY. It clears any pending/queued SIGCONT signals no matter the signal is blocked, ignored or caught. Though it may or may not stop the process later on. By default, it stops the process.

SIGTTOU(22): Background process trying to write, to TTY. It clears any pending/queued SIGCONT signals no matter the signal is blocked, ignored or caught. Though it may or may not stop the process later on. By default, it stops the process.

SIGURG(23): Urgent condition on socket. By default, the signal is ignored by the process. It is sent to a process with async IO configured by fcntl system call on Linux when out-of-band data is available on a file descriptor connected to a socket.

SIGXCPU(24): CPU limit exceeded. By default, the process is terminated and a core dump is created. It is sent to a process when it has used the CPU for a duration that exceeds a certain predetermined user set value.

SIGXFSZ(25): File size limit exceeded. By default, the process is terminated and a core dump is created. It is sent to a process when it has created a file that exceeded the maximum allowed size.

SIGVTALRM(26): Virtual alarm clock. By default, the process is terminated. It is sent to a process when a time limit has reached. It counts only the time spent executing the process itself.

SIGPROF(27): Profiling alarm clock. By default, the process is terminated. It is sent to a process when a time limited has reached. It counts the time spent by the process and the system executing on behalf of the process.

SIGWINCH(28): Window size change. By default, the process is ignored. It is sent to a process when its controlling terminal size changes. It gives the process an opportunity to adjust its display.

SIGIO/SIGPOLL(29): I/O now possible. By default, the process is terminated. It is sent to a process when an async IO event occurs.

SIGPWR(30): Power failure restart. By default, the process is terminated. It is sent to a process when the system experiences power failure. This gives the process an opportunity to save its state.

SIGSYS/SIGUNUSED(31): Bad system call. By default, the process is terminated and a core dump is created. It is sent to a process when a bad argument is passed to a system call.

References:
1. Linux signal overview man page: http://linux.die.net/man/7/signal
2. Advanced Programming in the UNIX Environment, 2nd edition.
3. Linux signal.h and some other header source files.
4. SIGKILL wikipedia page: http://en.wikipedia.org/wiki/SIGILL
5. SIGFPE wikipedia page: http://en.wikipedia.org/wiki/SIGFPE
6. SIGBUS wikipedia page: http://en.wikipedia.org/wiki/SIGBUS

Programming FIFO/Named PIPE in Linux

Previous post covers pipe, an IPC mechanism for processes that have common parent process. We refer to processes that have common parent process as related process. But for unrelated processes, pipe cannot be used, because one process has no way of referring to pipes have been created by another process.

When two unrelated processes share some information, an identifier must be associated with the shared information. Therefore, one process can create the IPC object and other processes can refer to the IPC object by the identifier. Linux provides named pipe (also called FIFO) to communication through pipe in two unrelated processes.

A FIFO is created by mkfifo function,

#include <sys/types.h>

#include <sys/stat.h>

int mkfifo(const char *pathname, mode_t mode);

mkfifo creates a special file with name pathname. mode specifies the special FIFO file’s permissions. It is modified by the process’s umask: the created file will have the permission (mode & ~umask).

Once a FIFO is created, it can be operated like a usual file with the exception that both ends of FIFO need to be open first before reading and writing can begin. In other words, opening a file for reading blocks until another process open it for writing, and vice versa.

Below are two programs that uses named pipe to communicate with each other. It is modified from the pipe post example.

Firstly, the two programs need to agree on the pipe names. This is defined in a header file included by both of them.

#ifndef TEST_FIFO_H

#define TEST_FIFO_H

 

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <fcntl.h>

#include <sys/types.h>

#include <sys/stat.h>

 

#define CS_FIFO_NAME "/tmp/cs"

#define SC_FIFO_NAME "/tmp/sc"

 

//user read, user write, group read and other read

#define FIFO_MODE (S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH)

 

#endif

Next, the server program creates two pipes for communication and opens one for read and one for write,

#include "fifo.h"

 

void createfifo() {

    int rv;

    if ((rv = mkfifo(CS_FIFO_NAME, FIFO_MODE)) == -1) {

        perror("error making cs fifo: ");

        return;

    }

    if ((rv = mkfifo(SC_FIFO_NAME, FIFO_MODE)) == -1) {

        perror("error making sc fifo: ");

        return;

    }

}

 

void server(int readfd, int writefd) {

    char msg[100];

    int n;

    if ((n = read(readfd, msg, 100)) < 0) {

        perror("error reading...");

    }

    msg[n] = '';

    printf("%d server received from client: %sn", getpid(), msg);

    printf("server enter something: ");

    fgets(msg, 100, stdin);

    write(writefd, msg, strlen(msg)-1); //-1: not send the newline

}

 

int main(int argc, char **argv) {

    int readfd, writefd;

    createfifo();

    readfd = open(CS_FIFO_NAME, O_RDONLY, 0);

    writefd = open(SC_FIFO_NAME, O_WRONLY, 0);

    server(readfd, writefd);

    close(readfd);

    close(writefd);

    return 0; 

}

The client program also opens the two fifos for write and read.

#include "fifo.h"

 

void client(int readfd, int writefd) {

    char msg[100];

    int n;

    char eof = EOF;

    printf("%d client enter something: ", getpid());

    fgets(msg, 100, stdin);

    if (write(writefd, msg, strlen(msg)-1) == -1) { //-1: not send the newline

        perror("error writing...");

        exit(0);

    }

    if((n = read(readfd, msg, 100)) < 0){

        perror("error reading...");

        exit(0);

    }

    msg[n] = '';

    printf("client received from server: %sn", msg);

}

 

void removefifo() {

    unlink(SC_FIFO_NAME);

    unlink(CS_FIFO_NAME);

}

 

int main(int argc, char **argv) {

    int readfd, writefd;

 

    writefd = open(CS_FIFO_NAME, O_WRONLY, 0);

    readfd = open(SC_FIFO_NAME, O_RDONLY, 0);

    client(readfd, writefd);

    close(readfd);

    close(writefd);

    removefifo();

    return 0;

}

Note that if we swap the order of the two lines in the client source code, a deadlock will occur.

writefd = open(CS_FIFO_NAME, O_WRONLY, 0);

readfd = open(SC_FIFO_NAME, O_RDONLY, 0);

This is because the server blocks at opening CS_FIFO_NAME for reading and waits for client to open it for writing, if the client opens SC_FIFO_NAME for reading first, it also blocks and waits for server open it for writing. Both programs stuck forever.

Note that similar to pipes, it is also possible to create NONBLOCK FIFOs. Also the PIPE_BUF limits apply to FIFOs.

For a thorough explanation on these subjects, one can refer to reference 2.

References:

1. mkfifo man page: http://linux.die.net/man/3/mkfifo

2. Unix Network Programming, Volumn 2.

Programming Pipe in Linux

Pipe is one of the message passing IPC techniques. A pipe is a one way communication channel. It is created by the pipe function,

#include <unistd.h>
int pipe(int fd[2]);

The array fd[2] returns two file descriptors referring to two ends of the created pipe: fd[0] for reading and fd[1] for writing. Data written to the pipe is buffered by the kernel until it is read by the read end.

Below is a program use pipe to create a two way communication channel between two processes.

#include <stdio.h>

#include <stdlib.h>

#include <unistd.h>

#include <string.h>

 

void client(int readfd, int writefd) {

    char msg[100];

    int n;

    char eof = EOF;

    printf("%d client enter something: ", getpid());

    fgets(msg, 100, stdin);

    if (write(writefd, msg, strlen(msg)-1) == -1) { //-1: not send the newline

        perror("error writing...");

        exit(0);

    }

    if((n = read(readfd, msg, 100)) < 0){

        perror("error reading...");

        exit(0);

    }

    msg[n] = '';

    printf("client received from server: %sn", msg);

}

 

void server(int readfd, int writefd) {

    char msg[100];

    int n;

    if ((n = read(readfd, msg, 100)) < 0) {

        perror("error reading...");

    }

    msg[n] = '';

    printf("%d server received from client: %sn", getpid(), msg);

    printf("server enter something: ");

    fgets(msg, 100, stdin);

    write(writefd, msg, strlen(msg)-1); //-1: not send the newline

}

 

int main(int argc, char **argv) {

    int pipe1fd[2], pipe2fd[2];

    int pid;

    if (pipe(pipe1fd) == -1) {

        perror("pipe:");

        return 0;

    }

    if (pipe(pipe2fd) == -1) {

        perror("pipe:");

        return 0;

    }

    pid = fork();

    if (pid == 0) {

        //child process

        close(pipe1fd[1]);  //close the write end

        close(pipe2fd[0]);  //close the read end

        client(pipe1fd[0], pipe2fd[1]);

        exit(0);

    } else {

        //parent process

        close(pipe1fd[0]);

        close(pipe2fd[1]);

        server(pipe2fd[0], pipe1fd[1]);

        wait(NULL);       //wait for child process to finish

        return 0;

    }

}

Compile the code using the command,

gcc -o pipe pipe.c

And below is a screenshot of a sample run,

Figure 1. Execution of Pipe Sample Program

The program above illustrates the steps of creating two pipes and use them for two way communication.

1. Create two pipes, pipe1 (pipe1fd[2]) and pipe2 (pipe2fd[2]).
2. call fork to create a child process
3. child process close write end of pipe1 (pipe1fd[1]) and read end of pipe2 (pipe2fd[0])
4. parent process close read end of pipe1 (pipe1fd[0]) and write end of pipe2 (pipe2fd[1])

POSIX defines “half-duplex” pipes, where communication is one way in pipes. System V Release 4 (SVR4) Unix implements “full-duplex” manner, which allows the two file descriptors for both reading and writing.

GNU/Linux implements “half-duplex” in a unique manner. One does not have to close one of them in order to use the other as shown in the above program. Comment out the close statements give the program below,

#include <stdio.h>

#include <stdlib.h>

#include <unistd.h>

#include <string.h>

 

void client(int readfd, int writefd) {

    char msg[100];

    int n;

    char eof = EOF;

    printf("client enter something: ");

    fgets(msg, 100, stdin);

    if (write(writefd, msg, strlen(msg)-1) == -1) { //-1: not send the newline

        perror("error writing...");

        exit(0);

    }

    if((n = read(readfd, msg, 100)) < 0){

        perror("error reading...");

        exit(0);

    }

    msg[n] = '';

    printf("client received from server: %sn", msg);

}

 

void server(int readfd, int writefd) {

    char msg[100];

    int n;

    if ((n = read(readfd, msg, 100)) < 0) {

        perror("error reading...");

    }

    msg[n] = '';

    printf("server received from client: %sn", msg);

    printf("server enter something: ");

    fgets(msg, 100, stdin);

    write(writefd, msg, strlen(msg)-1); //-1: not send the newline

}

 

int main(int argc, char **argv) {

    int pipe1fd[2], pipe2fd[2];

    int pid;

    if (pipe(pipe1fd) == -1) {

        perror("pipe:");

        return 0;

    }

    if (pipe(pipe2fd) == -1) {

        perror("pipe:");

        return 0;

    }

    pid = fork();

    if (pid == 0) {

        //child process

//        close(pipe1fd[1]);  //close the write end

//        close(pipe2fd[0]);  //close the read end

        client(pipe1fd[0], pipe2fd[1]);

        exit(0);

    } else {

        //parent process

//        close(pipe1fd[0]);

//        close(pipe2fd[1]);

        server(pipe2fd[0], pipe1fd[1]);

        wait(NULL);       //wait for child process to finish

        return 0;

    }

}

Compile the code using,

gcc -o pipe2 pipe2.c

And below is a screenshot of a run,

Figure 2. Execution of Modified Pipe Sample Program

popen and pclose

Linux provides two functions (popen and pclose) to create pipe to or from a process. It avoids the usage of pipe, fork, close, and wait.

#include <stdio.h>
FILE *popen(const char *command, const char *type);
int pclose(FILE *stream);

popen function starts a process by creating a pipe, forking and invoking a shell. The command argument accepts a shell command line, and the command is passed to /bin/sh using -c flag. The type argument expects either “r” or “w”. Note that since a pipe is unidirectional, we cannot pass both “r” and “w”, and the returned pipe stream is either read-only or write-only.

If the returned pipe stream is read-only, the calling process reads from the standard output.

If the returned pipe stream is write-only, the calling process writes to the standard input.
pclose function waits for the associated process to terminate and returns the exit status of the command passed to popen.

Note that with popen and pclose, it is not convenient to create two pipes for two way communication. Below is a sample program using popen and pclose,

#include <stdlib.h>

#include <unistd.h>

 

void createTest() {

    FILE *f;

    f = fopen("test.txt", "w");

    fprintf(f, "hello worldn");

    fclose(f);

}

 

int main(int argc, char **argv) {

    char buf[100], command[100];

    FILE *pf;

    createTest();

    sprintf(command, "cat test.txt");

    pf = popen(command, "r");   //read from cat

    while (fgets(buf, 100, pf) != NULL) {

        printf("%s", buf);

    }

    pclose(pf);

    return 0;

}

Compile the code using the command below,

gcc -o pipe3 pipe3.c

And then run the program,

./pipe3

The program should give output as “hello world”. The code creates a new process for “cat test.txt” command, and reads from the pipe stream between cat command and itself.

PIPE_BUF

Writing less than PIPE_BUF bytes is atomic. Writing more than PIPE_BUT bytes may not be atomic: the kernel may interleave the data with data written by other processes. For example, if two processes trying to write “aaaaaa” and “bbbbb” respectively to the same pipe. If the writes are atomic, the content is either “aaaaaabbbbb” or “bbbbbaaaaaa”. But if it’s not atomic, the content can be something like “aaabbaabbba”.

In Linux, PIPE_BUF is set at limits.h, and the program below can print it out.

#include <stdio.h>

#include <limits.h>

 

int main(void) {

    printf("%dn", PIPE_BUF);

    return 0;

}

On my Ubuntu Linux, the value is 4096.

Non Block Pipe

The default pipe blocks both reads and writes. A non blocking pipe can be created by fcntl with O_NONBLOCK flag. The code block below illustrates how to create a non block pipe,

#include <fcntl.h>

 

void nonblock(int fd) {

    int flags;

    if ((flags = fcntl(fd, F_GETFL, 0)) < 0) {

        perror("F_GETFL error:");

        return;

    }

    flags |= O_NONBLOCK;

    if (fcntl(fd, F_SETFL, flags) < 0) {

        perror("F_SETFL error:");

    }

}

More details about pipes can be referred at reference 4.

References:
1. POSIX wikipedia page: http://en.wikipedia.org/wiki/POSIX
2. Understanding Linux Kernel Inter-process Communication: Pipes, FIFO & IPC: http://linux.omnipotent.net/article.php?article_id=12504&page=2
3. popen Linux man page.
4. Unix Network Programming, Volume 2.
5. http://manpages.courier-mta.org/htmlman7/pipe.7.html

Computing Fibonacci Numbers

Fibonacci series is a number sequence of 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 44, … It can be expressed by the formula below,

1

Fibonacci series grows almost as fast as the power of 2. Fn ~= 2^(0.694n).

Computation of Fibonacci Numbers
The naive approach to compute Fibonacci numbers is to recursion according to the formula of Fabonacci numbers. A C implementation is given below,

/*

a naive approach of computing fibonacci series, the computation time also satisfy fabonacci series, since

fabonacci series grows almost as fast as power of 2, the computation is expoential time

*/

#include <stdio.h>

#include <stdlib.h>

 

unsigned long fibo(int n) {

    if (n == 0) {

       return 0;

    } else if (n == 1) {

       return 1;

    } else {

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

    }

}

 

int main(int argc, char **argv) {

    int n;

    unsigned long res;

    struct timeval stTime, edTime;

    n = atoi(argv[1]);

    gettimeofday(&stTime, NULL);

    res = fibo(n);

    gettimeofday(&edTime, NULL);

    printf("%lu, time takes (%u:%u)n", res, (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec));

}

Suppose the time takes to compute fibo(n-1) and fibo(n-2) are O(n-1) and O(n-2) respectively, then the time complexity of computing fibo(n) is O(n) = O(n-1) + O(n-2). The time complexity itself follows the Fibonacci series. Since we know Fn~= 2^(0.694n), then O(n) ~=2^(0.694n). The computation is exponential. And since O(n)~=(2^0.694)^n ~=1.62^n, compute the Fn+1 takes around 1.6 times longer than computing Fn.

Compile the code using the command,

gcc -o fibo fibo.c

Below is a sample run of the program,

Figure 1. Execution of Recursive Version Program

It is not easy to observe that a Fibonacci number depends on its two previous numbers. The computation above computes some small Fibonacci numbers for many times. We can use an array to remember the previous computation. A C implementation is shown as below,

/*

use dynamic programming to remember previous computation results

this computation runs in linear time

*/

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

 

int main(int argc, char **argv) {

    int n;

    int i;

    struct timeval stTime, edTime;

    unsigned long* fibo;

    n = atoi(argv[1]);

    fibo = (unsigned long*) malloc((n+1)*sizeof(unsigned long));

    memset(fibo, 0x00, sizeof(unsigned long)*n);

    fibo[1] = 1;

    gettimeofday(&stTime, NULL);

    for (i = 2; i <= n; ++i) {

        fibo[i] = fibo[i-1] + fibo[i-2];

    }

    gettimeofday(&edTime, NULL);

    printf("%lu, time takes (%u:%u)n", fibo[n], (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec));

}

Compile the code using the command below,

gcc -o fibo1 fibo1.c

A sample run is as below,

Figure 2. Execution of Dynamic Programming Version

It is not difficult to tell the code above runs in linear time. However, the code takes a lot of memory space and if we’re to compute a Fibonacci number Fn for a very large n, then we may run out of memory space.
To compute Fn, we only need values for Fn-1 and Fn-2, this leads to the code below,

/*

use dynamic programming to remember previous computation results

this computation runs in linear time

*/

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

 

int main(int argc, char **argv) {

    int n;

    int i;

    struct timeval stTime, edTime;

    unsigned long fibo[2];

    n = atoi(argv[1]);

    fibo[0] = 0;

    fibo[1] = 1;

    gettimeofday(&stTime, NULL);

    for (i = 2; i <= n; ++i) {

        fibo[i%2] = fibo[0] + fibo[1];

    }

    gettimeofday(&edTime, NULL);

    if (n&0x01) {

        printf("%lu, time takes (%u:%u)n", fibo[1], (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec));

    } else {

        printf("%lu, time takes (%u:%u)n", fibo[0], (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec));

    }

}

Compile the code using the command below,

gcc -o fibo2 fibo2.c

A sample run is as below,

Figure 3. Execution of Dynamic Programming Version Improved

References:
1. Algorithms. Dasgupta, C.H. Papadimitriou, and U. V. Vazirani, 2006
2. Wikipedia Fibonacci Number: http://en.wikipedia.org/wiki/Fibonacci_number