psat/cnf.c

#include "cnf.h"

cnf* readDIMACS(char* path) {
    cnf* c = malloc(sizeof(cnf));
    CHECK(c, "Failed to alloc CNF struct\n");

    FILE* f = fopen(path, "r");
    CHECK(f, "Failed to open file\n");

    u32 bufsize = 1000;
    char* buf = malloc(sizeof(char) * bufsize);
    CHECK(buf, "Failed to alloc read buffer\n");

    CHECK(fgets(buf, sizeof(char) * bufsize, f), "Failed to read file\n");

    // 'c' marks the beginning of a comment line in DIMACs. Skip this line
    while (buf[0] == 'c') {
        CHECK(fgets(buf, sizeof(char) * bufsize, f), "Failed to read file\n");
    }

    // The first non-comment line is the header, should have format 'p cnf [varcnt] [clausecnt]
    // We skip the first bit by iterating and skipping any character that's p, c, n, f, or space
    char* temp = buf;
    while (*temp == ' ' || *temp == 'p' || *temp == 'c' || *temp == 'n' || *temp == 'f') {
        temp++;
    }

    // Char by char we read in the number of variables
    c->varcnt = 0;
    while (((u8) (*temp - '0')) < 10) {
        c->varcnt *= 10;
        c->varcnt += (*temp - '0');
        temp++;
    }

    // Skip any trailing whitespace
    while (*temp == ' ') temp++;

    // Read in clausecnt
    c->clausecnt = 0;
    while (((u8) (*temp - '0')) < 10) {
        c->clausecnt *= 10;
        c->clausecnt += (*temp - '0');
        temp++;
    }

    // Resize line buffer: the maximum clause size is varcnt, the maximum size of a u32 in decimal
    // is 10, and conservatively given 2 extra characters for the - symbol and space delimiting,
    // resizing the buffer to this should mean we never encounter a clause we fail to read.
    u32 nsize = sizeof(char) * (c->varcnt * 12LU) + 10;
    char* nbuf = realloc(buf, nsize);
    CHECK(nbuf, "Failed to realloc read buffer\n");

    buf = nbuf;

    u32 cnt = 0;
    u32 cap = 32;

    // We store separate arrays: literal i of the CNF has a parity at index i in one array
    // the variable it corresponds to at the same index in another
    // And the clause it's a member of in a third
    c->variables = calloc( cap, sizeof(u32));
    CHECK(c->variables, "Failed to allocate literal variables\n");
    c->clauses = calloc( cap, sizeof(u32));
    CHECK(c->clauses, "Failed to allocate literal clauses\n");
    c->pars = calloc( cap, sizeof(u8));
    CHECK(c->pars, "Failed to allocate literal parities\n");
    c->lastvars = calloc( c->clausecnt, sizeof(u32));
    CHECK(c->lastvars, "Failed to allocate clause lastvars\n");


    for (u32 i = 0; i < c->clausecnt; ++i) {
        CHECK(fgets(buf, nsize, f), "Failed to read clause\n");

        temp = buf;
        while (*temp == ' ') temp++;

        // Iterate through char by char
        bool empty = true;
        bool tr = true;
        while (*temp != '\n') {
            if (cnt == cap) {
                // Out of space to add more literals, realloc arrays
                u32 ncap = cap << 1U;
                c->variables = realloc(c->variables, sizeof(u32) * ncap);
                memset(c->variables + cap, 0, sizeof(u32) * cap);
                CHECK(c->variables, "Failed to realloc variable array\n");
                c->clauses = realloc(c->clauses, sizeof(u32) * ncap);
                memset(c->clauses + cap, 0, sizeof(u32) * cap);
                CHECK(c->clauses, "Failed to realloc clause array\n");
                c->pars = realloc(c->pars, sizeof(u8) * ncap);
                memset(c->pars + cap, 0, sizeof(u8) * cap);
                CHECK(c->pars, "Failed to realloc parity array\n");
                cap = ncap;
            }

            if (*temp == '-') {
                // Mark that the literal currently being read is negated
                tr = false;
            } else if (((u8) (*temp - '0'))  < 10) {
                // Read in the literal's digits
                c->variables[cnt] *= 10;
                c->variables[cnt] += (*temp - '0');
            } else {
                // Skip any whitespace and pack the read value into the arrays, along with
                // any additional data
                while (temp[1] == ' ') temp++;
                c->pars[cnt] = tr;
                c->clauses[cnt] = i;
                c->variables[cnt] -= 1;
                if (c->lastvars[i] < c->variables[cnt]) c->lastvars[i] = c->variables[cnt];
                empty = false;
                tr = true;
                cnt++;
                if (temp[1] == '0') {
                    break;
                }
            }
            temp++;
        }

        if (empty) {
            printf("UNSAT: Empty clause %u\n", i);
            return NULL;
        }
    }

    // Lastvars is set to the index of the last variable of each clause. However, because in the
    // counter, the last variable has the lsb position, we flip the value of lastvars
    for (u32 i = 0; i < c->clausecnt; ++i) {
        c->lastvars[i] += 1U;
        c->lastvars[i] = c->varcnt - c->lastvars[i];
    }

    c->litcnt = cnt;

    // Realloc the arrays to exactly match the number of literals
    c->variables = realloc(c->variables, sizeof(u32) * c->litcnt);
    c->clauses = realloc(c->clauses, sizeof(u32) * c->litcnt);
    c->pars = realloc(c->pars, sizeof(u8) * c->litcnt);
    free(buf);
    if (fclose(f)) {
        printf("Failed to close file\n");
        return NULL;
    }
    return c;
}

void printcnf(cnf* c) {
    printf("p cnf %u %u\n", c->varcnt, c->clausecnt);
    u32 pclause = 0;
    for (u32 i = 0; i < c->litcnt; ++i) {
        if (c->clauses[i] != pclause) printf("0\n");
        if (c->pars[i] == 0) printf("-");
        printf("%u ", c->variables[i] + 1);
        pclause = c->clauses[i];
    }
    printf("0\n");
}

void freecnf(cnf* c) {
    free(c->pars);
    free(c->clauses);
    free(c->variables);
    free(c->lastvars);
    free(c);
}

const u64 RBITS = 8;
const u64 RSIZE = 1LU << RBITS;
const u64 RLVLS = (63LU / RBITS) + 1;
const u64 RMASK = RSIZE - 1;

void sortlastnum(cnf* c, u64 N) {
    // Radix sort on the literals themselves based on their clause's
    bool v = false;
    u32* d = malloc(sizeof(u32) * c->litcnt);
    u32* d2 = malloc(sizeof(u32) * c->litcnt);
    u8* d3 = malloc(sizeof(u8) * c->litcnt);

    u64 qb[RLVLS * RSIZE];
    u64 qe[RLVLS * RSIZE];
    memset(qb, 0, sizeof(u64) * RLVLS * RSIZE);
    memset(qe, 0, sizeof(u64) * RLVLS * RSIZE);

    for (u64 pass = 0; pass < RLVLS; ++pass) {
        u64 shift = pass * RBITS;

        for (u64 i = 0; i < N; ++i) {
            // U32MAX - because we want to sort highest to lowest, not lowest to highest
            u32 val = UINT32_MAX - c->lastvars[c->clauses[i]];
            u64 ind = (val >> shift) & RMASK;
            qe[pass * RSIZE + ind]++;
        }

        u64 uc = 0;
        for (u64 i = 0; i < RSIZE; ++i) if (qe[pass * RSIZE + i]) uc++;
        if (uc == 1) continue;

        qb[pass * RSIZE] = 0;
        for (u64 i = 1; i < RSIZE; ++i) {
            qb[pass * RSIZE + i] = qb[pass * RSIZE + i - 1] + qe[pass * RSIZE + i - 1];
        }

        for (u64 i = 0; i < N; ++i) {
            u32 val = UINT32_MAX - c->lastvars[c->clauses[i]];
            u64 ind = (val >> shift) & RMASK;
            d[qb[pass * RSIZE + ind]] = c->variables[i];
            d2[qb[pass * RSIZE + ind]] = c->clauses[i];
            d3[qb[pass * RSIZE + ind]] = c->pars[i];
            qb[pass * RSIZE + ind]++;
            __builtin_prefetch(d + qb[pass * RSIZE + ind] + 1);
            __builtin_prefetch(d2 + qb[pass * RSIZE + ind] + 1);
            __builtin_prefetch(d3 + qb[pass * RSIZE + ind] + 1);
        }

        // Every iteration we swap pointers to avoid doing a memcpy
        u32* tptr = c->variables;
        c->variables = d;
        d = tptr;
        tptr = c->clauses;
        c->clauses = d2;
        d2 = tptr;
        u8* tptr2 = c->pars;
        c->pars = d3;
        d3 = tptr2;
        v = !v;
    }

    // If the pointers are still swapped at the end, we swap them to their original locations and
    // copy over the result
    if (v) {
        u32* tptr = c->variables;
        c->variables = d;
        d = tptr;
        tptr = c->clauses;
        c->clauses = d2;
        d2 = tptr;
        u8* tptr2 = c->pars;
        c->pars = d3;
        d3 = tptr2;
        memcpy(c->variables, d, N * sizeof(u32));
        memcpy(c->clauses, d2, N * sizeof(u32));
        memcpy(c->pars, d3, N * sizeof(u8));
    }
    free(d);
    free(d2);
    free(d3);

    // Initialize an array with all the clause indices in their current order
    u32* swaparr = malloc(sizeof(u32) * c->clausecnt);
    for (u32 i = 0; i < c->clausecnt; ++i) {
        swaparr[i] = i;
    }

    // Perform another radix sort!
    v = false;
    d = malloc(sizeof(u32) * c->clausecnt);

    memset(qb, 0, sizeof(u64) * RLVLS * RSIZE);
    memset(qe, 0, sizeof(u64) * RLVLS * RSIZE);

    for (u64 pass = 0; pass < RLVLS; ++pass) {
        u64 shift = pass * RBITS;

        for (u64 i = 0; i < c->clausecnt; ++i) {
            // Sort highest to lowest by lastvar again, but this time sorting the indices in swaparr
            // instead of sorting literals
            u32 val = UINT32_MAX - c->lastvars[swaparr[i]];
            u64 ind = (val >> shift) & RMASK;
            qe[pass * RSIZE + ind]++;
        }

        u64 uc = 0;
        for (u64 i = 0; i < RSIZE; ++i) if (qe[pass * RSIZE + i]) uc++;
        if (uc == 1) continue;

        qb[pass * RSIZE] = 0;
        for (u64 i = 1; i < RSIZE; ++i) {
            qb[pass * RSIZE + i] = qb[pass * RSIZE + i - 1] + qe[pass * RSIZE + i - 1];
        }

        for (u64 i = 0; i < c->clausecnt; ++i) {
            u32 val = UINT32_MAX - c->lastvars[swaparr[i]];
            u64 ind = (val >> shift) & RMASK;
            d[qb[pass * RSIZE + ind]] = swaparr[i];
            qb[pass * RSIZE + ind]++;
            __builtin_prefetch(d + qb[pass * RSIZE + ind] + 1);
        }

        u32* tptr = swaparr;
        swaparr = d;
        d = tptr;
        v = !v;
    }

    if (v) {
        u32* tptr = swaparr;
        swaparr = d;
        d = tptr;
        memcpy(swaparr, d, c->clausecnt * sizeof(u32));
    }
    //free(d);

    // Swaparr now has each clause index value at the correct index for where
    // it appears in the CNF now that it has been sorted
    // We'll now invert the permutation in swaparr
    for (uint i = 0; i < c->clausecnt; ++i) {
        d[swaparr[i]] = i;
    }

    // giving us a lookup table from oldindex to newindex
    // Which we then correct in all the literals
    for (uint i = 0; i < c->litcnt; ++i) {
        c->clauses[i] = d[c->clauses[i]];
    }

    // and then finally reorder the lastvars
    for (uint i = 0; i < c->clausecnt; ++i) {
        swaparr[d[i]] = c->lastvars[i];
    }

    // Copy the lastvars back into the correct array, clean up, and we're done!
    memcpy(c->lastvars, swaparr, sizeof(u32) * c->clausecnt);
    free(swaparr);
    free(d);
}