Lagrange polynomial and coefficient calculations

lagint.c

/*
 *  lagint.c
 *  2019-2020  Sun Dro (f4tb0y@protonmail.com)
 *
 * Calculate Lagrange interpolated polynomials and 
 * displays the coefficient of each polynomials
 */

#include <stdio.h>
#include <errno.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
#include <signal.h>
#include <fcntl.h>
#include <math.h>

#include <sys/types.h>
#include <sys/wait.h>
#include <sys/stat.h>
#include <sys/file.h>

/*
    Example input:

    6.0,8.0,1.0,12.0,5.0,5.0,9.0,7.0,8.0,2.0,7.0,4.0,4.0,6.0,5.5,2.0
    4.0,9.0,6.0,3.0,2.0,9.0,7.0,3.0,9.0,2.0,1.0,6.0,8.0,7.0,5.8,4.0
    1.0,8.0,7.0,7.0,5.0,9.0,2.0,3.0,9.0,1.0,17.0,0.0,3.0,3.0,8.0,5.0
    9.0,10.0,7.0,7.0,8.0,9.0,5.0,3.0,2.0,1.0,3.0,1.0,4.0,8.0,6.0,7.0
    1.0,1.0,3.0,4.0,5.0,8.0,7.0,3.0,15.0,1.0,4.0,4.0,8.0,10.0,9.0,7.0
    8.0,7.0,1.0,0.0,18.0,8.0,2.0,9.0,4.0,4.0,3.0,6.0,6.0,8.0,10.0,9.0
    5.0,8.0,8.0,3.0,7.0,7.0,3.0,2.0,15.0,1.0,10.0,6.0,6.0,6.0,4.1,6.0
    1.0,9.0,3.0,3.0,4.0,4.0,9.0,4.0,6.0,2.0,7.0,1.0,5.0,3.0,3.5,3.0
*/

#define BUFF_MAX 8192
#define SIZE_ARR 32

struct Node {
    double coeff;
    int power;
    struct Node* next;
};

struct polynomial {
    struct Node *p[10];
};

void exit_failure()
{
    /* Terminate every process before exit */
    killpg(getpgrp(), SIGUSR2);
    exit(EXIT_FAILURE);
}

int read_file(int fd, char *buffer, size_t size)
{
    lseek(fd, SEEK_SET, 0);

    int bytes = read(fd, buffer, size - 1);
    int length = (bytes > 0) ? bytes : 0;

    if (length == 0 && errno == 0)
        return read_file(fd, buffer, size);

    buffer[length] = '\0';
    return length;
}

int get_line_count(char *buffer, size_t length)
{
    size_t i, lastp = 0;
    int count = 0;

    /* Count new line characters */
    for (i = 0; i < length; i++)
    {
        if (buffer[i] == '\n')
        {
            lastp = i;
            count++;
        }
    }

    /* Check if last line doesnot contain new line character */
    if ((i - lastp > 1) && buffer[i] != '\n') count++;

    return count;
}

int get_line(const char *path, int num, char *output, size_t size)
{
    int fd = open(path, O_RDONLY, 644);
    if (fd < 0)
    {
        fprintf(stderr, "Can not open file: %s\n", strerror(errno));
        return 0;
    }

    char buffer[BUFF_MAX];
    struct flock lock;

    lock.l_type = F_RDLCK;
    lock.l_whence = SEEK_SET;
    lock.l_start = 0;
    lock.l_len = 0;

    if (fcntl(fd, F_SETLKW, &lock) == -1) 
    {
        fprintf(stderr, "Can not lock the file: %s\n", strerror(errno));
        exit_failure();
    }

    int bytes = read_file(fd, buffer, sizeof(buffer));
    if (bytes <= 0)
    {
        fprintf(stderr, "Can not read file: %s\n", strerror(errno));
        return 0;
    }

    /* Also releases all locks */
    close(fd);

    int i, len = 0, line_num = 1;

    for (i = 0; i < bytes; i++)
    {
        if (num == line_num)
        {
            output[len++] = buffer[i];
            if (len == size) break;
        }

        if (buffer[i] == '\n') line_num++;
    }

    output[len] = '\0';
    return len;
}

int put_line(const char *path, int num, const char *line, size_t length)
{
    int fd = open(path, O_RDWR, 644);
    if (fd < 0)
    {
        fprintf(stderr, "Can not open file: %s\n", strerror(errno));
        return 1;
    }

    char buffer[BUFF_MAX];
    struct flock lock;

    lock.l_type = F_WRLCK;
    lock.l_whence = SEEK_SET;
    lock.l_start = 0;
    lock.l_len = 0;

    if (fcntl(fd, F_SETLKW, &lock) == -1) 
    {
        fprintf(stderr, "Can not lock the file: %s\n", strerror(errno));
        exit_failure();
    }

    int bytes = read_file(fd, buffer, sizeof(buffer));
    if (bytes <= 0)
    {
        fprintf(stderr, "Can not read file: %s\n", strerror(errno));
        return 0;
    }

    char front_buffer[BUFF_MAX];
    char back_buffer[BUFF_MAX];

    int front_len = 0;
    int front_pos = 0;
    int back_len = 0;
    int i, line_num = 1;

    for (i = 0; i < bytes; i++)
    {
        if (line_num < num)
        {
            front_buffer[front_len++] = buffer[i];
            front_pos = i;
        }
        else if (line_num > num)
        {
            back_buffer[back_len++] = buffer[i];
        }

        if (buffer[i] == '\n') line_num++;
    }

    int full_len = front_len + back_len + length;
    int front_offset = 0, back_offset = 0, line_offset = 0;
    int insert_done = 0;

    for (i = 0; i < full_len; i++)
    {
        if (front_pos && i <= front_pos)
        {
            buffer[i] = front_buffer[front_offset++];
        }
        else if (!insert_done)
        {
            buffer[i] = line[line_offset++];
            if (line_offset == length) insert_done = 1;
        }
        else
        {
            buffer[i] = back_buffer[back_offset++];
        }
    }

    buffer[full_len] = '\0';

    /* 
        It is necessary to truncate the file manually because the O_TRUNC flag 
        will damage the result of read() function and due to the synchronization 
        the reading and writing must be done in one call of the open() function.
    */
    ftruncate(fd, 0);

    /* Move offset in the begining of the file */
    lseek(fd, SEEK_SET, 0);

    bytes = write(fd, buffer, full_len);
    close(fd); /* Also releases all locks */

    return bytes;
}

int get_values(const char *line, int len, float (*x)[SIZE_ARR], float (*y)[SIZE_ARR])
{
    char value[16];
    int i, offset = 0, isx = 0;

    int ix = 0, iy = 0;

    for (i = 0; i < len; i++)
    {
        value[offset] = line[i];

        if (value[offset] == ',' || value[offset] == '\n')
        {
            value[offset] = '\0';
            offset = 0;
            isx = !isx;

            if (isx) (*x)[ix++] = atof(value);
            else (*y)[iy++] = atof(value);

            if (ix >= SIZE_ARR || iy >= SIZE_ARR)
            {
                fprintf(stderr, "Too big input\n");
                return 0;
            }

            continue;
        }

        offset++;
    }

    return ix;
}
 
/* Function add a new node at the end of list */
struct Node* add_node(struct Node* start, double coeff, int power)
{
    /* Create a new node */
    struct Node* newnode = (struct Node*)malloc(sizeof(struct Node)); /* new Node; */
    newnode->coeff = coeff;
    newnode->power = power;
    newnode->next = NULL;
 
    /* If linked list is empty */
    if (start == NULL)
        return newnode;
 
    /* If linked list has nodes */
    struct Node* ptr = start;
    while (ptr->next != NULL)
        ptr = ptr->next;
    ptr->next = newnode;
 
    return start;
}
 
/* Functionn to display coefficients from linked list */
void print_coeff(struct Node* ptr, int id)
{
    printf("Polynomial %d: ", id);

    while (ptr->next != NULL) 
    {
        printf("%s%0.1f,", 
            (ptr->coeff >= 0) ? 
            "+" : "", ptr->coeff);

        ptr = ptr->next;
    }

    printf("%0.1f\n", ptr->coeff);
}

void clear_list(struct Node* head)
{
    while (head != NULL)
    { 
       struct Node *tmp = head;
       head = head->next;
       free(tmp);
    }
}

/* Function to add coefficients of two elements having same powerer */
void remove_duplicates(struct Node* start)
{
    struct Node *ptr1; 
    struct Node *ptr2;
    struct Node *dup;
    ptr1 = start;
 
    /* Pick elements one by one */
    while (ptr1 != NULL && ptr1->next != NULL) 
    {
        ptr2 = ptr1;
 
        /* 
            Compare the picked element
            with rest of the elements 
        */
        while (ptr2->next != NULL) 
        {
 
            /* If powerer of two elements are same */
            if (ptr1->power == ptr2->next->power) 
            {
 
                /* Add their coefficients and put it in 1st element */
                ptr1->coeff = ptr1->coeff + ptr2->next->coeff;
                dup = ptr2->next;
                ptr2->next = ptr2->next->next;
 
                /* remove the 2nd element */
                free(dup);
            }
            else
                ptr2 = ptr2->next;
        }

        ptr1 = ptr1->next;
    }
}
 
/* Function to multiply two polynomial numbers */
struct Node* multiply(struct Node* poly1,struct Node* poly2,
               struct Node* poly3)
{
    /*
        Create two pointer and store the
        address of 1st and 2nd polynomials 
    */
    struct Node *ptr1 = poly1;
    struct Node *ptr2 = poly2;

    while (ptr1 != NULL) 
    {
        while (ptr2 != NULL) 
        {
            double coeff;
            int power;
 
            /* 
                Multiply the coefficient of both
                polynomials and store it in coeff
            */
            coeff = ptr1->coeff * ptr2->coeff;         
            power = ptr1->power + ptr2->power;

            /* 
                Invoke add_node function to create a 
                newnode by passing three parameters 
            */
            poly3 = add_node(poly3, coeff, power);
 
            /* 
                move the pointer of 2nd polynomial
                two get its next term
            */
            ptr2 = ptr2->next;
        }
 
        /* 
            Move the 2nd pointer to the
            starting point of 2nd polynomial 
        */
        ptr2 = poly2;
 
        /* move the pointer of 1st polynomial */
        ptr1 = ptr1->next;
    }
 
    /* 
        this function will be invoke to add
        the coefficient of the elements
        having same powerer from the resultant linked list 
    */
    remove_duplicates(poly3);
    return poly3;
}

/*
    This is the function to Calculate Lagrange and displays the coefficient of lagrange
*/
struct Node* func(struct Node* poly1, struct Node* poly2, double val)
{	
    /* 
        Create two pointer and store the
        address of 1st and 2nd polynomials 
    */
    struct Node *ptr1 = poly1;

    while (ptr1 != NULL) 
    {
        double coeff;
        int power;

        coeff = ptr1->coeff * val;         
        power = ptr1->power;

        poly2 = add_node(poly2, coeff, power);
        ptr1 = ptr1->next;
    }

    return poly2;
}

struct Node* func1(struct Node* poly1, struct Node* poly2)
{	
    /* 
        Create two pointer and store the 
        address of 1st and 2nd polynomials 
    */
    struct Node *ptr1 = poly1;

    while (ptr1 != NULL) 
    {
        double coeff;
        int power;

        coeff = ptr1->coeff;         
        power = ptr1->power;

        poly2 = add_node(poly2, coeff, power);
        ptr1 = ptr1->next;
    }

    return poly2;
}

void multiplay_polynomial(struct polynomial *poly[], int count, double* arr, int id)
{
    struct Node *h[count];
    struct Node *h1[count];
    int i, val = 0, counter = 0;

    for (i = 0; i < count; i++)
    {
        h[i] = NULL;
        h1[i] = NULL;
    }

    for (i = 0; i < count; i++)
    {
        int j;

        for (j = 0; j < count-2; j++)
        {
            struct Node *hold = NULL;
            hold = multiply(poly[i]->p[j], poly[i]->p[j+1], hold);

            clear_list(poly[i]->p[j+1]);
            poly[i]->p[j+1] = hold;
        }

        h[i] = poly[i]->p[j];
        val++;
    }

    for (i = 0; i < val; i++) h1[i] = func(h[i], h1[i], arr[i]);
    for (i = 1; i < val; i++) h1[0] = func1(h1[i],h1[0]);

    remove_duplicates(h1[0]);
    print_coeff(h1[0], id);

    for (i = 0; i < count; i++)
        clear_list(h1[i]);
}

double calculate_lagrange(float x[SIZE_ARR], float y[SIZE_ARR], int count, double xp, int display_id)
{
    /* Store the Result */
    double result_sum = 0;

    /* Array to store value of Lagrange coefficients */
    double *coeff_array = (double*)malloc(sizeof(double)*count);
    double *denominator = (double*)malloc(sizeof(double)*count);

    struct polynomial *poly[count];
    int i, j, counter[count];

    for (i = 0; i < count; i++)
        poly[i] = (struct polynomial*)malloc(sizeof(struct polynomial));

    for (i = 0; i < count; i++)
    {
        counter[i] = 0;
        double result = 1;
        double c1 = 1.0;

        for (j = 0; j < count; j++)
        {
            /* Using Formula */
            if (j != i)
            {
                struct Node *poly1 = NULL;
                poly1 = add_node(poly1, 1, 1);
                poly1 = add_node(poly1, -x[j], 0);
                poly[i]->p[counter[i]++] = poly1;

                result = result * (double)((xp - x[j]) / (x[i] - x[j]));
                c1 *= (x[i] - x[j]);
            }
        }

        denominator[i] = y[i] / c1;
        coeff_array[i] = result;
        result_sum += y[i] * coeff_array[i];
    }

    if (display_id >= 0)
        multiplay_polynomial(poly, count, denominator, display_id);

    free(coeff_array);
    free(denominator);

    for (i = 0; i < count; i++)
    {
        for (j = 0; j < counter[i]; j++)
            clear_list(poly[i]->p[j]);
    
        free(poly[i]);
    }

    return result_sum;
}

float estimate_error(float y_real, float y_estimated) {
    return (y_real-y_estimated) / y_real * 100;
}

void signal_handler(int sig)
{
    if (sig == SIGINT)
    {
        printf("Interrupt signal: %d (CTRL+C)\n", sig);
        killpg(getpgrp(), SIGINT); /* Interrupt every child before exit */
        exit(EXIT_SUCCESS); /* Interrupt is not an error */
    }
    else if (sig == SIGUSR2) 
    {
        printf("Terminated by process: %d\n", sig);
        exit(EXIT_FAILURE);
    }
    else if (sig == SIGUSR1)
    {
        /*printf("Received SIGUSR1: %d\n", getpid()); /* debug */

        /* 
            Do nothing, Signal handler for SIGUSR1 is registered only to avoid 
            the default action for this signal which will terminate the program.
        */
    }
}

int register_signals()
{
    struct sigaction sact;

    sigemptyset(&sact.sa_mask);
    sact.sa_flags = 0;
    sact.sa_handler = signal_handler;

    if (sigaction(SIGUSR1, &sact, NULL) != 0)
    {
        fprintf(stderr, "Failed to setup SIGUSR1 handler: %s\n", strerror(errno));
        return 0;
    }

    if (sigaction(SIGUSR2, &sact, NULL) != 0)
    {
        fprintf(stderr, "Failed to setup SIGUSR2 handler: %s\n", strerror(errno));
        return 0;
    }

    if (sigaction(SIGINT, &sact, NULL) != 0)
    {
        fprintf(stderr, "Failed to setup SIGINT handler: %s\n", strerror(errno));
        return 0;
    }

    /*
        Explicitly setting the disposition of SIGCHLD to SIG_IGN causes 
        any child process that subsequently terminates to be immediately 
        removed from the system instead of being converted into a zombie.
    */
    sact.sa_handler = SIG_IGN;

    if (sigaction(SIGCHLD, &sact, NULL) != 0)
    {
        fprintf(stderr, "Failed to setup SIGCHLD handler: %s\n", strerror(errno));
        return 0;
    }

    return 1;
}

void mask_and_suspend()
{
    sigset_t mask, oldmask;
    sigemptyset(&mask);
    sigaddset(&mask, SIGUSR1);

    sigprocmask(SIG_BLOCK, &mask, &oldmask);
    sigsuspend(&oldmask);
    sigprocmask(SIG_UNBLOCK, &mask, NULL);
}

void notify_process(pid_t pid)
{
    if (kill(pid, SIGUSR1) == -1)
    {
        fprintf(stderr, "Failed to send signal: %s\n", strerror(errno));
        exit_failure();
    }
}

void start_calculations(const char *path, int num)
{
    /* Wait access to start calculation */
    mask_and_suspend();

    char line[BUFF_MAX];
    float x[SIZE_ARR], y[SIZE_ARR];

    int length = get_line(path, num, line, sizeof(line));
    if (length <= 0) return;

    int count = get_values(line, length, &x, &y);
    if (count == 0)
    {
        fprintf(stderr, "Invalid line: %d\n", num);
        return;
    }

    /* First round */
    float p = calculate_lagrange(x, y, 6, x[6], -1);

    snprintf(line, sizeof(line), 
        "%.1f,%.1f,%.1f,%.1f,%.1f,%.1f,%.1f,%.1f,"
        "%.1f,%.1f,%.1f,%.1f,%.1f,%.1f,%.1f,%.1f,%.1f\n",
        x[0],y[0],x[1],y[1],x[2],y[2],x[3],y[3],
        x[4],y[4],x[5],y[5],x[6],y[6],x[7],y[7],p);

    if (!put_line(path, num, line, strlen(line)))
    {
        fprintf(stderr, "Can not put line: %d (%s)", num, strerror(errno));
        return;
    }

    /* Notify parent process about finish of first round */
    notify_process(getppid());

    /* Wait parent to accept second round */
    mask_and_suspend();

    /* Second round (also print polynomial’s 7 coefficients) */
    float r = calculate_lagrange(x, y, 7, x[6], num-1);

    snprintf(line, sizeof(line), 
        "%.1f,%.1f,%.1f,%.1f,%.1f,%.1f,%.1f,%.1f,%.1f,"
        "%.1f,%.1f,%.1f,%.1f,%.1f,%.1f,%.1f,%.1f,%.1f\n",
        x[0],y[0],x[1],y[1],x[2],y[2],x[3],y[3],
        x[4],y[4],x[5],y[5],x[6],y[6],x[7],y[7],p,r);

    if (!put_line(path, num, line, strlen(line))) 
        fprintf(stderr, "Can not put line: %d (%s)", num, strerror(errno));
}

pid_t spawn_process(const char *path, int num)
{
    /*
        On success, the PID of the child process is returned in the parent, 
        and 0 is returned in the child. On failure, -1 is returned in the 
        parent, no child process is created and errno is set appropriately. 
    */
    pid_t pid = fork();

    if (pid == -1)
    {
        fprintf(stderr, "Failed to spawn child process: %d (%s)\n", 
            num, strerror(errno));

        exit_failure();
    }
    else if (pid == 0)
    {
        /* We are in a child process */
        start_calculations(path, num);
        exit(EXIT_SUCCESS);
    }

    return pid;
}

int main(int argc, char* argv[]) 
{
    if (argc < 2)
    {
        fprintf(stderr, "Please specify file name\n");
        return 1;
    }

    const char *path = argv[1];
    char buffer[BUFF_MAX];

    int fd = open(path, O_RDWR, 644);
    if (fd < 0)
    {
        fprintf(stderr, "Can not open file: %s\n", strerror(errno));
        return 1;
    }

    int length = read_file(fd, buffer, sizeof(buffer));
    if (length <= 0)
    {
        fprintf(stderr, "Can not read file: %s\n", strerror(errno));
        return 1;
    }

    close(fd);

    int count = get_line_count(buffer, length);
    if (count == 0)
    {
        fprintf(stderr, "Invalid or empty file: %s\n", path);
        return 1;
    }

    if (!register_signals())
    {
        fprintf(stderr, "Failed to register signal handlers: %s\n", path);
        return 1;
    }

    pid_t pids[count];
    int i;

    /* Use only first 8 row */
    count = count > 8 ? 8 : count;

    /* Spawn all child processes */
    for (i = 1; i <= count; i++)
        pids[i-1] = spawn_process(path, i);

    /* Wait for first round to be completed */
    for (i = 0; i < count; i++)
    {
        notify_process(pids[i]);
        mask_and_suspend();
    }

    int used_count = 0;
    float err_sum = 0.;

    for (i = 1; i <= count; i++)
    {
        char line[BUFF_MAX];
        float x[SIZE_ARR], y[SIZE_ARR];

        int length = get_line(path, i, line, sizeof(line));
        if (length <= 0) return 1;

        int x_count = get_values(line, length, &x, &y);
        if (x_count == 0)
        {
            fprintf(stderr, "Invalid line: %d\n", i);
            return 1;
        }

        if (x_count > 8)
        {
            /* x[8] is calculated p(x_7) with degree of 5 */
            err_sum += estimate_error(y[6], x[8]);
            used_count++;
        }
    }

    float average_err = (err_sum / used_count);
    average_err = (average_err < 0) ? -average_err : average_err;
    printf("Error of polynomial of degree 5: %.1f\n", average_err);

    /* Run second round in childs */
    for (i = 0; i < count; i++) notify_process(pids[i]);

    /* Wait for second round to be completed */
    for (i = 0; i < count; i++) wait(NULL);

    used_count = 0;
    err_sum = 0.;

    for (i = 1; i <= count; i++)
    {
        char line[BUFF_MAX];
        float x[SIZE_ARR], y[SIZE_ARR];

        int length = get_line(path, i, line, sizeof(line));
        if (length <= 0) return 1;

        int x_count = get_values(line, length, &x, &y);
        if (x_count == 0)
        {
            fprintf(stderr, "Invalid line: %d\n", i);
            return 1;
        }

        if (x_count > 8)
        {
            /* Actully y[8] is calculated p(x_7) with degree of 6 */
            err_sum += estimate_error(y[6], y[8]);
            used_count++;
        }
    }

    average_err = (err_sum / (float)used_count);
    average_err = (average_err < 0) ? -average_err : average_err;
    printf("Error of polynomial of degree 6: %.1f\n", average_err);

    /* That's all */
    return 0;
}