Progress on early subsumption

Signed-off-by: gothictomato <gothictomato@pm.me>

CMakeLists.txt

@@ -2,8 +2,8 @@ cmake_minimum_required(VERSION 3.22)
 project(psat C)
 
 set(CMAKE_C_STANDARD 99)
-# set(CMAKE_C_FLAGS "-O3")
+# set(CMAKE_C_FLAGS "-mavx2 -O3 -ftree-loop-linear -ftree-loop-im -ftree-loop-ivcanon -fivopts -ftree-vectorize -ftracer -funroll-all-loops ")
 
-add_executable(psat main.c cnf.c cnf.h time.c time.h types.h gpusolver.c gpusolver.h cpusolver.c cpusolver.h tests/masterTest.c tests/masterTest.h)
+add_executable(psat main.c cnf.c cnf.h time.c time.h types.h gpusolver.c gpusolver.h cpusolver.c cpusolver.h tests/masterTest.c tests/masterTest.h ncnf.c ncnf.h rng.h rng.c)
 
-target_link_libraries(psat -lOpenCL)
+target_link_libraries(psat -lOpenCL -lgmp)

cnf.c

@@ -1,5 +1,6 @@
 #include "cnf.h"
 
+/*
 cnf* readDIMACS(char* path) {
     cnf* c = malloc(sizeof(cnf));
     CHECK(c, "Failed to alloc CNF struct\n");
@@ -325,3 +326,4 @@ void sortlastnum(cnf* c, u64 N) {
     free(swaparr);
     free(d);
 }
+ */

cnf.h

@@ -4,6 +4,7 @@
 #include <string.h>
 #include <stdio.h>
 
+
 #define CHECK(X, Y) if (X == NULL) { \
     printf(Y);                       \
     return NULL;                     \
@@ -20,6 +21,7 @@ typedef struct {
     u8* pars;
 } cnf;
 
+
 cnf* readDIMACS(char* path);
 
 void printcnf(cnf* c);

gpusolver.c

@@ -176,6 +176,11 @@ i32 gpusolve(cnf* c) {
         exit(1);
     }
 
+    size_t maxworkgrpu = 0;
+    res = clGetKernelWorkGroupInfo(kernel, deviceid, CL_KERNEL_WORK_GROUP_SIZE, sizeof(size_t), &maxworkgrpu, NULL);
+
+    printf("Max work group size: %lu\n", maxworkgrpu);
+
     res = clSetKernelArg(kernel, 0, sizeof(cl_mem), (void*) &gpuheader);
     res = clSetKernelArg(kernel, 1, sizeof(cl_mem), (void*) &gpulvars);
     res = clSetKernelArg(kernel, 2, sizeof(cl_mem), (void*) &gpuvariables);

main.c

@@ -1,10 +1,256 @@
 #include <stdio.h>
-#include "cnf.h"
+// #include "gpusolver.h"
-#include "gpusolver.h"
 #include "time.h"
-#include "tests/masterTest.h"
+// #include "tests/masterTest.h"
+#include "gmp.h"
+#include "rng.h"
+#include "ncnf.h"
 
+#define ADD (0)
+#define CMP (1)
+#define CHK (2)
+
+void ctrthings2(cnf* c, u32* state, u32* ctr, u32* max) {
+    u32 wcnt = 1U + (c->cnts[0] / 32U);
+    u32* mode = state;
+    u32* index = state + 1;
+    u32* addval = state + 2;
+
+    u32 varcnt = c->cnts[0] - 1;
+    u32 chkmsk = 0xFFFFFFFFU * (*mode == CHK);
+    u32 chkcls = *index & chkmsk;
+    u32 chkind = c->clausedat[3 * chkcls] + (*addval & chkmsk);
+    u32 var = c->variables[chkind];
+    u8 par = c->parities[chkind];
+    u32 vword = (varcnt - var) >> 5U;
+    u32 vbit = (varcnt - var) & 0b11111U;
+    u8 corpar = (ctr[vword] >> vbit) & 1U;
+    u8 isvalid = (par == corpar);
+    u8 islvar = ((*addval + 1) == c->clausedat[3 * chkcls + 1]);
+    if (*mode == CHK) {
+        // printf("> %u %u\n", *index, *addval);
+        // printf("%u %u %u\n", var, par, corpar);
+        // printf("%u %u\n", chkcls, chkind);
+        // printf("%u %u %u\n", islvar, isvalid, c->clausedat[3 * chkcls + 1]);
+
+        /*
+         * if last var
+         *  if valid, add
+         *  if invalid, iterate clause
+         * else
+         *  if valid, iterate addval up to len
+         *  if invalid, iterate claus
+         * If current var is valid:
+         *  if last var:
+         *
+         */
+        if (isvalid) {
+            if (islvar) {
+                u32 j = c->clausedat[3 * chkcls + 2];
+                *mode -= 2;
+                printf("j: %u\n", j);
+                *index = j >> 5U;
+                *addval = 1U << (j & 0b11111U);
+            } else {
+                *addval += 1U;
+            }
+        } else {
+            *addval = 0;
+            *index += 1U;
+            if (*index == c->cnts[1]) {
+                printf("SAT\n");
+                *mode = 4;
+                return;
+            }
+        }
+    } /* else {
+
+        //printf("YEET\n");
+        if (*index >= wcnt) printf("FUCK\n");
+        u32 nval = ctr[*index] + *addval; // Find the result of the current step if it was addition
+        *addval = (nval < ctr[*index]) * (*mode == ADD); // If in add mode, set addval to carry, else set 0
+        ctr[*index] = nval * (*mode == ADD) + ctr[*index] * (*mode != ADD); // If in add mode, set new ctr val, otherwise leave unchanged
+        *addval -= (ctr[*index] < max[*index]) * (*mode == CMP); // If in comparison mode, decrement addval if less than
+        *addval += (ctr[*index] > max[*index]) * (*mode == CMP); // If in comparison mode, increment addval if greater than
+        bool addcond = (*addval == 0) | (*index == (wcnt - 1)); // Exit condition for the ADD state: If addval is zero (no carry) or we're at the last word
+        bool cmpcond = (*addval != 0) | (*index == 0); // Exit condition for the CMP state: if addval is nonzero (lt or gt) or we're at the least significant word
+        if (*mode == CMP && cmpcond && *addval != -1) { // If in cmpmode and the comparison result is not less than, unsat
+            printf("UNSAT\n");
+            *mode = 4;
+            return;
+        }
+        bool cmpdone = cmpcond * (*mode == CMP); // if comparison completion conditions are satisfied and in CMP mode
+        u32 addindex = (*index + 1) * !addcond + (wcnt - 1) * addcond; // if add completion is satisfied, set index to most significant word, else increment by 1
+        *index = addindex * (*mode == ADD) + (*index - 1) * (*mode == CMP); // If in add mode, use addindex; if in cmp mode, decrement index by 1
+        *index *= !cmpdone;
+        // Leave adval alone if:
+        //  not in add mode
+        //      add mode isn't done
+        //  not in cmp mode
+        //      cmp mode isn't done
+        //
+        *addval *= !(((addcond) & (*mode != ADD)) & cmpdone); // If add is complete, zero addval, else leave unchanged
+        *mode += addcond * (*mode == ADD) + cmpdone; // If in add mode and add completion is reached, increment mode. If in cmp mode and cmp completion reached, increment mode.
+
+    }
+    */
+
+    u32 cmpaddind = *index * (*mode != CHK);
+    if (cmpaddind >= wcnt) printf("FUCK\n");
+    u32 nval = ctr[cmpaddind] + *addval; // Find the result of the current step if it was addition
+    *addval = (nval < ctr[cmpaddind]) * (*mode == ADD) + (*addval) * (*mode == CHK); // If in add mode, set addval to carry. If in cmp mode, set to 0. If in check mode, leave alone.
+    ctr[cmpaddind] = nval * (*mode == ADD) + ctr[cmpaddind] * (*mode != ADD); // If in add mode, set new ctr val, otherwise leave unchanged
+    *addval -= (ctr[cmpaddind] < max[cmpaddind]) * (*mode == CMP); // If in comparison mode, decrement addval if less than
+    *addval += (ctr[cmpaddind] > max[cmpaddind]) * (*mode == CMP); // If in comparison mode, increment addval if greater than
+    bool addcond = (*addval == 0) | (cmpaddind == (wcnt - 1)); // Exit condition for the ADD state: If addval is zero (no carry) or we're at the last word
+    bool cmpcond = (*addval != 0) | (cmpaddind == 0); // Exit condition for the CMP state: if addval is nonzero (lt or gt) or we're at the least significant word
+    if (*mode == CMP && cmpcond && *addval != -1) { // If in cmpmode and the comparison result is not less than, unsat
+        printf("UNSAT\n");
+        *mode = 4;
+        return;
+    }
+    bool cmpdone = cmpcond & (*mode == CMP); // if comparison completion conditions are satisfied and in CMP mode
+    u32 addindex = (cmpaddind + 1) * !addcond + (wcnt - 1) * addcond; // if add completion is satisfied, set index to most significant word, else increment by 1
+    *index = addindex * (*mode == ADD) + (*index - (*mode == CMP)) * (*mode != ADD); // If in add mode, use addindex; if in cmp mode, decrement index by 1
+    *index *= !cmpdone;
+    // Leave adval alone if:
+    //  not in add mode
+    //      add mode isn't done
+    //  not in cmp mode
+    //      cmp mode isn't done
+    //
+    *addval *= !(((addcond) & (*mode == ADD)) | cmpdone); // If add is complete, or cmp is complete, zero. Else leave unchanged.
+    *mode += addcond * (*mode == ADD) + cmpdone; // If in add mode and add completion is reached, increment mode. If in cmp mode and cmp completion reached, increment mode.
+
+}
+
+void printbits(unsigned a) {
+    for (unsigned i = 0; i < 32; ++i) {
+        unsigned ind = 31 - i;
+        printf("%u", (a >> ind) & 1U);
+    }
+}
+
+#define TESTS (274877906944LU >> 10U)
+#define CSZE (83LU)
+#define eqprob (0.01f)
 int main() {
+    /*
+    printf("Tests: %lu\n", TESTS);
+    rngstate rng;
+    u64 rseed = utime();
+    seed(&rng, rseed);
+    printf("Seed: %lu\n", rseed);
+
+    mpz_t a, b, c, d, e;
+
+    mpz_inits(a, b, c, d, e, NULL);
+
+    u32 ctrp[CSZE];
+    u32 maxp[CSZE];
+
+    u32* ctr = ctrp;
+    u32* max = maxp;
+
+    u32 state[3];
+    u32 hdr[2];
+
+    char buf[4096];
+
+
+    for (u64 i = 0; i < TESTS; ++i) {
+        memset(ctr, 0, sizeof(u32) * CSZE);
+        memset(max, 0, sizeof(u32) * CSZE);
+
+        //u32 lenval = ru32(&rng) % CSZE;
+        u32 lenval = CSZE;
+        if (lenval < 2) lenval = 2;
+        u32 jval = ru32(&rng) % ((lenval - 1) * 32);
+
+        for (u32 j = 0; j < lenval - 2; ++j) {
+            ctr[j] = ru32(&rng);
+            max[j] = ru32(&rng);
+        }
+
+        mpz_import(a, lenval, -1, sizeof(u32), 0, 0, ctr);
+        mpz_import(b, lenval, -1, sizeof(u32), 0, 0, max);
+
+        state[1] = jval >> 5U;
+        state[2] = jval & 0b11111U;
+        state[2] = 1U << state[2];
+        mpz_ui_pow_ui(c, 2, jval);
+
+
+        if (rf32(&rng) < eqprob && mpz_cmp) {
+            mpz_sub(a, b, c);
+            if (mpz_sgn(a) != -1) {
+                mpz_export(ctr, NULL, -1, sizeof(u32), 0, 0, a);
+            } else {
+                mpz_import(a, lenval, -1, sizeof(u32), 0, 0, ctr);
+            }
+        }
+
+        mpz_add(d, a, c);
+
+        state[0] = 0;
+
+        while (state[0] < 2U) {
+            ctrthings2(hdr, lenval, state, ctr, max);
+        }
+
+        mpz_import(c, lenval, -1, sizeof(u32), 0, 0, ctr);
+
+
+
+        i32 res = mpz_cmp(d, b);
+
+        if (res == -1) {
+            if (state[2] != -1) {
+                printf("Fuck2 %lu\n", i);
+                printf("d: ");
+                mpz_out_str(stdout, 10, d);
+
+                printf("\nb: ");
+                mpz_out_str(stdout, 10, b);
+                printf("\n");
+                printf("mode: %u\n", state[2]);
+                exit(0);
+            }
+        } else {
+            if (state[2] != 2) {
+                printf("Fuck3 %lu\n", i);
+                printf("d: ");
+                mpz_out_str(stdout, 10, d);
+                printf("\nc: ");
+                mpz_out_str(stdout, 10, c);
+                printf("\nb: ");
+                mpz_out_str(stdout, 10, b);
+                printf("\n");
+                printf("mode: %u\n", state[2]);
+                exit(0);
+            }
+        }
+
+        res = mpz_cmp(d, c);
+        if (res != 0) {
+            printf("Fuck %lu\na: ", i);
+            mpz_out_str(stdout, 10, a);
+            printf("\nd: ");
+            mpz_out_str(stdout, 10, d);
+            printf("\nc: ");
+            mpz_out_str(stdout, 10, c);
+            printf("\n%u %u\n", lenval, jval);
+            exit(0);
+        }
+    }
+
+
+
+    mpz_clears(a, b, c, d, e, NULL);
+     */
+
+
+
     /*
     srand( utime());
     u32 cnt = 4096;
@@ -28,7 +274,43 @@ int main() {
      */
 
     /* Expects a path to a DIMACS file */
-    cnf* c = readDIMACS("/home/lev/Downloads/uf50/uf50-088.cnf");
+
+
+    cnf* c = readDIMACS("/home/lev/Downloads/uf20/uf20-022.cnf");
+    sortlastnum(c);
+
+    // printcnf(c);
+    u32 wcnt = 1U + (c->cnts[0] / 32U);
+
+
+    u32* ctr = calloc(wcnt, sizeof(u32));
+    u32* max = calloc(wcnt, sizeof(u32));
+
+    max[c->cnts[0] >> 5U] = 1U << (c->cnts[0] & 0b11111U);
+
+    u32 state[3];
+
+    state[0] = 2;
+    state[1] = state[2] = 0;
+
+
+    u32 mtr = 0;
+    while (state[0] < 3U) {
+        u32 cmd = state[0];
+        ctrthings2(c, state, ctr, max);
+        if (state[0] != 2 && cmd == 2) {
+            //printf("\n");
+            //for (unsigned i = wcnt - 1; i < wcnt; --i) printbits(ctr[i]);
+            //printf("\n");
+            //mtr++;
+            //if (mtr == 10) exit(15);
+            // printf("%u %u %u\n", state[0], state[1], state[2]);
+        }
+    }
+
+    return 0;
+    /*
+
     // printf("%u\n", c->litcnt);
 
     // for (u32 i = 0; i < c->clausecnt; ++i) printf("%u ", c->lastvars[i]);
@@ -46,5 +328,6 @@ int main() {
     // if (res == 1) break;
     freecnf(c);
 
-    return 0;
+    return 0;*/
+
 }

ncnf.c

@@ -0,0 +1,210 @@
+#include "ncnf.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
+    u32* varcnt = &(c->cnts[0]);
+    *varcnt = 0;
+    while (((u8) (*temp - '0')) < 10) {
+        *varcnt *= 10;
+        *varcnt += (*temp - '0');
+        temp++;
+    }
+
+    // Skip any trailing whitespace
+    while (*temp == ' ') temp++;
+
+    // Read in clausecnt
+    u32* clausecnt = &(c->cnts[1]);
+    *clausecnt = 0;
+    while (((u8) (*temp - '0')) < 10) {
+        *clausecnt *= 10;
+        *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) * ((*varcnt) * 12LU) + 10;
+    char* nbuf = realloc(buf, nsize);
+    CHECK(nbuf, "Failed to realloc read buffer\n")
+
+    buf = nbuf;
+
+    u32 cnt = 0;
+    u32 cap = 32;
+
+    c->clausedat = calloc(*clausecnt, sizeof(u32) * 3);
+    CHECK(c->clausedat, "Failed to allocate clause data\n")
+    c->variables = calloc(cap, sizeof(u32));
+    CHECK(c->variables, "Failed to allocate literal variables\n")
+    c->parities = calloc(cap, sizeof(u8));
+    CHECK(c->parities, "Failed to allocate literal parities\n");
+
+    for (u32 i = 0; i < *clausecnt; ++i) {
+        c->clausedat[3 * i] = cnt;
+        CHECK(fgets(buf, nsize, f), "Failed to read clause\n");
+
+        temp = buf;
+        while (*temp == ' ') temp++;
+
+        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->parities = realloc(c->parities, sizeof(u8) * ncap);
+                memset(c->parities + cap, 0, sizeof(u8) * cap);
+                CHECK(c->parities, "Failed to realloc parity array\n");
+                cap = ncap;
+            }
+
+            if (*temp == '-') {
+                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->parities[cnt] = tr;
+                c->variables[cnt] -= 1;
+                if (c->clausedat[3 * i + 2] < c->variables[cnt]) c->clausedat[3 * i + 2] = 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;
+        }
+        c->clausedat[3 * i + 1] = cnt - c->clausedat[3 * i];
+    }
+
+    // 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 < *clausecnt; ++i) {
+        c->clausedat[3 * i + 2] += 1U;
+        c->clausedat[3 * i + 2] = *varcnt - c->clausedat[3 * i + 2];
+    }
+
+    // Realloc the arrays to exactly match the number of literals
+    c->variables = realloc(c->variables, sizeof(u32) * cnt);
+    c->parities = realloc(c->parities, sizeof(u8) * cnt);
+    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->cnts[0], c->cnts[1]);
+    for (u32 i = 0; i < c->cnts[1]; ++i) {
+        for (u32 j = 0; j < c->clausedat[3 * i + 1]; ++j) {
+            if (c->parities[c->clausedat[3 * i] + j] == 0) printf("-");
+            printf("%u ", c->variables[c->clausedat[3 * i] + j] + 1);
+        }
+        printf("0\n");
+    }
+}
+
+void freecnf(cnf* c) {
+    free(c->clausedat);
+    free(c->variables);
+    free(c->parities);
+    free(c);
+}
+
+const u64 RBITS = 8;
+const u64 RSIZE = 1LU << RBITS;
+const u64 RLVLS = (31LU / RBITS) + 1;
+const u64 RMASK = RSIZE - 1;
+
+void sortlastnum(cnf* c) {
+    bool v = false;
+    u32* d = calloc(c->cnts[1], sizeof(u32) * 3);
+
+    u32 qb[RLVLS * RSIZE];
+    u32 qe[RLVLS * RSIZE];
+    memset(qb, 0, sizeof(u32) * RLVLS * RSIZE);
+    memset(qe, 0, sizeof(u32) * RLVLS * RSIZE);
+
+    for (u32 pass = 0; pass < RLVLS; ++pass) {
+        u32 shift = pass * RBITS;
+
+        for (u32 i = 0; i < c->cnts[1]; ++i) {
+            u32 val = UINT32_MAX - c->clausedat[3 * i + 2];
+            u32 ind = (val >> shift) & RMASK;
+            qe[pass * RSIZE + ind]++;
+        }
+
+        u32 uc = 0;
+        for (u32 i = 0; i < RSIZE; ++i) if (qe[pass * RSIZE + i]) uc++;
+        if (uc == 1) continue;
+
+        qb[pass * RSIZE] = 0;
+        for (u32 i = 1; i < RSIZE; ++i) {
+            qb[pass * RSIZE + i] = qb[pass * RSIZE + i - 1] + qe[pass * RSIZE + i - 1];
+        }
+
+        for (u32 i = 0; i < c->cnts[1]; ++i) {
+            u32 val = UINT32_MAX - c->clausedat[3 * i + 2];
+            u32 ind = (val >> shift) & RMASK;
+            memcpy(d + 3 * qb[pass * RSIZE + ind], c->clausedat + 3 * i, sizeof(u32) * 3);
+            qb[pass * RSIZE + ind]++;
+            __builtin_prefetch(d + 3 * qb[pass * RSIZE + ind] + 1);
+        }
+
+        u32* tptr = c->clausedat;
+        c->clausedat = d;
+        d = tptr;
+        v = !v;
+    }
+
+    if (v) {
+        u32* tptr = c->clausedat;
+        c->clausedat = d;
+        d = tptr;
+        memcpy(c->clausedat, d, sizeof(u32) * 3 * c->cnts[1]);
+    }
+
+    free(d);
+}

ncnf.h

@@ -0,0 +1,25 @@
+#pragma once
+#include "types.h"
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+
+#define CHECK(X, Y) if (X == NULL) { \
+    printf(Y);                       \
+    return NULL;                     \
+}
+
+typedef struct {
+    u32 cnts[2];    // { varcnt, clausecnt }
+    u32* clausedat; // { ind, len, jval }
+    u32* variables;
+    u8* parities;
+} cnf;
+
+cnf* readDIMACS(char* path);
+
+void printcnf(cnf* c);
+
+void sortlastnum(cnf* c);
+
+void freecnf(cnf* c);

rng.c

@@ -0,0 +1,121 @@
+#include "rng.h"
+#include <string.h>
+#include "time.h"
+#include "stdlib.h"
+
+static inline u64 rotl(const u64 x, i32 k) {
+    return (x << k) | (x >> (64 - k));
+}
+
+void seed(rngstate* s, u64 seed) {
+    /* Unrolled & modified version of splitmix64
+     * https://prng.di.unimi.it/splitmix64.c */
+    u64 z = (seed += 0x9e3779b97f4a7c15);
+    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+    z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+    s->state[0] = z ^ (z >> 31);
+
+    z = (seed += 0x9e3779b97f4a7c15);
+    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+    z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+    s->state[1] = z ^ (z >> 31);
+
+    z = (seed += 0x9e3779b97f4a7c15);
+    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+    z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+    s->state[2] = z ^ (z >> 31);
+
+    z = (seed += 0x9e3779b97f4a7c15);
+    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+    z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+    s->state[3] = z ^ (z >> 31);
+
+    z = (seed += 0x9e3779b97f4a7c15);
+    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+    z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+    s->state[4] = z ^ (z >> 31);
+
+    z = (seed += 0x9e3779b97f4a7c15);
+    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+    z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+    s->state[5] = z ^ (z >> 31);
+
+    z = (seed += 0x9e3779b97f4a7c15);
+    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+    z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+    s->state[6] = z ^ (z >> 31);
+
+    z = (seed + 0x9e3779b97f4a7c15);
+    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+    z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+    s->state[7] = z ^ (z >> 31);
+
+    s->ready = true;
+}
+
+u64 ru64(rngstate* s) {
+    if (!s->ready) seed(s, utime());
+    const u64 res = rotl(s->state[1] * 5, 7) * 9;
+    const u64 t = s->state[1] << 11U;
+
+    s->state[2] ^= s->state[0];
+    s->state[5] ^= s->state[1];
+    s->state[1] ^= s->state[2];
+    s->state[7] ^= s->state[3];
+    s->state[3] ^= s->state[4];
+    s->state[4] ^= s->state[5];
+    s->state[0] ^= s->state[6];
+    s->state[6] ^= s->state[7];
+
+    s->state[6] ^= t;
+
+    s->state[7] = rotl(s->state[7], 21);
+
+    return res;
+}
+
+u64 ru64pow(rngstate* s, u64 max) {
+    u64 v = ru64(s) % max;
+    v++;
+    return max / v;
+}
+
+u32 ru32(rngstate* s) {
+    return (u32) ru64(s);
+}
+
+f64 rf64(rngstate* s) {
+    return ((f64) ru64(s)) / ((f64) UINT64_MAX);
+}
+
+f32 rf32(rngstate* s) {
+    return ((f32) ru32(s)) / ((f32) UINT32_MAX);
+}
+
+void jump(rngstate* s) {
+    static const u64 JUMP[] = { 0x33ed89b6e7a353f9, 0x760083d7955323be, 0x2837f2fbb5f22fae, 0x4b8c5674d309511c, 0xb11ac47a7ba28c25, 0xf1be7667092bcc1c, 0x53851efdb6df0aaf, 0x1ebbc8b23eaf25db };
+    u64 t[8];
+    memset(t, 0, sizeof(u64) * 8);
+    for (i32 i = 0; i < 8; ++i) {
+        for (i32 b = 0; b < 64; ++b) {
+            if (JUMP[i] & 1LU << b) for (i32 w = 0; w < 8; ++w) t[w] ^= s->state[w];
+            ru64(s);
+        }
+    }
+
+    memcpy(s->state, t, sizeof(u64) * 8);
+}
+
+void longjump(rngstate* s) {
+    static const u64 LONGJUMP[] = { 0x11467fef8f921d28, 0xa2a819f2e79c8ea8, 0xa8299fc284b3959a, 0xb4d347340ca63ee1, 0x1cb0940bedbff6ce, 0xd956c5c4fa1f8e17, 0x915e38fd4eda93bc, 0x5b3ccdfa5d7daca5 };
+    u64 t[8];
+    memset(t, 0, sizeof(u64) * 8);
+    for (i32 i = 0; i < 8; ++i) {
+        for (i32 b = 0; b < 64; ++b) {
+            if (LONGJUMP[i] & 1LU << b) for (i32 w = 0; w < 8; ++w) t[w] ^= s->state[w];
+            ru64(s);
+        }
+    }
+
+    memcpy(s->state, t, sizeof(u64) * 8);
+}

rng.h

@@ -0,0 +1,15 @@
+#pragma once
+#include "types.h"
+
+typedef struct {
+    bool ready;
+    u64 state[8];
+} rngstate;
+
+void seed(rngstate* s, u64 seed);
+u64 ru64(rngstate* s);
+u32 ru32(rngstate* s);
+f64 rf64(rngstate* s);
+f32 rf32(rngstate* s);
+void jump(rngstate* s);
+u64 ru64pow(rngstate* s, u64 max);

tests/masterTest.c

@@ -33,7 +33,8 @@ i32 runuf20() {
 
         cnf* c = readDIMACS(buf);
 
-        sortlastnum(c, c->litcnt);
+        // TODO: Uncomment
+        // sortlastnum(c, c->litcnt);
 
         u64 start = utime();
         i32 res = gpusolve(c);
@@ -59,7 +60,8 @@ i32 runuf50() {
 
         cnf* c = readDIMACS(buf);
 
-        sortlastnum(c, c->litcnt);
+        // TODO: Uncomment
+        // sortlastnum(c, c->litcnt);
 
         u64 start = utime();
         i32 res = gpusolve(c);
@@ -85,7 +87,8 @@ i32 runuuf50() {
 
         cnf* c = readDIMACS(buf);
 
-        sortlastnum(c, c->litcnt);
+        // TODO: Uncomment
+        // sortlastnum(c, c->litcnt);
 
         u64 start = utime();
         i32 res = gpusolve(c);