summaryrefslogtreecommitdiff
path: root/vere/ext/h2o/patches/h2o-2.2.6/deps/klib
diff options
context:
space:
mode:
Diffstat (limited to 'vere/ext/h2o/patches/h2o-2.2.6/deps/klib')
-rw-r--r--vere/ext/h2o/patches/h2o-2.2.6/deps/klib/kopen.c344
-rw-r--r--vere/ext/h2o/patches/h2o-2.2.6/deps/klib/ksw.c637
-rw-r--r--vere/ext/h2o/patches/h2o-2.2.6/deps/klib/test/kbit_test.c141
3 files changed, 1122 insertions, 0 deletions
diff --git a/vere/ext/h2o/patches/h2o-2.2.6/deps/klib/kopen.c b/vere/ext/h2o/patches/h2o-2.2.6/deps/klib/kopen.c
new file mode 100644
index 0000000..095fb4e
--- /dev/null
+++ b/vere/ext/h2o/patches/h2o-2.2.6/deps/klib/kopen.c
@@ -0,0 +1,344 @@
+#include <stdio.h>
+#include <fcntl.h>
+#include <errno.h>
+#include <ctype.h>
+#include <unistd.h>
+#include <string.h>
+#include <stdlib.h>
+#include <signal.h>
+#include <sys/types.h>
+#include <sys/wait.h>
+#ifndef _WIN32
+#include <netdb.h>
+#include <arpa/inet.h>
+#include <sys/socket.h>
+#endif
+
+#ifdef _WIN32
+#define _KO_NO_NET
+#endif
+
+#ifndef _KO_NO_NET
+static int socket_wait(int fd, int is_read)
+{
+ fd_set fds, *fdr = 0, *fdw = 0;
+ struct timeval tv;
+ int ret;
+ tv.tv_sec = 5; tv.tv_usec = 0; // 5 seconds time out
+ FD_ZERO(&fds);
+ FD_SET(fd, &fds);
+ if (is_read) fdr = &fds;
+ else fdw = &fds;
+ ret = select(fd+1, fdr, fdw, 0, &tv);
+ if (ret == -1) perror("select");
+ return ret;
+}
+
+static int socket_connect(const char *host, const char *port)
+{
+#define __err_connect(func) do { perror(func); freeaddrinfo(res); return -1; } while (0)
+
+ int on = 1, fd;
+ struct linger lng = { 0, 0 };
+ struct addrinfo hints, *res = 0;
+ memset(&hints, 0, sizeof(struct addrinfo));
+ hints.ai_family = AF_UNSPEC;
+ hints.ai_socktype = SOCK_STREAM;
+ if (getaddrinfo(host, port, &hints, &res) != 0) __err_connect("getaddrinfo");
+ if ((fd = socket(res->ai_family, res->ai_socktype, res->ai_protocol)) == -1) __err_connect("socket");
+ if (setsockopt(fd, SOL_SOCKET, SO_REUSEADDR, &on, sizeof(on)) == -1) __err_connect("setsockopt");
+ if (setsockopt(fd, SOL_SOCKET, SO_LINGER, &lng, sizeof(lng)) == -1) __err_connect("setsockopt");
+ if (connect(fd, res->ai_addr, res->ai_addrlen) != 0) __err_connect("connect");
+ freeaddrinfo(res);
+ return fd;
+#undef __err_connect
+}
+
+static int http_open(const char *fn)
+{
+ char *p, *proxy, *q, *http_host, *host, *port, *path, *buf;
+ int fd, ret, l;
+
+ /* parse URL; adapted from khttp_parse_url() in knetfile.c */
+ if (strstr(fn, "http://") != fn) return 0;
+ // set ->http_host
+ for (p = (char*)fn + 7; *p && *p != '/'; ++p);
+ l = p - fn - 7;
+ http_host = calloc(l + 1, 1);
+ strncpy(http_host, fn + 7, l);
+ http_host[l] = 0;
+ for (q = http_host; *q && *q != ':'; ++q);
+ if (*q == ':') *q++ = 0;
+ // get http_proxy
+ proxy = getenv("http_proxy");
+ // set host, port and path
+ if (proxy == 0) {
+ host = strdup(http_host); // when there is no proxy, server name is identical to http_host name.
+ port = strdup(*q? q : "80");
+ path = strdup(*p? p : "/");
+ } else {
+ host = (strstr(proxy, "http://") == proxy)? strdup(proxy + 7) : strdup(proxy);
+ for (q = host; *q && *q != ':'; ++q);
+ if (*q == ':') *q++ = 0;
+ port = strdup(*q? q : "80");
+ path = strdup(fn);
+ }
+
+ /* connect; adapted from khttp_connect() in knetfile.c */
+ l = 0;
+ fd = socket_connect(host, port);
+ buf = calloc(0x10000, 1); // FIXME: I am lazy... But in principle, 64KB should be large enough.
+ l += sprintf(buf + l, "GET %s HTTP/1.0\r\nHost: %s\r\n", path, http_host);
+ l += sprintf(buf + l, "\r\n");
+ write(fd, buf, l);
+ l = 0;
+ while (read(fd, buf + l, 1)) { // read HTTP header; FIXME: bad efficiency
+ if (buf[l] == '\n' && l >= 3)
+ if (strncmp(buf + l - 3, "\r\n\r\n", 4) == 0) break;
+ ++l;
+ }
+ buf[l] = 0;
+ if (l < 14) { // prematured header
+ close(fd);
+ fd = -1;
+ }
+ ret = strtol(buf + 8, &p, 0); // HTTP return code
+ if (ret != 200) {
+ close(fd);
+ fd = -1;
+ }
+ free(buf); free(http_host); free(host); free(port); free(path);
+ return fd;
+}
+
+typedef struct {
+ int max_response, ctrl_fd;
+ char *response;
+} ftpaux_t;
+
+static int kftp_get_response(ftpaux_t *aux)
+{
+ unsigned char c;
+ int n = 0;
+ char *p;
+ if (socket_wait(aux->ctrl_fd, 1) <= 0) return 0;
+ while (read(aux->ctrl_fd, &c, 1)) { // FIXME: this is *VERY BAD* for unbuffered I/O
+ if (n >= aux->max_response) {
+ aux->max_response = aux->max_response? aux->max_response<<1 : 256;
+ aux->response = realloc(aux->response, aux->max_response);
+ }
+ aux->response[n++] = c;
+ if (c == '\n') {
+ if (n >= 4 && isdigit(aux->response[0]) && isdigit(aux->response[1]) && isdigit(aux->response[2])
+ && aux->response[3] != '-') break;
+ n = 0;
+ continue;
+ }
+ }
+ if (n < 2) return -1;
+ aux->response[n-2] = 0;
+ return strtol(aux->response, &p, 0);
+}
+
+static int kftp_send_cmd(ftpaux_t *aux, const char *cmd, int is_get)
+{
+ if (socket_wait(aux->ctrl_fd, 0) <= 0) return -1; // socket is not ready for writing
+ write(aux->ctrl_fd, cmd, strlen(cmd));
+ return is_get? kftp_get_response(aux) : 0;
+}
+
+static int ftp_open(const char *fn)
+{
+ char *p, *host = 0, *port = 0, *retr = 0;
+ char host2[80], port2[10];
+ int v[6], l, fd = -1, ret, pasv_port, pasv_ip[4];
+ ftpaux_t aux;
+
+ /* parse URL */
+ if (strstr(fn, "ftp://") != fn) return 0;
+ for (p = (char*)fn + 6; *p && *p != '/'; ++p);
+ if (*p != '/') return 0;
+ l = p - fn - 6;
+ port = strdup("21");
+ host = calloc(l + 1, 1);
+ strncpy(host, fn + 6, l);
+ retr = calloc(strlen(p) + 8, 1);
+ sprintf(retr, "RETR %s\r\n", p);
+
+ /* connect to ctrl */
+ memset(&aux, 0, sizeof(ftpaux_t));
+ aux.ctrl_fd = socket_connect(host, port);
+ if (aux.ctrl_fd == -1) goto ftp_open_end; /* fail to connect ctrl */
+
+ /* connect to the data stream */
+ kftp_get_response(&aux);
+ kftp_send_cmd(&aux, "USER anonymous\r\n", 1);
+ kftp_send_cmd(&aux, "PASS kopen@\r\n", 1);
+ kftp_send_cmd(&aux, "TYPE I\r\n", 1);
+ kftp_send_cmd(&aux, "PASV\r\n", 1);
+ for (p = aux.response; *p && *p != '('; ++p);
+ if (*p != '(') goto ftp_open_end;
+ ++p;
+ sscanf(p, "%d,%d,%d,%d,%d,%d", &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
+ memcpy(pasv_ip, v, 4 * sizeof(int));
+ pasv_port = (v[4]<<8&0xff00) + v[5];
+ kftp_send_cmd(&aux, retr, 0);
+ sprintf(host2, "%d.%d.%d.%d", pasv_ip[0], pasv_ip[1], pasv_ip[2], pasv_ip[3]);
+ sprintf(port2, "%d", pasv_port);
+ fd = socket_connect(host2, port2);
+ if (fd == -1) goto ftp_open_end;
+ ret = kftp_get_response(&aux);
+ if (ret != 150) {
+ close(fd);
+ fd = -1;
+ }
+ close(aux.ctrl_fd);
+
+ftp_open_end:
+ free(host); free(port); free(retr); free(aux.response);
+ return fd;
+}
+#endif /* !defined(_KO_NO_NET) */
+
+static char **cmd2argv(const char *cmd)
+{
+ int i, beg, end, argc;
+ char **argv, *p, *q, *str;
+ end = strlen(cmd);
+ for (i = end - 1; i >= 0; --i)
+ if (!isspace(cmd[i])) break;
+ end = i + 1;
+ for (beg = 0; beg < end; ++beg)
+ if (!isspace(cmd[beg])) break;
+ if (beg == end) return 0;
+ for (i = beg + 1, argc = 0; i < end; ++i)
+ if (isspace(cmd[i]) && !isspace(cmd[i-1]))
+ ++argc;
+ argv = (char**)calloc(argc + 2, sizeof(void*));
+ argv[0] = str = (char*)calloc(end - beg + 1, 1);
+ strncpy(argv[0], cmd + beg, end - beg);
+ for (i = argc = 1, q = p = str; i < end - beg; ++i)
+ if (isspace(str[i])) str[i] = 0;
+ else if (str[i] && str[i-1] == 0) argv[argc++] = &str[i];
+ return argv;
+}
+
+#define KO_STDIN 1
+#define KO_FILE 2
+#define KO_PIPE 3
+#define KO_HTTP 4
+#define KO_FTP 5
+
+typedef struct {
+ int type, fd;
+ pid_t pid;
+} koaux_t;
+
+void *kopen(const char *fn, int *_fd)
+{
+ koaux_t *aux = 0;
+ *_fd = -1;
+ if (strstr(fn, "http://") == fn) {
+ aux = calloc(1, sizeof(koaux_t));
+ aux->type = KO_HTTP;
+ aux->fd = http_open(fn);
+ } else if (strstr(fn, "ftp://") == fn) {
+ aux = calloc(1, sizeof(koaux_t));
+ aux->type = KO_FTP;
+ aux->fd = ftp_open(fn);
+ } else if (strcmp(fn, "-") == 0) {
+ aux = calloc(1, sizeof(koaux_t));
+ aux->type = KO_STDIN;
+ aux->fd = STDIN_FILENO;
+ } else {
+ const char *p, *q;
+ for (p = fn; *p; ++p)
+ if (!isspace(*p)) break;
+ if (*p == '<') { // pipe open
+ int need_shell, pfd[2];
+ pid_t pid;
+ // a simple check to see if we need to invoke a shell; not always working
+ for (q = p + 1; *q; ++q)
+ if (ispunct(*q) && *q != '.' && *q != '_' && *q != '-' && *q != ':')
+ break;
+ need_shell = (*q != 0);
+ pipe(pfd);
+ pid = vfork();
+ if (pid == -1) { /* vfork() error */
+ close(pfd[0]); close(pfd[1]);
+ return 0;
+ }
+ if (pid == 0) { /* the child process */
+ char **argv; /* FIXME: I do not know if this will lead to a memory leak */
+ close(pfd[0]);
+ dup2(pfd[1], STDOUT_FILENO);
+ close(pfd[1]);
+ if (!need_shell) {
+ argv = cmd2argv(p + 1);
+ execvp(argv[0], argv);
+ free(argv[0]); free(argv);
+ } else execl("/bin/sh", "sh", "-c", p + 1, NULL);
+ exit(1);
+ } else { /* parent process */
+ close(pfd[1]);
+ aux = calloc(1, sizeof(koaux_t));
+ aux->type = KO_PIPE;
+ aux->fd = pfd[0];
+ aux->pid = pid;
+ }
+ } else {
+#ifdef _WIN32
+ *_fd = open(fn, O_RDONLY | O_BINARY);
+#else
+ *_fd = open(fn, O_RDONLY);
+#endif
+ if (*_fd) {
+ aux = calloc(1, sizeof(koaux_t));
+ aux->type = KO_FILE;
+ aux->fd = *_fd;
+ }
+ }
+ }
+ *_fd = aux->fd;
+ return aux;
+}
+
+int kclose(void *a)
+{
+ koaux_t *aux = (koaux_t*)a;
+ if (aux->type == KO_PIPE) {
+ int status;
+ pid_t pid;
+ pid = waitpid(aux->pid, &status, WNOHANG);
+ if (pid != aux->pid) kill(aux->pid, 15);
+ }
+ return 0;
+}
+
+#ifdef _KO_MAIN
+#define BUF_SIZE 0x10000
+int main(int argc, char *argv[])
+{
+ void *x;
+ int l, fd;
+ unsigned char buf[BUF_SIZE];
+ FILE *fp;
+ if (argc == 1) {
+ fprintf(stderr, "Usage: kopen <file>\n");
+ return 1;
+ }
+ x = kopen(argv[1], &fd);
+ fp = fdopen(fd, "r");
+ if (fp == 0) {
+ fprintf(stderr, "ERROR: fail to open the input\n");
+ return 1;
+ }
+ do {
+ if ((l = fread(buf, 1, BUF_SIZE, fp)) != 0)
+ fwrite(buf, 1, l, stdout);
+ } while (l == BUF_SIZE);
+ fclose(fp);
+ kclose(x);
+ return 0;
+}
+#endif
diff --git a/vere/ext/h2o/patches/h2o-2.2.6/deps/klib/ksw.c b/vere/ext/h2o/patches/h2o-2.2.6/deps/klib/ksw.c
new file mode 100644
index 0000000..05af32f
--- /dev/null
+++ b/vere/ext/h2o/patches/h2o-2.2.6/deps/klib/ksw.c
@@ -0,0 +1,637 @@
+/* The MIT License
+
+ Copyright (c) 2011 by Attractive Chaos <attractor@live.co.uk>
+
+ Permission is hereby granted, free of charge, to any person obtaining
+ a copy of this software and associated documentation files (the
+ "Software"), to deal in the Software without restriction, including
+ without limitation the rights to use, copy, modify, merge, publish,
+ distribute, sublicense, and/or sell copies of the Software, and to
+ permit persons to whom the Software is furnished to do so, subject to
+ the following conditions:
+
+ The above copyright notice and this permission notice shall be
+ included in all copies or substantial portions of the Software.
+
+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+ EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
+ BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+ ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+ CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ SOFTWARE.
+*/
+
+#include <stdlib.h>
+#include <stdint.h>
+#if defined(URBIT_RUNTIME_CPU_AARCH64)
+#include "sse2neon.h"
+#else
+#include <emmintrin.h>
+#endif
+#include "ksw.h"
+
+#ifdef __GNUC__
+#define LIKELY(x) __builtin_expect((x),1)
+#define UNLIKELY(x) __builtin_expect((x),0)
+#else
+#define LIKELY(x) (x)
+#define UNLIKELY(x) (x)
+#endif
+
+const kswr_t g_defr = { 0, -1, -1, -1, -1, -1, -1 };
+
+struct _kswq_t {
+ int qlen, slen;
+ uint8_t shift, mdiff, max, size;
+ __m128i *qp, *H0, *H1, *E, *Hmax;
+};
+
+/**
+ * Initialize the query data structure
+ *
+ * @param size Number of bytes used to store a score; valid valures are 1 or 2
+ * @param qlen Length of the query sequence
+ * @param query Query sequence
+ * @param m Size of the alphabet
+ * @param mat Scoring matrix in a one-dimension array
+ *
+ * @return Query data structure
+ */
+kswq_t *ksw_qinit(int size, int qlen, const uint8_t *query, int m, const int8_t *mat)
+{
+ kswq_t *q;
+ int slen, a, tmp, p;
+
+ size = size > 1? 2 : 1;
+ p = 8 * (3 - size); // # values per __m128i
+ slen = (qlen + p - 1) / p; // segmented length
+ q = (kswq_t*)malloc(sizeof(kswq_t) + 256 + 16 * slen * (m + 4)); // a single block of memory
+ q->qp = (__m128i*)(((size_t)q + sizeof(kswq_t) + 15) >> 4 << 4); // align memory
+ q->H0 = q->qp + slen * m;
+ q->H1 = q->H0 + slen;
+ q->E = q->H1 + slen;
+ q->Hmax = q->E + slen;
+ q->slen = slen; q->qlen = qlen; q->size = size;
+ // compute shift
+ tmp = m * m;
+ for (a = 0, q->shift = 127, q->mdiff = 0; a < tmp; ++a) { // find the minimum and maximum score
+ if (mat[a] < (int8_t)q->shift) q->shift = mat[a];
+ if (mat[a] > (int8_t)q->mdiff) q->mdiff = mat[a];
+ }
+ q->max = q->mdiff;
+ q->shift = 256 - q->shift; // NB: q->shift is uint8_t
+ q->mdiff += q->shift; // this is the difference between the min and max scores
+ // An example: p=8, qlen=19, slen=3 and segmentation:
+ // {{0,3,6,9,12,15,18,-1},{1,4,7,10,13,16,-1,-1},{2,5,8,11,14,17,-1,-1}}
+ if (size == 1) {
+ int8_t *t = (int8_t*)q->qp;
+ for (a = 0; a < m; ++a) {
+ int i, k, nlen = slen * p;
+ const int8_t *ma = mat + a * m;
+ for (i = 0; i < slen; ++i)
+ for (k = i; k < nlen; k += slen) // p iterations
+ *t++ = (k >= qlen? 0 : ma[query[k]]) + q->shift;
+ }
+ } else {
+ int16_t *t = (int16_t*)q->qp;
+ for (a = 0; a < m; ++a) {
+ int i, k, nlen = slen * p;
+ const int8_t *ma = mat + a * m;
+ for (i = 0; i < slen; ++i)
+ for (k = i; k < nlen; k += slen) // p iterations
+ *t++ = (k >= qlen? 0 : ma[query[k]]);
+ }
+ }
+ return q;
+}
+
+kswr_t ksw_u8(kswq_t *q, int tlen, const uint8_t *target, int _gapo, int _gape, int xtra) // the first gap costs -(_o+_e)
+{
+ int slen, i, m_b, n_b, te = -1, gmax = 0, minsc, endsc;
+ uint64_t *b;
+ __m128i zero, gapoe, gape, shift, *H0, *H1, *E, *Hmax;
+ kswr_t r;
+
+#define __max_16(ret, xx) do { \
+ (xx) = _mm_max_epu8((xx), _mm_srli_si128((xx), 8)); \
+ (xx) = _mm_max_epu8((xx), _mm_srli_si128((xx), 4)); \
+ (xx) = _mm_max_epu8((xx), _mm_srli_si128((xx), 2)); \
+ (xx) = _mm_max_epu8((xx), _mm_srli_si128((xx), 1)); \
+ (ret) = _mm_extract_epi16((xx), 0) & 0x00ff; \
+ } while (0)
+
+ // initialization
+ r = g_defr;
+ minsc = (xtra&KSW_XSUBO)? xtra&0xffff : 0x10000;
+ endsc = (xtra&KSW_XSTOP)? xtra&0xffff : 0x10000;
+ m_b = n_b = 0; b = 0;
+ zero = _mm_set1_epi32(0);
+ gapoe = _mm_set1_epi8(_gapo + _gape);
+ gape = _mm_set1_epi8(_gape);
+ shift = _mm_set1_epi8(q->shift);
+ H0 = q->H0; H1 = q->H1; E = q->E; Hmax = q->Hmax;
+ slen = q->slen;
+ for (i = 0; i < slen; ++i) {
+ _mm_store_si128(E + i, zero);
+ _mm_store_si128(H0 + i, zero);
+ _mm_store_si128(Hmax + i, zero);
+ }
+ // the core loop
+ for (i = 0; i < tlen; ++i) {
+ int j, k, cmp, imax;
+ __m128i e, h, f = zero, max = zero, *S = q->qp + target[i] * slen; // s is the 1st score vector
+ h = _mm_load_si128(H0 + slen - 1); // h={2,5,8,11,14,17,-1,-1} in the above example
+ h = _mm_slli_si128(h, 1); // h=H(i-1,-1); << instead of >> because x64 is little-endian
+ for (j = 0; LIKELY(j < slen); ++j) {
+ /* SW cells are computed in the following order:
+ * H(i,j) = max{H(i-1,j-1)+S(i,j), E(i,j), F(i,j)}
+ * E(i+1,j) = max{H(i,j)-q, E(i,j)-r}
+ * F(i,j+1) = max{H(i,j)-q, F(i,j)-r}
+ */
+ // compute H'(i,j); note that at the beginning, h=H'(i-1,j-1)
+ h = _mm_adds_epu8(h, _mm_load_si128(S + j));
+ h = _mm_subs_epu8(h, shift); // h=H'(i-1,j-1)+S(i,j)
+ e = _mm_load_si128(E + j); // e=E'(i,j)
+ h = _mm_max_epu8(h, e);
+ h = _mm_max_epu8(h, f); // h=H'(i,j)
+ max = _mm_max_epu8(max, h); // set max
+ _mm_store_si128(H1 + j, h); // save to H'(i,j)
+ // now compute E'(i+1,j)
+ h = _mm_subs_epu8(h, gapoe); // h=H'(i,j)-gapo
+ e = _mm_subs_epu8(e, gape); // e=E'(i,j)-gape
+ e = _mm_max_epu8(e, h); // e=E'(i+1,j)
+ _mm_store_si128(E + j, e); // save to E'(i+1,j)
+ // now compute F'(i,j+1)
+ f = _mm_subs_epu8(f, gape);
+ f = _mm_max_epu8(f, h);
+ // get H'(i-1,j) and prepare for the next j
+ h = _mm_load_si128(H0 + j); // h=H'(i-1,j)
+ }
+ // NB: we do not need to set E(i,j) as we disallow adjecent insertion and then deletion
+ for (k = 0; LIKELY(k < 16); ++k) { // this block mimics SWPS3; NB: H(i,j) updated in the lazy-F loop cannot exceed max
+ f = _mm_slli_si128(f, 1);
+ for (j = 0; LIKELY(j < slen); ++j) {
+ h = _mm_load_si128(H1 + j);
+ h = _mm_max_epu8(h, f); // h=H'(i,j)
+ _mm_store_si128(H1 + j, h);
+ h = _mm_subs_epu8(h, gapoe);
+ f = _mm_subs_epu8(f, gape);
+ cmp = _mm_movemask_epi8(_mm_cmpeq_epi8(_mm_subs_epu8(f, h), zero));
+ if (UNLIKELY(cmp == 0xffff)) goto end_loop16;
+ }
+ }
+end_loop16:
+ //int k;for (k=0;k<16;++k)printf("%d ", ((uint8_t*)&max)[k]);printf("\n");
+ __max_16(imax, max); // imax is the maximum number in max
+ if (imax >= minsc) { // write the b array; this condition adds branching unfornately
+ if (n_b == 0 || (int32_t)b[n_b-1] + 1 != i) { // then append
+ if (n_b == m_b) {
+ m_b = m_b? m_b<<1 : 8;
+ b = (uint64_t*)realloc(b, 8 * m_b);
+ }
+ b[n_b++] = (uint64_t)imax<<32 | i;
+ } else if ((int)(b[n_b-1]>>32) < imax) b[n_b-1] = (uint64_t)imax<<32 | i; // modify the last
+ }
+ if (imax > gmax) {
+ gmax = imax; te = i; // te is the end position on the target
+ for (j = 0; LIKELY(j < slen); ++j) // keep the H1 vector
+ _mm_store_si128(Hmax + j, _mm_load_si128(H1 + j));
+ if (gmax + q->shift >= 255 || gmax >= endsc) break;
+ }
+ S = H1; H1 = H0; H0 = S; // swap H0 and H1
+ }
+ r.score = gmax + q->shift < 255? gmax : 255;
+ r.te = te;
+ if (r.score != 255) { // get a->qe, the end of query match; find the 2nd best score
+ int max = -1, low, high, qlen = slen * 16;
+ uint8_t *t = (uint8_t*)Hmax;
+ for (i = 0; i < qlen; ++i, ++t)
+ if ((int)*t > max) max = *t, r.qe = i / 16 + i % 16 * slen;
+ //printf("%d,%d\n", max, gmax);
+ if (b) {
+ i = (r.score + q->max - 1) / q->max;
+ low = te - i; high = te + i;
+ for (i = 0; i < n_b; ++i) {
+ int e = (int32_t)b[i];
+ if ((e < low || e > high) && (int)(b[i]>>32) > r.score2)
+ r.score2 = b[i]>>32, r.te2 = e;
+ }
+ }
+ }
+ free(b);
+ return r;
+}
+
+kswr_t ksw_i16(kswq_t *q, int tlen, const uint8_t *target, int _gapo, int _gape, int xtra) // the first gap costs -(_o+_e)
+{
+ int slen, i, m_b, n_b, te = -1, gmax = 0, minsc, endsc;
+ uint64_t *b;
+ __m128i zero, gapoe, gape, *H0, *H1, *E, *Hmax;
+ kswr_t r;
+
+#define __max_8(ret, xx) do { \
+ (xx) = _mm_max_epi16((xx), _mm_srli_si128((xx), 8)); \
+ (xx) = _mm_max_epi16((xx), _mm_srli_si128((xx), 4)); \
+ (xx) = _mm_max_epi16((xx), _mm_srli_si128((xx), 2)); \
+ (ret) = _mm_extract_epi16((xx), 0); \
+ } while (0)
+
+ // initialization
+ r = g_defr;
+ minsc = (xtra&KSW_XSUBO)? xtra&0xffff : 0x10000;
+ endsc = (xtra&KSW_XSTOP)? xtra&0xffff : 0x10000;
+ m_b = n_b = 0; b = 0;
+ zero = _mm_set1_epi32(0);
+ gapoe = _mm_set1_epi16(_gapo + _gape);
+ gape = _mm_set1_epi16(_gape);
+ H0 = q->H0; H1 = q->H1; E = q->E; Hmax = q->Hmax;
+ slen = q->slen;
+ for (i = 0; i < slen; ++i) {
+ _mm_store_si128(E + i, zero);
+ _mm_store_si128(H0 + i, zero);
+ _mm_store_si128(Hmax + i, zero);
+ }
+ // the core loop
+ for (i = 0; i < tlen; ++i) {
+ int j, k, imax;
+ __m128i e, h, f = zero, max = zero, *S = q->qp + target[i] * slen; // s is the 1st score vector
+ h = _mm_load_si128(H0 + slen - 1); // h={2,5,8,11,14,17,-1,-1} in the above example
+ h = _mm_slli_si128(h, 2);
+ for (j = 0; LIKELY(j < slen); ++j) {
+ h = _mm_adds_epi16(h, *S++);
+ e = _mm_load_si128(E + j);
+ h = _mm_max_epi16(h, e);
+ h = _mm_max_epi16(h, f);
+ max = _mm_max_epi16(max, h);
+ _mm_store_si128(H1 + j, h);
+ h = _mm_subs_epu16(h, gapoe);
+ e = _mm_subs_epu16(e, gape);
+ e = _mm_max_epi16(e, h);
+ _mm_store_si128(E + j, e);
+ f = _mm_subs_epu16(f, gape);
+ f = _mm_max_epi16(f, h);
+ h = _mm_load_si128(H0 + j);
+ }
+ for (k = 0; LIKELY(k < 16); ++k) {
+ f = _mm_slli_si128(f, 2);
+ for (j = 0; LIKELY(j < slen); ++j) {
+ h = _mm_load_si128(H1 + j);
+ h = _mm_max_epi16(h, f);
+ _mm_store_si128(H1 + j, h);
+ h = _mm_subs_epu16(h, gapoe);
+ f = _mm_subs_epu16(f, gape);
+ if(UNLIKELY(!_mm_movemask_epi8(_mm_cmpgt_epi16(f, h)))) goto end_loop8;
+ }
+ }
+end_loop8:
+ __max_8(imax, max);
+ if (imax >= minsc) {
+ if (n_b == 0 || (int32_t)b[n_b-1] + 1 != i) {
+ if (n_b == m_b) {
+ m_b = m_b? m_b<<1 : 8;
+ b = (uint64_t*)realloc(b, 8 * m_b);
+ }
+ b[n_b++] = (uint64_t)imax<<32 | i;
+ } else if ((int)(b[n_b-1]>>32) < imax) b[n_b-1] = (uint64_t)imax<<32 | i; // modify the last
+ }
+ if (imax > gmax) {
+ gmax = imax; te = i;
+ for (j = 0; LIKELY(j < slen); ++j)
+ _mm_store_si128(Hmax + j, _mm_load_si128(H1 + j));
+ if (gmax >= endsc) break;
+ }
+ S = H1; H1 = H0; H0 = S;
+ }
+ r.score = gmax; r.te = te;
+ {
+ int max = -1, low, high, qlen = slen * 8;
+ uint16_t *t = (uint16_t*)Hmax;
+ for (i = 0, r.qe = -1; i < qlen; ++i, ++t)
+ if ((int)*t > max) max = *t, r.qe = i / 8 + i % 8 * slen;
+ if (b) {
+ i = (r.score + q->max - 1) / q->max;
+ low = te - i; high = te + i;
+ for (i = 0; i < n_b; ++i) {
+ int e = (int32_t)b[i];
+ if ((e < low || e > high) && (int)(b[i]>>32) > r.score2)
+ r.score2 = b[i]>>32, r.te2 = e;
+ }
+ }
+ }
+ free(b);
+ return r;
+}
+
+static void revseq(int l, uint8_t *s)
+{
+ int i, t;
+ for (i = 0; i < l>>1; ++i)
+ t = s[i], s[i] = s[l - 1 - i], s[l - 1 - i] = t;
+}
+
+kswr_t ksw_align(int qlen, uint8_t *query, int tlen, uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int xtra, kswq_t **qry)
+{
+ int size;
+ kswq_t *q;
+ kswr_t r, rr;
+ kswr_t (*func)(kswq_t*, int, const uint8_t*, int, int, int);
+
+ q = (qry && *qry)? *qry : ksw_qinit((xtra&KSW_XBYTE)? 1 : 2, qlen, query, m, mat);
+ if (qry && *qry == 0) *qry = q;
+ func = q->size == 2? ksw_i16 : ksw_u8;
+ size = q->size;
+ r = func(q, tlen, target, gapo, gape, xtra);
+ if (qry == 0) free(q);
+ if ((xtra&KSW_XSTART) == 0 || ((xtra&KSW_XSUBO) && r.score < (xtra&0xffff))) return r;
+ revseq(r.qe + 1, query); revseq(r.te + 1, target); // +1 because qe/te points to the exact end, not the position after the end
+ q = ksw_qinit(size, r.qe + 1, query, m, mat);
+ rr = func(q, tlen, target, gapo, gape, KSW_XSTOP | r.score);
+ revseq(r.qe + 1, query); revseq(r.te + 1, target);
+ free(q);
+ if (r.score == rr.score)
+ r.tb = r.te - rr.te, r.qb = r.qe - rr.qe;
+ return r;
+}
+
+/********************
+ *** SW extension ***
+ ********************/
+
+typedef struct {
+ int32_t h, e;
+} eh_t;
+
+int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int w, int h0, int *_qle, int *_tle)
+{
+ eh_t *eh; // score array
+ int8_t *qp; // query profile
+ int i, j, k, gapoe = gapo + gape, beg, end, max, max_i, max_j, max_gap;
+ if (h0 < 0) h0 = 0;
+ // allocate memory
+ qp = malloc(qlen * m);
+ eh = calloc(qlen + 1, 8);
+ // generate the query profile
+ for (k = i = 0; k < m; ++k) {
+ const int8_t *p = &mat[k * m];
+ for (j = 0; j < qlen; ++j) qp[i++] = p[query[j]];
+ }
+ // fill the first row
+ eh[0].h = h0; eh[1].h = h0 > gapoe? h0 - gapoe : 0;
+ for (j = 2; j <= qlen && eh[j-1].h > gape; ++j)
+ eh[j].h = eh[j-1].h - gape;
+ // adjust $w if it is too large
+ k = m * m;
+ for (i = 0, max = 0; i < k; ++i) // get the max score
+ max = max > mat[i]? max : mat[i];
+ max_gap = (int)((double)(qlen * max - gapo) / gape + 1.);
+ max_gap = max_gap > 1? max_gap : 1;
+ w = w < max_gap? w : max_gap;
+ // DP loop
+ max = h0, max_i = max_j = -1;
+ beg = 0, end = qlen;
+ for (i = 0; LIKELY(i < tlen); ++i) {
+ int f = 0, h1, m = 0, mj = -1;
+ int8_t *q = &qp[target[i] * qlen];
+ // compute the first column
+ h1 = h0 - (gapo + gape * (i + 1));
+ if (h1 < 0) h1 = 0;
+ // apply the band and the constraint (if provided)
+ if (beg < i - w) beg = i - w;
+ if (end > i + w + 1) end = i + w + 1;
+ if (end > qlen) end = qlen;
+ for (j = beg; LIKELY(j < end); ++j) {
+ // At the beginning of the loop: eh[j] = { H(i-1,j-1), E(i,j) }, f = F(i,j) and h1 = H(i,j-1)
+ // Similar to SSE2-SW, cells are computed in the following order:
+ // H(i,j) = max{H(i-1,j-1)+S(i,j), E(i,j), F(i,j)}
+ // E(i+1,j) = max{H(i,j)-gapo, E(i,j)} - gape
+ // F(i,j+1) = max{H(i,j)-gapo, F(i,j)} - gape
+ eh_t *p = &eh[j];
+ int h = p->h, e = p->e; // get H(i-1,j-1) and E(i-1,j)
+ p->h = h1; // set H(i,j-1) for the next row
+ h += q[j];
+ h = h > e? h : e;
+ h = h > f? h : f;
+ h1 = h; // save H(i,j) to h1 for the next column
+ mj = m > h? mj : j;
+ m = m > h? m : h; // m is stored at eh[mj+1]
+ h -= gapoe;
+ h = h > 0? h : 0;
+ e -= gape;
+ e = e > h? e : h; // computed E(i+1,j)
+ p->e = e; // save E(i+1,j) for the next row
+ f -= gape;
+ f = f > h? f : h; // computed F(i,j+1)
+ }
+ eh[end].h = h1; eh[end].e = 0;
+ if (m == 0) break;
+ if (m > max) max = m, max_i = i, max_j = mj;
+ // update beg and end for the next round
+ for (j = mj; j >= beg && eh[j].h; --j);
+ beg = j + 1;
+ for (j = mj + 2; j <= end && eh[j].h; ++j);
+ end = j;
+ //beg = 0; end = qlen; // uncomment this line for debugging
+ }
+ free(eh); free(qp);
+ if (_qle) *_qle = max_j + 1;
+ if (_tle) *_tle = max_i + 1;
+ return max;
+}
+
+/********************
+ * Global alignment *
+ ********************/
+
+#define MINUS_INF -0x40000000
+
+static inline uint32_t *push_cigar(int *n_cigar, int *m_cigar, uint32_t *cigar, int op, int len)
+{
+ if (*n_cigar == 0 || op != (cigar[(*n_cigar) - 1]&0xf)) {
+ if (*n_cigar == *m_cigar) {
+ *m_cigar = *m_cigar? (*m_cigar)<<1 : 4;
+ cigar = realloc(cigar, (*m_cigar) << 2);
+ }
+ cigar[(*n_cigar)++] = len<<4 | op;
+ } else cigar[(*n_cigar)-1] += len<<4;
+ return cigar;
+}
+
+int ksw_global(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int w, int *n_cigar_, uint32_t **cigar_)
+{
+ eh_t *eh;
+ int8_t *qp; // query profile
+ int i, j, k, gapoe = gapo + gape, score, n_col;
+ uint8_t *z; // backtrack matrix; in each cell: f<<4|e<<2|h; in principle, we can halve the memory, but backtrack will be a little more complex
+ if (n_cigar_) *n_cigar_ = 0;
+ // allocate memory
+ n_col = qlen < 2*w+1? qlen : 2*w+1; // maximum #columns of the backtrack matrix
+ z = malloc(n_col * tlen);
+ qp = malloc(qlen * m);
+ eh = calloc(qlen + 1, 8);
+ // generate the query profile
+ for (k = i = 0; k < m; ++k) {
+ const int8_t *p = &mat[k * m];
+ for (j = 0; j < qlen; ++j) qp[i++] = p[query[j]];
+ }
+ // fill the first row
+ eh[0].h = 0; eh[0].e = MINUS_INF;
+ for (j = 1; j <= qlen && j <= w; ++j)
+ eh[j].h = -(gapo + gape * j), eh[j].e = MINUS_INF;
+ for (; j <= qlen; ++j) eh[j].h = eh[j].e = MINUS_INF; // everything is -inf outside the band
+ // DP loop
+ for (i = 0; LIKELY(i < tlen); ++i) { // target sequence is in the outer loop
+ int32_t f = MINUS_INF, h1, beg, end;
+ int8_t *q = &qp[target[i] * qlen];
+ uint8_t *zi = &z[i * n_col];
+ beg = i > w? i - w : 0;
+ end = i + w + 1 < qlen? i + w + 1 : qlen; // only loop through [beg,end) of the query sequence
+ h1 = beg == 0? -(gapo + gape * (i + 1)) : MINUS_INF;
+ for (j = beg; LIKELY(j < end); ++j) {
+ // This loop is organized in a similar way to ksw_extend() and ksw_sse2(), except:
+ // 1) not checking h>0; 2) recording direction for backtracking
+ eh_t *p = &eh[j];
+ int32_t h = p->h, e = p->e;
+ uint8_t d; // direction
+ p->h = h1;
+ h += q[j];
+ d = h > e? 0 : 1;
+ h = h > e? h : e;
+ d = h > f? d : 2;
+ h = h > f? h : f;
+ h1 = h;
+ h -= gapoe;
+ e -= gape;
+ d |= e > h? 1<<2 : 0;
+ e = e > h? e : h;
+ p->e = e;
+ f -= gape;
+ d |= f > h? 2<<4 : 0; // if we want to halve the memory, use one bit only, instead of two
+ f = f > h? f : h;
+ zi[j - beg] = d; // z[i,j] keeps h for the current cell and e/f for the next cell
+ }
+ eh[end].h = h1; eh[end].e = MINUS_INF;
+ }
+ score = eh[qlen].h;
+ if (n_cigar_ && cigar_) { // backtrack
+ int n_cigar = 0, m_cigar = 0, which = 0;
+ uint32_t *cigar = 0, tmp;
+ i = tlen - 1; k = (i + w + 1 < qlen? i + w + 1 : qlen) - 1; // (i,k) points to the last cell
+ while (i >= 0 && k >= 0) {
+ which = z[i * n_col + (k - (i > w? i - w : 0))] >> (which<<1) & 3;
+ if (which == 0) cigar = push_cigar(&n_cigar, &m_cigar, cigar, 0, 1), --i, --k;
+ else if (which == 1) cigar = push_cigar(&n_cigar, &m_cigar, cigar, 2, 1), --i;
+ else cigar = push_cigar(&n_cigar, &m_cigar, cigar, 1, 1), --k;
+ }
+ if (i >= 0) cigar = push_cigar(&n_cigar, &m_cigar, cigar, 2, i + 1);
+ if (k >= 0) cigar = push_cigar(&n_cigar, &m_cigar, cigar, 1, k + 1);
+ for (i = 0; i < n_cigar>>1; ++i) // reverse CIGAR
+ tmp = cigar[i], cigar[i] = cigar[n_cigar-1-i], cigar[n_cigar-1-i] = tmp;
+ *n_cigar_ = n_cigar, *cigar_ = cigar;
+ }
+ free(eh); free(qp); free(z);
+ return score;
+}
+
+/*******************************************
+ * Main function (not compiled by default) *
+ *******************************************/
+
+#ifdef _KSW_MAIN
+
+#include <unistd.h>
+#include <stdio.h>
+#include <zlib.h>
+#include "kseq.h"
+KSEQ_INIT(gzFile, gzread)
+
+unsigned char seq_nt4_table[256] = {
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4
+};
+
+int main(int argc, char *argv[])
+{
+ int c, sa = 1, sb = 3, i, j, k, forward_only = 0, max_rseq = 0;
+ int8_t mat[25];
+ int gapo = 5, gape = 2, minsc = 0, xtra = KSW_XSTART;
+ uint8_t *rseq = 0;
+ gzFile fpt, fpq;
+ kseq_t *kst, *ksq;
+
+ // parse command line
+ while ((c = getopt(argc, argv, "a:b:q:r:ft:1")) >= 0) {
+ switch (c) {
+ case 'a': sa = atoi(optarg); break;
+ case 'b': sb = atoi(optarg); break;
+ case 'q': gapo = atoi(optarg); break;
+ case 'r': gape = atoi(optarg); break;
+ case 't': minsc = atoi(optarg); break;
+ case 'f': forward_only = 1; break;
+ case '1': xtra |= KSW_XBYTE; break;
+ }
+ }
+ if (optind + 2 > argc) {
+ fprintf(stderr, "Usage: ksw [-1] [-f] [-a%d] [-b%d] [-q%d] [-r%d] [-t%d] <target.fa> <query.fa>\n", sa, sb, gapo, gape, minsc);
+ return 1;
+ }
+ if (minsc > 0xffff) minsc = 0xffff;
+ xtra |= KSW_XSUBO | minsc;
+ // initialize scoring matrix
+ for (i = k = 0; i < 4; ++i) {
+ for (j = 0; j < 4; ++j)
+ mat[k++] = i == j? sa : -sb;
+ mat[k++] = 0; // ambiguous base
+ }
+ for (j = 0; j < 5; ++j) mat[k++] = 0;
+ // open file
+ fpt = gzopen(argv[optind], "r"); kst = kseq_init(fpt);
+ fpq = gzopen(argv[optind+1], "r"); ksq = kseq_init(fpq);
+ // all-pair alignment
+ while (kseq_read(ksq) > 0) {
+ kswq_t *q[2] = {0, 0};
+ kswr_t r;
+ for (i = 0; i < (int)ksq->seq.l; ++i) ksq->seq.s[i] = seq_nt4_table[(int)ksq->seq.s[i]];
+ if (!forward_only) { // reverse
+ if ((int)ksq->seq.m > max_rseq) {
+ max_rseq = ksq->seq.m;
+ rseq = (uint8_t*)realloc(rseq, max_rseq);
+ }
+ for (i = 0, j = ksq->seq.l - 1; i < (int)ksq->seq.l; ++i, --j)
+ rseq[j] = ksq->seq.s[i] == 4? 4 : 3 - ksq->seq.s[i];
+ }
+ gzrewind(fpt); kseq_rewind(kst);
+ while (kseq_read(kst) > 0) {
+ for (i = 0; i < (int)kst->seq.l; ++i) kst->seq.s[i] = seq_nt4_table[(int)kst->seq.s[i]];
+ r = ksw_align(ksq->seq.l, (uint8_t*)ksq->seq.s, kst->seq.l, (uint8_t*)kst->seq.s, 5, mat, gapo, gape, xtra, &q[0]);
+ if (r.score >= minsc)
+ printf("%s\t%d\t%d\t%s\t%d\t%d\t%d\t%d\t%d\n", kst->name.s, r.tb, r.te+1, ksq->name.s, r.qb, r.qe+1, r.score, r.score2, r.te2);
+ if (rseq) {
+ r = ksw_align(ksq->seq.l, rseq, kst->seq.l, (uint8_t*)kst->seq.s, 5, mat, gapo, gape, xtra, &q[1]);
+ if (r.score >= minsc)
+ printf("%s\t%d\t%d\t%s\t%d\t%d\t%d\t%d\t%d\n", kst->name.s, r.tb, r.te+1, ksq->name.s, (int)ksq->seq.l - r.qb, (int)ksq->seq.l - 1 - r.qe, r.score, r.score2, r.te2);
+ }
+ }
+ free(q[0]); free(q[1]);
+ }
+ free(rseq);
+ kseq_destroy(kst); gzclose(fpt);
+ kseq_destroy(ksq); gzclose(fpq);
+ return 0;
+}
+#endif
diff --git a/vere/ext/h2o/patches/h2o-2.2.6/deps/klib/test/kbit_test.c b/vere/ext/h2o/patches/h2o-2.2.6/deps/klib/test/kbit_test.c
new file mode 100644
index 0000000..ece2a45
--- /dev/null
+++ b/vere/ext/h2o/patches/h2o-2.2.6/deps/klib/test/kbit_test.c
@@ -0,0 +1,141 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <time.h>
+#if defined(URBIT_RUNTIME_CPU_AARCH64)
+#include "sse2neon.h"
+#else
+#include <emmintrin.h>
+#endif
+#include "kbit.h"
+
+// from bowtie-0.9.8.1
+inline static int bt1_pop64(uint64_t x) // the kbi_popcount64() equivalence; similar to popcount_2() in wiki
+{
+ x -= ((x >> 1) & 0x5555555555555555llu);
+ x = (x & 0x3333333333333333llu) + ((x >> 2) & 0x3333333333333333llu);
+ x = (x + (x >> 4)) & 0x0F0F0F0F0F0F0F0Fllu;
+ x = x + (x >> 8);
+ x = x + (x >> 16);
+ x = x + (x >> 32);
+ return x & 0x3F;
+}
+
+inline static int bt1_countInU64(uint64_t dw, int c) // the kbi_DNAcount64() equivalence
+{
+ uint64_t dwA = dw & 0xAAAAAAAAAAAAAAAAllu;
+ uint64_t dwNA = dw & ~0xAAAAAAAAAAAAAAAAllu;
+ uint64_t tmp;
+ switch (c) {
+ case 0: tmp = (dwA >> 1) | dwNA; break;
+ case 1: tmp = ~(dwA >> 1) & dwNA; break;
+ case 2: tmp = (dwA >> 1) & ~dwNA; break;
+ default: tmp = (dwA >> 1) & dwNA;
+ }
+ tmp = bt1_pop64(tmp);
+ if (c == 0) tmp = 32 - tmp;
+ return (int)tmp;
+}
+
+// from bigmagic
+static uint32_t sse2_bit_count32(const __m128i* block, const __m128i* block_end)
+{
+ const unsigned mu1 = 0x55555555;
+ const unsigned mu2 = 0x33333333;
+ const unsigned mu3 = 0x0F0F0F0F;
+ const unsigned mu4 = 0x0000003F;
+
+ uint32_t tcnt[4];
+
+ // Loading masks
+ __m128i m1 = _mm_set_epi32 (mu1, mu1, mu1, mu1);
+ __m128i m2 = _mm_set_epi32 (mu2, mu2, mu2, mu2);
+ __m128i m3 = _mm_set_epi32 (mu3, mu3, mu3, mu3);
+ __m128i m4 = _mm_set_epi32 (mu4, mu4, mu4, mu4);
+ __m128i mcnt;
+ mcnt = _mm_xor_si128(m1, m1); // cnt = 0
+
+ __m128i tmp1, tmp2;
+ do
+ {
+ __m128i b = _mm_load_si128(block);
+ ++block;
+
+ // b = (b & 0x55555555) + (b >> 1 & 0x55555555);
+ tmp1 = _mm_srli_epi32(b, 1); // tmp1 = (b >> 1 & 0x55555555)
+ tmp1 = _mm_and_si128(tmp1, m1);
+ tmp2 = _mm_and_si128(b, m1); // tmp2 = (b & 0x55555555)
+ b = _mm_add_epi32(tmp1, tmp2); // b = tmp1 + tmp2
+
+ // b = (b & 0x33333333) + (b >> 2 & 0x33333333);
+ tmp1 = _mm_srli_epi32(b, 2); // (b >> 2 & 0x33333333)
+ tmp1 = _mm_and_si128(tmp1, m2);
+ tmp2 = _mm_and_si128(b, m2); // (b & 0x33333333)
+ b = _mm_add_epi32(tmp1, tmp2); // b = tmp1 + tmp2
+
+ // b = (b + (b >> 4)) & 0x0F0F0F0F;
+ tmp1 = _mm_srli_epi32(b, 4); // tmp1 = b >> 4
+ b = _mm_add_epi32(b, tmp1); // b = b + (b >> 4)
+ b = _mm_and_si128(b, m3); // & 0x0F0F0F0F
+
+ // b = b + (b >> 8);
+ tmp1 = _mm_srli_epi32 (b, 8); // tmp1 = b >> 8
+ b = _mm_add_epi32(b, tmp1); // b = b + (b >> 8)
+
+ // b = (b + (b >> 16)) & 0x0000003F;
+ tmp1 = _mm_srli_epi32 (b, 16); // b >> 16
+ b = _mm_add_epi32(b, tmp1); // b + (b >> 16)
+ b = _mm_and_si128(b, m4); // (b >> 16) & 0x0000003F;
+
+ mcnt = _mm_add_epi32(mcnt, b); // mcnt += b
+
+ } while (block < block_end);
+
+ _mm_store_si128((__m128i*)tcnt, mcnt);
+
+ return tcnt[0] + tcnt[1] + tcnt[2] + tcnt[3];
+}
+
+int main(void)
+{
+ int i, N = 100000000;
+ uint64_t *x, cnt;
+ clock_t t;
+ int c = 1;
+
+ x = (uint64_t*)calloc(N, 8);
+ srand48(11);
+ for (i = 0; i < N; ++i)
+ x[i] = (uint64_t)lrand48() << 32 ^ lrand48();
+
+ fprintf(stderr, "\n===> Calculate # of 1 in an integer (popcount) <===\n");
+
+ t = clock(); cnt = 0;
+ for (i = 0; i < N; ++i) cnt += kbi_popcount64(x[i]);
+ fprintf(stderr, "%20s\t%20ld\t%10.6f\n", "kbit", (long)cnt, (double)(clock() - t) / CLOCKS_PER_SEC);
+
+ t = clock(); cnt = 0;
+ for (i = 0; i < N; ++i) cnt += bt1_pop64(x[i]);
+ fprintf(stderr, "%20s\t%20ld\t%10.6f\n", "wiki-popcount_2", (long)cnt, (double)(clock() - t) / CLOCKS_PER_SEC);
+
+ t = clock(); cnt = 0;
+ for (i = 0; i < N; ++i) cnt += __builtin_popcountl(x[i]);
+ fprintf(stderr, "%20s\t%20ld\t%10.6f\n", "__builtin_popcountl", (long)cnt, (double)(clock() - t) / CLOCKS_PER_SEC);
+
+ t = clock(); cnt = 0;
+ cnt += sse2_bit_count32((__m128i*)x, (__m128i*)(x+N));
+ fprintf(stderr, "%20s\t%20ld\t%10.6f\n", "SSE2-32bit", (long)cnt, (double)(clock() - t) / CLOCKS_PER_SEC);
+
+ fprintf(stderr, "\n===> Count '%c' in 2-bit encoded integers <===\n", "ACGT"[c]);
+
+ t = clock(); cnt = 0;
+ for (i = 0; i < N; ++i) cnt += kbi_DNAcount64(x[i], c);
+ fprintf(stderr, "%20s\t%20ld\t%10.6f\n", "kbit", (long)cnt, (double)(clock() - t) / CLOCKS_PER_SEC);
+
+ t = clock(); cnt = 0;
+ for (i = 0; i < N; ++i) cnt += bt1_countInU64(x[i], c);
+ fprintf(stderr, "%20s\t%20ld\t%10.6f\n", "bowtie1", (long)cnt, (double)(clock() - t) / CLOCKS_PER_SEC);
+
+ fprintf(stderr, "\n");
+ free(x);
+ return 0;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2025-10-16 11:10:03 +0000