problem 184 partial

gistfile1.c

// (C) Mason Larobina <[email protected]>
// http://projecteuler.net/problem=184

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>
#include <sys/param.h>

void dump(char **m, int r) {
    int w = (r * 2) + 1;

    putchar('\n');

    for (int i = 0; i < w; i++) {
        printf("%4d: ", i - r);     /* relative y */

        for (int j = 0; j < w; j++) {
            if (i == r && j == r) 
                putchar('*');       /* center point */
            else 
                putchar(m[i][j] ? '0' + m[i][j] : '.');
        }

        putchar('\n');
    }

    putchar('\n');
}

void clear(char **m, int r) {
    int w = (r * 2) + 1;
    memset(m[0], 0, w * w);
}

typedef struct point { 
    int x, y;
} point;

void paint_points(char **m, int r, point *points, int len) {
    for (int i = 0; i < len; i++) { 
        point *p = &points[i];
        m[ r + p->y ][ r + p->x ] += 1;
    }
}

// -10: ...............................
//  -9: ...........1111................
//  -8: ..........1....................
//  -7: ........11.....................
//  -6: ........1......................
//  -5: .......1.......................
//  -4: ......1........................
//  -3: ......1........................
//  -2: ......1........................
//  -1: ......1........................
//   0: ......1........*...............
//   1: ...............................
point *arc_1(int r, int n, int *len) { 
    assert(n <= r);
    assert(len != NULL);

    static point points[200];
    point *p = points;

    int r_min = r - n, 
        pn_cube = (n-1) * (n-1),
        n_cube = n * n;

    //printf("r %d n %d r_min %d\n", r, n, r_min);

    for (int i = r; i > r_min; i--) {
        int y = r - i, 
            count = 0, 
            j = sqrt((double)(pn_cube - (y * y))); 
        
        j = r - (j ? j : 1);

        for (; j > r_min; j--) {
            int x = r - j;

            int xx = x * x, 
                yy = y * y, 
                xxyy = xx + yy;

            if (xxyy < pn_cube) /* already calculated */
                continue;

            if (xxyy < n_cube) {
                count++;

                p->x = 0 - x;
                p->y = 0 - y;
                p++;

                continue; 
            } 

            break;
        }

        if (!count)
            break;
    }

    assert(p != points); 

    size_t size = sizeof(point) * (p - points);
    point *ret = malloc(size);
    memcpy(ret, points, size);
    *len = (p - points);
    return ret;
}

point *fill_2(int r, point *p1, int *len, char **m) {
    assert(len != NULL);

    static point points[20000]; 
    point *p = points;
    
    int r2 = r * r;

    double p1_angle = 0;
    if (p1->y != 0) 
        p1_angle = atan((double)p1->y / (double)p1->x);
    
    for (int i = r - 1; i > 0; i--) { 
        int y = r - i, xmin, xmax;

        xmax = xmin = (int)floor(sqrt((double)(r2 - (y * y))));

        if (p1_angle != 0)
            xmin = MIN(xmin, (int)(floor((double)y / tan(p1_angle))));

        if (-y == p1->y && -xmin == p1->x) 
            xmin--;

        xmin = r - xmin;
        xmax = r + 1 + xmax;

        for (int x = xmin; x < xmax; x++) {
            p->x = x - r;
            p->y = -y;
            p++;

            m[i][x] += 1;
        }
    }

    size_t size = sizeof(point) * (p - points);
    point *ret = malloc(size);
    memcpy(ret, points, size);
    *len = (p - points);
    return ret;
}

long fill_3(int r, point *p1, point *p2, char **m) {
    //printf("p1 (%d, %d) p2 (%d, %d)\n", p1->x, p1->y, p2->x, p2->y);

    double p1_angle = 0, p2_angle = 0;
    if (p1->y != 0) 
        p1_angle = atan((double)p1->y / (double)p1->x);

    if (p2->x != 0) 
        p2_angle = atan((double)p2->y / (double)p2->x);

    //printf("p1 %f p2 %f\n", p1_angle, p2_angle);

    int r2 = r * r;
    long total_count = 0;
    
    for (int i = r + 1; i < r*2; i++) { 
        int y = i - r, xmax, xmin, count = 0;

        xmax = xmin = (int)floor(sqrt((double)(r2 - (y * y))));

        if (p1_angle != 0)
            xmax = MIN(xmax, (int)(floor((double)y / tan(p1_angle))));

        if (p2_angle != 0)
            xmin = MIN(xmin, (int)(floor((double)y / tan(p2_angle))));

        xmin = r + MAX(0, xmin);
        xmax = r + 1 + xmax;

        //printf("xmin %d xmax %d\n", xmin, xmax);

        for (int x = xmin; x < xmax; x++) {
            //m[i][x] += 1;
            count++;
        }

        if (count) 
            total_count += count;
        else
            return total_count;
    }
}

int main(int argc, char **argv) {
    assert(argc >= 2);

    int r = atoi(argv[1]);  /* max Ir */

    printf("Ir = %d\n", r);
    
    // pointer buffer size[200]
    assert(r <= 120);

    int w = (r * 2) + 1;
    char **m = malloc(sizeof(char*) * w); /* pointers to sub arrays */
    m[0] = malloc(sizeof(char) * w * w);  /* data chunk */

    clear(m, r);

    for (int i = 1; i < w; i ++)
        m[i] = m[0] + i * w;        /* set pointers to sub arrays */
    
    int n = argc >= 3 ? atoi(argv[2]) : r;
     
    int p_len;
    point *points1 = arc_1(r, n, &p_len);

    long count = 0;
    
    for (int i = 0; i < p_len; i++) { 
        point *p1 = &points1[i];
        paint_points(m, r, p1, 1);
        
        int p2_len;
        point *points2 = fill_2(r, p1, &p2_len, m);
        //printf("len2 %d\n", p2_len);

        for (int j = 0; j < p2_len; j++) {
            point *p2 = &points2[j];

     //       clear(m, r);
            paint_points(m, r, p1, 1);
            paint_points(m, r, p2, 1);
            count += fill_3(r, p1, p2, m);
    //        dump(m, r);
        }
    }

    printf("count %ld\n", count * 4);

    return 0;
}