diff --git a/cbits/blob.c b/cbits/blob.c new file mode 100644 --- /dev/null +++ b/cbits/blob.c @@ -0,0 +1,275 @@ + +#include +#include + +// ----------------------------------------------------------------------------- + +void identity(int n, const uint64_t *src, int* pm, uint64_t *tgt) { + memcpy(tgt, src, n<<3); + *pm = n; +} + +void cons(uint64_t x, int n, const uint64_t *src, int* pm, uint64_t *tgt) { + memcpy(tgt+1, src, n<<3); + tgt[0] = x; + *pm = n+1; +} + +void snoc(uint64_t x, int n, const uint64_t *src, int* pm, uint64_t *tgt) { + memcpy(tgt, src, n<<3); + tgt[n] = x; + *pm = n+1; +} + +void tail(int n, const uint64_t *src, int* pm, uint64_t *tgt) { + if (n>1) { + memcpy(tgt, src+1, (n-1)<<3); + *pm = n-1; + } + else { + tgt[0] = 0; + *pm = 1; + } +} + +// ----------------------------------------------------------------------------- + +#define IDX_RIGHT(i,shift) ( ( (i) + (shift) ) % n ) +#define IDX_LEFT( i,shift) ( ( (i) + n - (shift) ) % n ) + +void rotate_left_words(int k0, int n, const uint64_t *src, int* pm, uint64_t *tgt) +{ + int k = k0 % n; + for(int i=0;i> r); + x = src[i]; + } + *pm = n; +} + +// we assume that 0 <= k < 64 +void rotate_right_bits(int k, int n, const uint64_t *src, int* pm, uint64_t *tgt) +{ + if (k==0) { identity(n,src,pm,tgt); return; } + + int r = 64 - k; + uint64_t x = src[n-1]; + for(int i=0;i> k) | (src[i+1] << r); + } + tgt[n-1] = (src[n-1] >> k) | (src[0] << r); + *pm = n; +} + +void rotate_left(int k0, int n, const uint64_t *src, int* pm, uint64_t *tgt) +{ + if (k0==0) { identity(n,src,pm,tgt); return; } + + int k = k0 & 63; + int s = (k0 >> 6) % n; + int r = 64 - k; + uint64_t x = src[ IDX_LEFT(n-1,s) ]; + + if (k==0) { rotate_left_words(s,n,src,pm,tgt); return; } + if (s==0) { rotate_left_bits (k,n,src,pm,tgt); return; } + + for(int i=0;i> r); + x = src[o]; + } + *pm = n; +} + +void rotate_right(int k0, int n, const uint64_t *src, int* pm, uint64_t *tgt) +{ + if (k0==0) { identity(n,src,pm,tgt); return; } + + int k = k0 & 63; + int s = (k0 >> 6) % n; + int r = 64 - k; + uint64_t x = src[ IDX_RIGHT(n-1,s) ]; + + if (k==0) { rotate_right_words(s,n,src,pm,tgt); return; } + if (s==0) { rotate_right_bits (k,n,src,pm,tgt); return; } + + for(int i=0;i> k) + | ( src[ IDX_RIGHT(i+1,s)] << r); + } + *pm = n; +} + +// ----------------------------------------------------------------------------- + +void shift_left_words(int k, int n, const uint64_t *src, int* pm, uint64_t *tgt) +{ + for(int i=0;i= n) { tgt[0]=0; *pm=1; return; } + + for(int i=0;i> r); + x = src[i]; + } + // always increase the size + tgt[n] = (x >> r); + *pm = n+1; +} + +// we assume that 0 <= k < 64 +// non-strict version: we adopt the convention here that blobs are +// extended with zeros to infinity, and only increase the size if necessary +void shift_left_bits_nonstrict(int k, int n, const uint64_t *src, int* pm, uint64_t *tgt) +{ + if (k==0) { identity(n,src,pm,tgt); return; } + + int r = 64 - k; + uint64_t x = 0; + for(int i=0;i> r); + x = src[i]; + } + // we adopt the convention here that blobs are extended with zeros to infinity + // and only increase the size if necessary + uint64_t y = (x >> r); + if (y==0) { + *pm = n; + } + else { + tgt[n] = y; + *pm = n+1; + } + +// tgt[n] = (x >> r); +// *pm = n+1; +} + +// we assume that 0 <= k < 64 +void shift_right_bits(int k, int n, const uint64_t *src, int* pm, uint64_t *tgt) +{ + if (k==0) { identity(n,src,pm,tgt); return; } + + int r = 64 - k; + for(int i=0;i> k) | (src[i+1] << r); + } + tgt[n-1] = src[n-1] >> k; + *pm = n; +} + +// strict version: we _always_ increase the size +void shift_left_strict(int k0, int n, const uint64_t *src, int* pm, uint64_t *tgt) +{ + if (k0==0) { identity(n,src,pm,tgt); return; } + + int k = k0 & 63; + int s = (k0 >> 6); + int r = 64 - k; + + if (k==0) { shift_left_words(s,n,src,pm,tgt); return; } + if (s==0) { shift_left_bits_strict(k,n,src,pm,tgt); return; } + + for(int i=0;i> r); + x = src[i]; + } + // we always always increase + tgt[n+s] = (x >> r); + *pm = n+s+1; +} + +// non-strict version: +// we adopt the convention here that blobs are extended with +// zeros to infinity and only increase the size if necessary +void shift_left_nonstrict(int k0, int n, const uint64_t *src, int* pm, uint64_t *tgt) +{ + if (k0==0) { identity(n,src,pm,tgt); return; } + + int k = k0 & 63; + int s = (k0 >> 6); + int r = 64 - k; + + if (k==0) { shift_left_words (s,n,src,pm,tgt); return; } + if (s==0) { shift_left_bits_nonstrict (k,n,src,pm,tgt); return; } + + for(int i=0;i> r); + x = src[i]; + } + // we adopt the convention here that blobs are extended with zeros to infinity + // and only increase the size if necessary + uint64_t y = (x >> r); + if (y==0) { + *pm = n+s; + } + else { + tgt[n+s] = y; + *pm = n+s+1; + } +} + +void shift_right(int k0, int n, const uint64_t *src, int* pm, uint64_t *tgt) +{ + if (k0==0) { identity(n,src,pm,tgt); return; } + + int k = k0 & 63; + int s = (k0 >> 6); + int r = 64 - k; + + if (k==0) { shift_right_words(s,n,src,pm,tgt); return; } + if (s==0) { shift_right_bits (k,n,src,pm,tgt); return; } + + for(int i=0;i> k) | (src[i+s+1] << r); + } + tgt[n-s-1] = src[n-1] >> k; + *pm = n-s; +} + +// ----------------------------------------------------------------------------- diff --git a/cbits/wordvec.c b/cbits/wordvec.c new file mode 100644 --- /dev/null +++ b/cbits/wordvec.c @@ -0,0 +1,884 @@ + +#include +#include + +#include // DEBUGGING ONLY + +#include "blob.h" + +// ----------------------------------------------------------------------------- + +#define LT (-1) +#define EQ 0 +#define GT 1 + +// ----------------------------------------------------------------------------- + +void vec_identity(int n, const uint64_t *src, int* pm, uint64_t *tgt) { + memcpy(tgt, src, n<<3); + *pm = n; +} + +// ----------------------------------------------------------------------------- + +// ones for the first n bit, zeros for the rest +inline uint64_t nbit_mask (int n) { + uint64_t mask = 1; // fucking c implicit conversions + if (n<64) { + mask = (mask << n) - 1; + } + else { + mask = 0; + mask = ~mask; // 0xfffff...f + } + return mask; +} + +// zeros for the first n bit, ones for the rest +inline uint64_t nbit_compl_mask(int n) { + uint64_t mask = 1; + if (n<64) { + mask = (mask << n) - 1; + return ~mask; + } + return 0; +} + +// the minimum required bits to a store a given number, rounded up to multiples of 4 +inline int required_bits_not_rounded(uint64_t x) +{ + int bits = 0; + while(x > 0) { x = (x>>1); bits++; } + if (bits == 0) { bits = 1; } + return bits; +} + +// the minimum required bits to a store a given number, rounded up to multiples of 4 +inline int required_bits(uint64_t x) +{ + int bits = 0; + while(x > 0) { x = (x>>1); bits++; } + if (bits == 0) { bits = 1; } + bits = (bits+3) & (~3); + return bits; +} + +int export_required_bits_not_rounded(uint64_t x) { return required_bits_not_rounded(x); } +int export_required_bits (uint64_t x) { return required_bits (x); } + +inline int bits2reso(int bits) +{ + return ( (bits >> 2) - 1 ); +} + +inline int required_reso(uint64_t x) +{ + return ( (required_bits(x) >> 2) - 1 ); +} + +// ----------------------------------------------------------------------------- + +#define MAX_SMALL_LENGTH 31 +#define MAX_SMALL_BITS 16 + +// header of the empty vector (small, 4 bits, 0 length) +#define EMPTY_HEADER 0 + +#define SMALL_HEADER(len,reso) ( ((reso) << 1) | (((uint64_t)(len)) << 3) ) +#define BIG_HEADER(len,reso) ( 1 | ((reso) << 1) | (((uint64_t)(len)) << 5) ) + +#define VEC_HEADER_CODE(src) \ + uint64_t head = src[0]; \ + int is_small = (head & 1) ^ 1; \ + \ + int header_bits, reso_bits, len_bits, len_ofs; \ + uint64_t reso_mask, len_mask, header_mask; \ + \ + if (is_small) \ + { \ + header_bits = 8; header_mask = 0xff; \ + reso_bits = 2; reso_mask = 0x03; \ + len_bits = 5; len_mask = 0x1f; \ + } \ + else \ + { \ + header_bits = 32; header_mask = 0xffffffff; \ + reso_bits = 4; reso_mask = 0x0f; \ + len_bits = 27; len_mask = 0x07ffffff; \ + } \ + \ + int reso,bits,len; \ + reso = (head >> 1) & reso_mask; \ + bits = ((reso + 1) << 2); \ + len = (head >> (1 + reso_bits)) & len_mask; + +// ------------------------------------- + +#define VEC_READ_LOOP \ + const uint64_t *p = src; \ + int p_ofs = header_bits; \ + uint64_t elem_mask = nbit_mask(bits); \ + for(int i=0;i> p_ofs); \ + } \ + else { \ + elem = ( p[0] >> p_ofs ) \ + | ( (p[1] & nbit_mask(p_new-64)) << (64-p_ofs)); \ + } \ + elem &= elem_mask; \ + if (p_new >= 64) { \ + p_ofs = p_new-64; \ + p++; \ + } \ + else { \ + p_ofs = p_new; \ + } + +// zipping, extended with zeros, length given by the use (zip_len) +#define VEC_ZIP_LOOP \ + const uint64_t *p1 = src1; \ + const uint64_t *p2 = src2; \ + int p_ofs1 = header_bits1; \ + int p_ofs2 = header_bits2; \ + uint64_t elem_mask1 = nbit_mask(bits1); \ + uint64_t elem_mask2 = nbit_mask(bits2); \ + for(int i=0;i> p_ofs1); \ + } \ + else { \ + elem1 = ( p1[0] >> p_ofs1 ) \ + | ( (p1[1] & nbit_mask(p_new1-64)) << (64-p_ofs1)); \ + } \ + elem1 &= elem_mask1; \ + if (p_new1 >= 64) { \ + p_ofs1 = p_new1-64; \ + p1++; \ + } \ + else { \ + p_ofs1 = p_new1; \ + } \ + } \ + if (i> p_ofs2); \ + } \ + else { \ + elem2 = ( p2[0] >> p_ofs2 ) \ + | ( (p2[1] & nbit_mask(p_new2-64)) << (64-p_ofs2)); \ + } \ + elem2 &= elem_mask2; \ + if (p_new2 >= 64) { \ + p_ofs2 = p_new2-64; \ + p2++; \ + } \ + else { \ + p_ofs2 = p_new2; \ + } \ + } + +// write next element +#define WRITE_ELEMENT(elem) \ + int q_new = q_ofs + tgt_bits; \ + if (q_new <= 64) { \ + uint64_t tmp; \ + tmp = q[0] & nbit_mask(q_ofs); \ + q[0] = tmp | (elem << q_ofs); \ + } \ + else { \ + uint64_t tmp; \ + tmp = q[0] & nbit_mask(q_ofs); \ + q[0] = tmp | (elem << q_ofs); \ + q[1] = elem >> (64-q_ofs); \ + } \ + if (q_new >= 64) { \ + q_ofs = q_new-64; \ + q++; \ + } \ + else { \ + q_ofs = q_new; \ + } + + +#define STORE_OUTPUT_LENGTH(tgt_len) \ + if (q_ofs == 0) { \ + *tgt_len = (q - tgt); \ + } \ + else { \ + *tgt_len = (q - tgt + 1); \ + } + +// ----------------------------------------------------------------------------- + +void copy_elements_into + ( int src_len , int src_bits , const uint64_t *src , int src_bit_ofs + , int *tgt_len , int tgt_bits , uint64_t *tgt , int tgt_bit_ofs + ) +{ + const uint64_t *p = src; + uint64_t *q = tgt; + + int p_ofs = src_bit_ofs; + int q_ofs = tgt_bit_ofs; + + p += (p_ofs >> 6); p_ofs &= 63; + q += (q_ofs >> 6); q_ofs &= 63; + + uint64_t elem_mask = nbit_mask(src_bits); + + for(int i=0;i> p_ofs); + } + else { + elem = ( p[0] >> p_ofs ) + | ( (p[1] & nbit_mask(p_new-64)) << (64-p_ofs)); + } + elem &= elem_mask; + if (p_new >= 64) { + p_ofs = p_new-64; + p++; + } + else { + p_ofs = p_new; + } + + // write next element + int q_new = q_ofs + tgt_bits; + if (q_new <= 64) { + tmp = q[0] & nbit_mask(q_ofs); + q[0] = tmp | (elem << q_ofs); + } + else { + tmp = q[0] & nbit_mask(q_ofs); + q[0] = tmp | (elem << q_ofs); + q[1] = elem >> (64-q_ofs); + } + if (q_new >= 64) { + q_ofs = q_new-64; + q++; + } + else { + q_ofs = q_new; + } + } + + if (q_ofs == 0) { + *tgt_len = (q - tgt); + } + else { + *tgt_len = (q - tgt + 1); + } + +} + +// ----------------------------------------------------------------------------- + +void vec_tail(int n, const uint64_t *src, int* pm, uint64_t *tgt) +{ + VEC_HEADER_CODE(src) + + if (len==0) { tgt[0] = EMPTY_HEADER; *pm = 1; return; } + + if (is_small) { + shift_right(bits, n, src, pm, tgt); + tgt[0] = (tgt[0] & (~header_mask)) | SMALL_HEADER(len-1,reso); + } + else { + shift_right(bits, n, src, pm, tgt); + tgt[0] = (tgt[0] & (~header_mask)) | BIG_HEADER(len-1,reso); + } +} + +uint64_t vec_head_tail(int n, const uint64_t *src, int* pm, uint64_t *tgt) +{ + VEC_HEADER_CODE(src) + + if (len==0) { tgt[0] = EMPTY_HEADER; *pm = 1; return 0; } + + if (is_small) { + uint64_t head; + head = (src[0] >> 8) & nbit_mask(bits); + shift_right(bits, n, src, pm, tgt); + tgt[0] = (tgt[0] & (~header_mask)) | SMALL_HEADER(len-1,reso); + return head; + } + else { + uint64_t head; + if (bits <= 32) { + head = (src[0] >> 32) & nbit_mask(bits); + } + else { + head = ((src[0] >> 32) | (src[1] << 32)) & nbit_mask(bits); + } + shift_right(bits, n, src, pm, tgt); + tgt[0] = (tgt[0] & (~header_mask)) | BIG_HEADER(len-1,reso); + return head; + } +} + +// ----------------------------------------------------------------------------- +// CONS + +void vec_cons(uint64_t x, int n, const uint64_t *src, int* pm, uint64_t *tgt) +{ + VEC_HEADER_CODE(src) + + int x_bits = required_bits(x); + int x_reso = required_reso(x); + + if (len==0) { + // cons to an empty vector + if (x_reso <= 2) { + tgt[0] = SMALL_HEADER(1,x_reso) | (x << 8); + *pm = 1; + return; + } + else { + tgt[0] = BIG_HEADER(1,x_reso) | (x << 32); + uint64_t y = (x >> 32); + if (y==0) { + *pm = 1; + } + else { + tgt[1] = y; + *pm = 2; + } + return; + } + } + + if (x_bits <= bits) { + // the new element fits without changing the resolution + if (is_small) { + // the old vector is small + if (len+1 <= MAX_SMALL_LENGTH) { + // the length fits, too + uint64_t mask = nbit_compl_mask(8+bits); + shift_left_strict(bits, n, src, pm, tgt); + tgt[0] = (tgt[0] & mask) | SMALL_HEADER(len+1,reso) | (x << 8); + } + else { + // the length does not fit + if (bits <= 32) { + // the new element fits into the first word + uint64_t mask = nbit_compl_mask(32+bits); + shift_left_strict(bits+24, n, src, pm, tgt); + tgt[0] = (tgt[0] & mask) | BIG_HEADER(len+1,reso) | (x << 32); + } + else { + // the new element does not fit into the first word + uint64_t mask = nbit_compl_mask(bits-32); + shift_left_strict(bits+24, n, src, pm, tgt); + tgt[0] = BIG_HEADER(len+1,reso) | (x << 32); + tgt[1] = (tgt[1] & mask) | (x >> 32); + } + } + } + else { + // the old vector is big + if (bits <= 32) { + // the new element fits into the first word + shift_left_strict(bits, n, src, pm, tgt); + if (bits < 32) { + uint64_t mask = nbit_compl_mask(32+bits); + tgt[0] = (tgt[0] & mask) | BIG_HEADER(len+1,reso) | (x << 32); + } + else { + tgt[0] = BIG_HEADER(len+1,reso) | (x << 32); + } + } + else { + // the new element does not fit into the first word + uint64_t mask = nbit_compl_mask(bits-32); + shift_left_strict(bits, n, src, pm, tgt); + tgt[0] = BIG_HEADER(len+1,reso) | (x << 32); + tgt[1] = (tgt[1] & mask) | (x >> 32); + } + } + } + else { + // the new element needs more bits + if ( (x_bits <= MAX_SMALL_BITS) && (len+1 <= MAX_SMALL_LENGTH) ) { + // but we still fit into a small vector + tgt[0] = SMALL_HEADER(len+1,x_reso) | (x << 8); + copy_elements_into + ( len , bits , src , header_bits + , pm , x_bits , tgt , 8 + x_bits + ); + } + else { + // we need a big vector + tgt[0] = BIG_HEADER(len+1,x_reso) | (x << 32); + uint64_t y = (x >> 32); + if (y > 0) { tgt[1] = y; } + copy_elements_into + ( len , bits , src , header_bits + , pm , x_bits , tgt , 32 + x_bits + ); + } + } +} + +// ----------------------------------------------------------------------------- +// SNOC + +#define SNOC_WRITE(tgt_ofs,y_bits) \ + int bit_ofs = (tgt_ofs); \ + int word_ofs = (bit_ofs) >> 6; \ + bit_ofs &= 63; \ + int new_ofs = bit_ofs + y_bits; \ + if (new_ofs <= 64) { \ + uint64_t mask = nbit_mask(bit_ofs); \ + tgt[word_ofs] = (tgt[word_ofs] & mask) | (x << bit_ofs); \ + } \ + else { \ + uint64_t mask = nbit_mask(bit_ofs); \ + tgt[word_ofs ] = (tgt[word_ofs] & mask) | (x << bit_ofs); \ + tgt[word_ofs+1] = (x >> (64 - bit_ofs)); \ + } \ + if (new_ofs <= 64) { \ + *pm = word_ofs + 1; \ + } \ + else { \ + *pm = word_ofs + 2; \ + } + + +void vec_snoc(uint64_t x, int n, const uint64_t *src, int* pm, uint64_t *tgt) +{ + VEC_HEADER_CODE(src) + + int x_bits = required_bits(x); + int x_reso = required_reso(x); + + if (len==0) { + // snoc to an empty vector + if (x_reso <= 2) { + tgt[0] = SMALL_HEADER(1,x_reso) | (x << 8); + *pm = 1; + return; + } + else { + tgt[0] = BIG_HEADER(1,x_reso) | (x << 32); + uint64_t y = (x >> 32); + if (y==0) { + *pm = 1; + } + else { + tgt[1] = y; + *pm = 2; + } + return; + } + } + + if (x_bits <= bits) { + // the new element fits without changing the resolution + if (is_small) { + // the old vector is small + if (len+1 <= MAX_SMALL_LENGTH) { + // the length fits, too + memcpy(tgt, src, n<<3); + uint64_t mask = nbit_compl_mask(8); + tgt[0] = (tgt[0] & mask) | SMALL_HEADER(len+1,reso); + SNOC_WRITE( 8+bits*len , bits ) + } + else { + // the length does not fit + shift_left_strict(24, n, src, pm, tgt); + uint64_t mask = nbit_compl_mask(32); + tgt[0] = (tgt[0] & mask) | BIG_HEADER(len+1,reso); + SNOC_WRITE( 32+bits*len , bits ) + } + } + else { + // the old vector is big + memcpy(tgt, src, n<<3); + uint64_t mask = nbit_compl_mask(32); + tgt[0] = (tgt[0] & mask) | BIG_HEADER(len+1,reso); + SNOC_WRITE( 32+bits*len , bits ) + } + } + else { + // the new element needs more bits + if ( (x_bits <= MAX_SMALL_BITS) && (len+1 <= MAX_SMALL_LENGTH) ) { + // but we still fit into a small vector + tgt[0] = SMALL_HEADER(len+1,x_reso); + copy_elements_into + ( len , bits , src , header_bits + , pm , x_bits , tgt , 8 + ); + SNOC_WRITE(8 + x_bits*len, x_bits) + } + else { + // we need a big vector + tgt[0] = BIG_HEADER(len+1,x_reso); + copy_elements_into + ( len , bits , src , header_bits + , pm , x_bits , tgt , 32 + ); + SNOC_WRITE(32 + x_bits*len, x_bits) + } + } +} + +// ----------------------------------------------------------------------------- +// folds + +uint64_t vec_max(int n, const uint64_t *src) +{ + VEC_HEADER_CODE(src) + uint64_t max = 0; + VEC_READ_LOOP + max = (elem > max) ? elem : max; + } + return max; +} + +uint64_t vec_sum(int n, const uint64_t *src) +{ + VEC_HEADER_CODE(src) + uint64_t sum = 0; + VEC_READ_LOOP + sum += elem; + } + return sum; +} + +// ----------------------------------------------------------------------------- +// zipping folds + +// strictly equal (as vectors) +uint64_t vec_equal_strict(int n1, const uint64_t *src1, int n2, const uint64_t *src2) +{ + int len1,bits1,header_bits1; + int len2,bits2,header_bits2; + + { VEC_HEADER_CODE(src1) ; len1 = len ; bits1 = bits ; header_bits1 = header_bits; } + { VEC_HEADER_CODE(src2) ; len2 = len ; bits2 = bits ; header_bits2 = header_bits; } + + if (len1 != len2) { + return 0; + } + + int bool = 1; + int zip_len = len1; + VEC_ZIP_LOOP + if (elem1 != elem2) { + bool = 0; + break; + } + } + return bool; +} + +// equal when extended by zeros (as monomials, partitions, etc) +uint64_t vec_equal_extzero(int n1, const uint64_t *src1, int n2, const uint64_t *src2) +{ + int len1,bits1,header_bits1; + int len2,bits2,header_bits2; + + { VEC_HEADER_CODE(src1) ; len1 = len ; bits1 = bits ; header_bits1 = header_bits; } + { VEC_HEADER_CODE(src2) ; len2 = len ; bits2 = bits ; header_bits2 = header_bits; } + + int bool = 1; + int zip_len = (len1>=len2) ? len1 : len2; + VEC_ZIP_LOOP + if (elem1 != elem2) { + bool = 0; + break; + } + } + return bool; +} + +// strict comparison (as vectors): +// first compares the length, then if equal, lexicographically the sentences +// returns: +// -1 = LT +// 0 = EQ +// +1 = GT +uint64_t vec_compare_strict(int n1, const uint64_t *src1, int n2, const uint64_t *src2) +{ + int len1,bits1,header_bits1; + int len2,bits2,header_bits2; + + { VEC_HEADER_CODE(src1) ; len1 = len ; bits1 = bits ; header_bits1 = header_bits; } + { VEC_HEADER_CODE(src2) ; len2 = len ; bits2 = bits ; header_bits2 = header_bits; } + + if (len1 < len2) { return LT; } + if (len1 > len2) { return GT; } + + int result = EQ; + int zip_len = len1; + VEC_ZIP_LOOP + if (elem1 < elem2) { result = LT; break; } + if (elem1 > elem2) { result = GT; break; } + } + return result; +} + +// lexicographically compare sequences extended to infinity with zeros +uint64_t vec_compare_extzero(int n1, const uint64_t *src1, int n2, const uint64_t *src2) +{ + int len1,bits1,header_bits1; + int len2,bits2,header_bits2; + + { VEC_HEADER_CODE(src1) ; len1 = len ; bits1 = bits ; header_bits1 = header_bits; } + { VEC_HEADER_CODE(src2) ; len2 = len ; bits2 = bits ; header_bits2 = header_bits; } + + int result = EQ; + int zip_len = (len1>=len2) ? len1 : len2; + VEC_ZIP_LOOP + if (elem1 < elem2) { result = LT; break; } + if (elem1 > elem2) { result = GT; break; } + } + return result; +} + +// pointwise less or equal, extended to infinity with zeros +uint64_t vec_less_or_equal(int n1, const uint64_t *src1, int n2, const uint64_t *src2) +{ + int len1,bits1,header_bits1; + int len2,bits2,header_bits2; + + { VEC_HEADER_CODE(src1) ; len1 = len ; bits1 = bits ; header_bits1 = header_bits; } + { VEC_HEADER_CODE(src2) ; len2 = len ; bits2 = bits ; header_bits2 = header_bits; } + + int bool = 1; + int zip_len = (len1>=len2) ? len1 : len2; + VEC_ZIP_LOOP + if (elem1 > elem2) { + bool = 0; + break; + } + if (i >= len1-1) { break; } // if we are over the first list, then 0 <= anything + } + return bool; +} + +// dominance order of partitions +uint64_t vec_partial_sums_less_or_equal(int n1, const uint64_t *src1, int n2, const uint64_t *src2) +{ + int len1,bits1,header_bits1; + int len2,bits2,header_bits2; + + { VEC_HEADER_CODE(src1) ; len1 = len ; bits1 = bits ; header_bits1 = header_bits; } + { VEC_HEADER_CODE(src2) ; len2 = len ; bits2 = bits ; header_bits2 = header_bits; } + + int bool = 1; + int zip_len = (len1>=len2) ? len1 : len2; + uint64_t sum1 = 0; + uint64_t sum2 = 0; + VEC_ZIP_LOOP + sum1 += elem1; + sum2 += elem2; + if (sum1 > sum2) { + bool = 0; + break; + } + if (i >= len1-1) { break; } // if we are over the first list, then sum(0) <= sum(anything) + } + return bool; +} + +// ----------------------------------------------------------------------------- + +void vec_add(int n1, const uint64_t *src1, int n2, const uint64_t *src2, int *pm, uint64_t *tgt) +{ + int len1,bits1,header_bits1; + int len2,bits2,header_bits2; + + { VEC_HEADER_CODE(src1) ; len1 = len ; bits1 = bits ; header_bits1 = header_bits; } + { VEC_HEADER_CODE(src2) ; len2 = len ; bits2 = bits ; header_bits2 = header_bits; } + + int zip_len = (len1>=len2) ? len1 : len2; + + // compute upper bound for the result + uint64_t bound = 0; + { VEC_ZIP_LOOP + uint64_t x = elem1 + elem2; + if (x > bound) { bound = x; } + } + } + + int tgt_bits = required_bits(bound); + int tgt_reso = bits2reso(tgt_bits); + + uint64_t *q = tgt; + int q_ofs; + + // write header + if ( (tgt_bits <= MAX_SMALL_BITS) && (zip_len <= MAX_SMALL_LENGTH) ) { + q[0] = SMALL_HEADER( zip_len , tgt_reso ); + q_ofs = 8; + } + else { + q[0] = BIG_HEADER( zip_len , tgt_reso ); + q_ofs = 32; + } + + // write result + VEC_ZIP_LOOP + uint64_t y = elem1 + elem2; + WRITE_ELEMENT(y) + } + + STORE_OUTPUT_LENGTH(pm) +} + +// ----------------------------------------------------------------------------- + +// subtraction with overflow indicator +int vec_sub_overflow(int n1, const uint64_t *src1, int n2, const uint64_t *src2, int *pm, uint64_t *tgt) +{ + int len1,bits1,header_bits1; + int len2,bits2,header_bits2; + + { VEC_HEADER_CODE(src1) ; len1 = len ; bits1 = bits ; header_bits1 = header_bits; } + { VEC_HEADER_CODE(src2) ; len2 = len ; bits2 = bits ; header_bits2 = header_bits; } + + int zip_len = (len1>=len2) ? len1 : len2; + + // compute upper bound for the result + uint64_t bound = 0; + { VEC_ZIP_LOOP + uint64_t x = elem1 + elem2; + if (x > bound) { bound = x; } + } + } + + int tgt_bits = required_bits(bound); + int tgt_reso = bits2reso(tgt_bits); + + uint64_t *q = tgt; + int q_ofs; + + // write header + if ( (tgt_bits <= MAX_SMALL_BITS) && (zip_len <= MAX_SMALL_LENGTH) ) { + q[0] = SMALL_HEADER( zip_len , tgt_reso ); + q_ofs = 8; + } + else { + q[0] = BIG_HEADER( zip_len , tgt_reso ); + q_ofs = 32; + } + + int overflow = 0; + // write result + VEC_ZIP_LOOP + if (elem2 > elem1) { overflow = 1; } + uint64_t y = elem1 - elem2; + WRITE_ELEMENT(y) + } + + STORE_OUTPUT_LENGTH(pm) + return overflow; +} + +// ----------------------------------------------------------------------------- +// maps + +void vec_scale(uint64_t s, int n, const uint64_t *src, int *pm, uint64_t *tgt) +{ + VEC_HEADER_CODE(src) + + uint64_t sbnd = 1; + sbnd = sbnd << (64-bits); + int tgt_bits; + if (s <= sbnd) { + uint64_t bound = 1; + bound = bound << bits; + bound = s * (bound - 1); + tgt_bits = required_bits(bound); + } + else { + tgt_bits = 64; + } + int tgt_reso = bits2reso(tgt_bits); + + uint64_t *q = tgt; + int q_ofs; + // write header + if ( (tgt_bits <= MAX_SMALL_BITS) && (len <= MAX_SMALL_LENGTH) ) { + q[0] = SMALL_HEADER( len , tgt_reso ); + q_ofs = 8; + } + else { + q[0] = BIG_HEADER( len , tgt_reso ); + q_ofs = 32; + } + + VEC_READ_LOOP + uint64_t y = s * elem; + WRITE_ELEMENT(y) + } + + STORE_OUTPUT_LENGTH(pm) +} + +// ----------------------------------------------------------------------------- + +// x1, x1+x2, x1+x2+x3 +uint64_t vec_partial_sums(int n, const uint64_t *src, int *pm, uint64_t *tgt) +{ + VEC_HEADER_CODE(src) + + int overflow = 0; + uint64_t sum = 0, tmp; + { + VEC_READ_LOOP + tmp = sum + elem; + if (tmp < sum) { overflow = 1; } + sum = tmp; + } + } + + int tgt_bits; + if (overflow) { + tgt_bits = 64; + } + else { + tgt_bits = required_bits(sum); + } + int tgt_reso = bits2reso(tgt_bits); + + uint64_t *q = tgt; + int q_ofs; + // write header + if ( (tgt_bits <= MAX_SMALL_BITS) && (len <= MAX_SMALL_LENGTH) ) { + q[0] = SMALL_HEADER( len , tgt_reso ); + q_ofs = 8; + } + else { + q[0] = BIG_HEADER( len , tgt_reso ); + q_ofs = 32; + } + + uint64_t acc = 0; + VEC_READ_LOOP + acc += elem; + WRITE_ELEMENT(acc) + } + + STORE_OUTPUT_LENGTH(pm) + return sum; +} + +// ----------------------------------------------------------------------------- diff --git a/compact-word-vectors.cabal b/compact-word-vectors.cabal --- a/compact-word-vectors.cabal +++ b/compact-word-vectors.cabal @@ -1,10 +1,14 @@ Name: compact-word-vectors -Version: 0.1 +Version: 0.2 Synopsis: Small vectors of small integers stored very compactly. Description: A data structure to store small vectors of small integers - with minimal memory overhead. For example the vector + with minimal memory overhead. For example the (word) vector corresponding to [1..14] only takes 16 bytes (2 machine words on 64 bit architectures) of heap memory. + + See the module "Data.Vector.Compact.WordVec" for more + details. + License: BSD3 License-file: LICENSE Author: Balazs Komuves @@ -31,12 +35,15 @@ Data.Vector.Compact.IntVec Data.Vector.Compact.Blob - Default-Extensions: CPP, BangPatterns - Other-Extensions: MagicHash, UnboxedTuples - Default-Language: Haskell2010 + Default-Extensions: CPP, BangPatterns + Other-Extensions: MagicHash, ForeignFunctionInterface Hs-Source-Dirs: src + + C-Sources: cbits/blob.c + cbits/wordvec.c + Include-Dirs: cbits ghc-options: -fwarn-tabs -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-unused-imports diff --git a/src/Data/Vector/Compact/Blob.hs b/src/Data/Vector/Compact/Blob.hs --- a/src/Data/Vector/Compact/Blob.hs +++ b/src/Data/Vector/Compact/Blob.hs @@ -21,16 +21,65 @@ -- * tables indexed by such things -- -{-# LANGUAGE CPP, BangPatterns, MagicHash, UnboxedTuples #-} -module Data.Vector.Compact.Blob where +{-# LANGUAGE CPP, BangPatterns, MagicHash, ForeignFunctionInterface #-} +module Data.Vector.Compact.Blob + ( + -- * The Blob type + Blob(..) + , blobTag + , blobSizeInWords + , blobSizeInBytes + , blobSizeInBits + -- * Conversion to\/from lists + , blobFromWordList , blobFromWordListN + , blobToWordList + -- * Conversion to\/from 'ByteArray'-s + , blobFromByteArray + , blobToByteArray + -- * Equality comparison + , eqBlob + -- * Head, tail, cons, etc + , head + , tail + , last + , consWord + , snocWord + -- * Indexing + , indexWord , indexByte + , extractSmallWord , extractSmallWord64 + -- * Resizing + , extendToSize + , cutToSize + , forceToSize + -- * Higher-order functions + , mapBlob + , shortZipWith + , longZipWith + , unsafeZipWith + -- * Hexadecimal printing + , Hex(..) + , hexWord64 , hexWord64_ + -- * (Indirect) access to the raw data + -- + -- $raw + , peekBlob + , pokeBlob + -- * Wrappers for C implementations + -- + -- $wrapper + , CFun10 , CFun20 , CFun11 , CFun21 , CFun11_ , CFun21_ + , wrapCFun10 , wrapCFun20 , wrapCFun11 , wrapCFun21 , wrapCFun11_ , wrapCFun21_ + ) + where -------------------------------------------------------------------------------- +import Prelude hiding ( head , tail , last ) import Data.Char import Data.Bits import Data.Int import Data.Word -import Data.List +import qualified Data.List as L import Control.Monad import Control.Monad.ST @@ -41,9 +90,14 @@ import GHC.Exts import GHC.IO +import Foreign.C.Types import Foreign.Ptr +import Foreign.Storable +import Foreign.Marshal import Foreign.Marshal.Array +import System.IO.Unsafe as Unsafe + import Control.Monad.Primitive import Data.Primitive.ByteArray @@ -61,7 +115,7 @@ | Blob4 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 | Blob5 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 | Blob6 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 - | BlobN !(ByteArray#) + | BlobN {-# UNPACK #-} !ByteArray -------------------------------------------------------------------------------- @@ -70,13 +124,13 @@ -- | Number of 'Word64'-s blobSizeInWords :: Blob -> Int -blobSizeInWords blob = case blob of - BlobN arr -> shiftR (I# (sizeofByteArray# arr)) 3 +blobSizeInWords !blob = case blob of + BlobN !arr -> shiftR (sizeofByteArray arr) 3 otherwise -> blobTag blob + 1 blobSizeInBytes :: Blob -> Int -blobSizeInBytes blob = case blob of - BlobN arr -> I# (sizeofByteArray# arr) +blobSizeInBytes !blob = case blob of + BlobN !arr -> sizeofByteArray arr otherwise -> shiftL (blobTag blob + 1) 3 blobSizeInBits :: Blob -> Int @@ -89,7 +143,7 @@ blobFromWordList ws = blobFromWordListN (length ws) ws blobFromWordListN :: Int -> [Word64] -> Blob -blobFromWordListN n ws = case n of +blobFromWordListN !n ws = case n of 0 -> Blob1 0 1 -> case ws of { (a:_) -> Blob1 a } 2 -> case ws of { (a:b:_) -> Blob2 a b } @@ -97,8 +151,7 @@ 4 -> case ws of { (a:b:c:d:_) -> Blob4 a b c d } 5 -> case ws of { (a:b:c:d:e:_) -> Blob5 a b c d e } 6 -> case ws of { (a:b:c:d:e:f:_) -> Blob6 a b c d e f } - _ -> case byteArrayFromListN n ws of - ByteArray ba# -> BlobN ba# + _ -> BlobN (byteArrayFromListN n ws) blobToWordList :: Blob -> [Word64] blobToWordList blob = case blob of @@ -108,18 +161,17 @@ Blob4 a b c d -> a:b:c:d:[] Blob5 a b c d e -> a:b:c:d:e:[] Blob6 a b c d e f -> a:b:c:d:e:f:[] - BlobN ba# -> foldrByteArray (:) [] (ByteArray ba#) + BlobN ba -> foldrByteArray (:) [] ba -------------------------------------------------------------------------------- -- * Conversion to\/from @ByteArray@-s -- | Note: we pad the input with zero bytes, assuming little-endian architecture. blobFromByteArray :: ByteArray -> Blob -blobFromByteArray ba@(ByteArray ba#) +blobFromByteArray !ba | nwords > 6 = if nwords1 == nwords - then BlobN ba# - else let ByteArray new# = byteArrayFromListN nwords words - in BlobN new# + then BlobN ba + else BlobN (byteArrayFromListN nwords words ) | nwords == 0 = Blob1 0 | otherwise = blobFromWordListN nwords words where @@ -129,23 +181,27 @@ words :: [Word64] words = if nwords1 == nwords - then foldrByteArray (:) [] (ByteArray ba#) + then foldrByteArray (:) [] ba else let !ofs = shiftL nwords1 3 !m = nbytes - ofs w8_to_w64 :: Word8 -> Word64 w8_to_w64 = fromIntegral - !lastWord = foldl' (.|.) 0 [ shiftL (w8_to_w64 (indexByteArray ba (ofs + i))) (shiftL i 3) | i<-[0..m-1] ] - in foldrByteArray (:) [] (ByteArray ba#) ++ [lastWord] + !lastWord = L.foldl' (.|.) 0 + [ shiftL (w8_to_w64 (indexByteArray ba (ofs + i))) (shiftL i 3) + | i<-[0..m-1] + ] + in foldrByteArray (:) [lastWord] ba blobToByteArray :: Blob -> ByteArray -blobToByteArray blob = case blob of - BlobN ba# -> ByteArray ba# - _ -> byteArrayFromListN (blobSizeInWords blob) (blobToWordList blob) +blobToByteArray !blob = case blob of + BlobN ba -> ba + _ -> byteArrayFromListN (blobSizeInWords blob) (blobToWordList blob) -------------------------------------------------------------------------------- +-- * Instances instance Show Blob where - showsPrec prec blob + showsPrec prec !blob = showParen (prec > 10) $ showString "blobFromWordList " . shows (map Hex $ blobToWordList blob) @@ -163,11 +219,11 @@ ( Blob4 a b c d , Blob4 p q r s ) -> a==p && b==q && c==r && d==s ( Blob5 a b c d e , Blob5 p q r s t ) -> a==p && b==q && c==r && d==s && e==t ( Blob6 a b c d e f , Blob6 p q r s t u ) -> a==p && b==q && c==r && d==s && e==t && f==u - ( BlobN one# , BlobN two# ) -> ByteArray one# == ByteArray two# + ( BlobN one , BlobN two ) -> one == two _ -> error "FATAL ERROR: should not happen" -------------------------------------------------------------------------------- --- * Hexadecimal +-- * Hexadecimal printing newtype Hex = Hex Word64 @@ -190,30 +246,30 @@ -------------------------------------------------------------------------------- -- * Peek -peekWord :: Blob -> Int -> Word64 -peekWord blob idx = case blob of +indexWord :: Blob -> Int -> Word64 +indexWord !blob !idx = case blob of Blob1 a | idx == 0 -> a - | otherwise -> error "Blob/peekWord: index out of bounds" + | otherwise -> error "Blob/indexWord: index out of bounds" Blob2 a b | idx == 0 -> a | idx == 1 -> b - | otherwise -> error "Blob/peekWord: index out of bounds" + | otherwise -> error "Blob/indexWord: index out of bounds" Blob3 a b c | idx == 0 -> a | idx == 1 -> b | idx == 2 -> c - | otherwise -> error "Blob/peekWord: index out of bounds" + | otherwise -> error "Blob/indexWord: index out of bounds" Blob4 a b c d | idx == 0 -> a | idx == 1 -> b | idx == 2 -> c | idx == 3 -> d - | otherwise -> error "Blob/peekWord: index out of bounds" + | otherwise -> error "Blob/indexWord: index out of bounds" Blob5 a b c d e | idx == 0 -> a @@ -221,7 +277,7 @@ | idx == 2 -> c | idx == 3 -> d | idx == 4 -> e - | otherwise -> error "Blob/peekWord: index out of bounds" + | otherwise -> error "Blob/indexWord: index out of bounds" Blob6 a b c d e f | idx == 0 -> a @@ -230,31 +286,78 @@ | idx == 3 -> d | idx == 4 -> e | idx == 5 -> f - | otherwise -> error "Blob/peekWord: index out of bounds" + | otherwise -> error "Blob/indexWord: index out of bounds" - BlobN arr# -> indexByteArray (ByteArray arr#) idx + BlobN arr -> indexByteArray arr idx -- | NOTE: We assume a little-endian architecture here. -- Though it seems that since GHC does not gives us direct access to the closure, -- it doesn\'t matter after all... -- -peekByte :: Blob -> Int -> Word8 -peekByte blob idx = - let w = peekWord blob (shiftR idx 3) +indexByte :: Blob -> Int -> Word8 +indexByte !blob !idx = + let !w = indexWord blob (shiftR idx 3) in fromIntegral $ shiftR w (8 * (idx .&. 7)) -blobHead :: Blob -> Word64 -blobHead blob = case blob of +-------------------------------------------------------------------------------- +-- * Head and last + +head :: Blob -> Word64 +head blob = case blob of Blob1 a -> a Blob2 a _ -> a Blob3 a _ _ -> a Blob4 a _ _ _ -> a Blob5 a _ _ _ _ -> a Blob6 a _ _ _ _ _ -> a - BlobN arr# -> indexByteArray (ByteArray arr#) 0 - + BlobN arr -> indexByteArray arr 0 + +last :: Blob -> Word64 +last blob = case blob of + Blob1 z -> z + Blob2 _ z -> z + Blob3 _ _ z -> z + Blob4 _ _ _ z -> z + Blob5 _ _ _ _ z -> z + Blob6 _ _ _ _ _ z -> z + BlobN arr -> indexByteArray arr (blobSizeInWords blob - 1) + -------------------------------------------------------------------------------- +-- * Cons, Snoc, tail +-- | Prepend a word at the start +consWord :: Word64 -> Blob -> Blob +consWord !y !blob = case blob of + Blob1 a -> Blob2 y a + Blob2 a b -> Blob3 y a b + Blob3 a b c -> Blob4 y a b c + Blob4 a b c d -> Blob5 y a b c d + Blob5 a b c d e -> Blob6 y a b c d e + _ -> wrapCFun11_ (c_cons y) (+1) blob + +-- | Append a word at the end +snocWord :: Blob -> Word64 -> Blob +snocWord !blob !z = case blob of + Blob1 a -> Blob2 a z + Blob2 a b -> Blob3 a b z + Blob3 a b c -> Blob4 a b c z + Blob4 a b c d -> Blob5 a b c d z + Blob5 a b c d e -> Blob6 a b c d e z + _ -> wrapCFun11_ (c_snoc z) (+1) blob + +-- | Remove the first word +tail :: Blob -> Blob +tail !blob = case blob of + Blob1 _ -> Blob1 0 + Blob2 _ b -> Blob1 b + Blob3 _ b c -> Blob2 b c + Blob4 _ b c d -> Blob3 b c d + Blob5 _ b c d e -> Blob4 b c d e + Blob6 _ b c d e f -> Blob5 b c d e f + _ -> wrapCFun11_ c_tail id blob + +-------------------------------------------------------------------------------- + -- | @extractSmallWord n blob ofs@ extracts a small word of @n@ bits starting from the -- @ofs@-th bit. This should satisfy -- @@ -264,27 +367,188 @@ -- than the size (in bits) of the blob. -- extractSmallWord :: Integral a => Int -> Blob -> Int -> a -extractSmallWord n blob ofs = fromIntegral (extractSmallWord64 n blob ofs) +extractSmallWord !n !blob !ofs = fromIntegral (extractSmallWord64 n blob ofs) extractSmallWord64 :: Int -> Blob -> Int -> Word64 extractSmallWord64 !n !blob !ofs - | q2 == q1 = mask .&. shiftR (peekWord blob q1) r1 - | q2 == q1 + 1 = mask .&. (shiftR (peekWord blob q1) r1 .|. shiftL (peekWord blob q2) (64-r1)) + | q2 == q1 = mask .&. shiftR (indexWord blob q1) r1 + | q2 == q1 + 1 = mask .&. (shiftR (indexWord blob q1) r1 .|. shiftL (indexWord blob q2) (64-r1)) | otherwise = error "Blob/extractSmallWord: FATAL ERROR" where - mask = shiftL 1 n - 1 - end = ofs + n - 1 - q1 = shiftR ofs 6 - q2 = shiftR end 6 - r1 = ofs .&. 63 + !mask = shiftL 1 n - 1 + !end = ofs + n - 1 + !q1 = shiftR ofs 6 + !q2 = shiftR end 6 + !r1 = ofs .&. 63 +{- -- | An alternate implementation using 'testBit', for testing purposes only extractSmallWord64_naive :: Int -> Blob -> Int -> Word64 extractSmallWord64_naive n blob ofs = sum [ shiftL 1 i | i<-[0..n-1] , testBit blob (ofs+i) ] +-} -------------------------------------------------------------------------------- --- * change size +-- * (Indirect) access to the raw data +-- +-- $raw +-- +-- Note: Because GHC does not support direct manipulation of heap data +-- (the garbage collector can move it anytime), these involve copying. +-- +pokeBlob :: Ptr Word64 -> Blob -> IO Int +pokeBlob !ptr !blob = case blob of + Blob1 a -> poke ptr a >> return 1 + Blob2 a b -> pokeArray ptr [a,b] >> return 2 + Blob3 a b c -> pokeArray ptr [a,b,c] >> return 3 + Blob4 a b c d -> pokeArray ptr [a,b,c,d] >> return 4 + Blob5 a b c d e -> pokeArray ptr [a,b,c,d,e] >> return 5 + Blob6 a b c d e f -> pokeArray ptr [a,b,c,d,e,f] >> return 6 + BlobN ba -> let !nbytes = sizeofByteArray ba + in copyByteArrayToPtr ba 0 ptr nbytes >> return (shiftR nbytes 3) +peekBlob :: Int -> Ptr Word64 -> IO Blob +peekBlob !n !ptr = + case n of + 0 -> return (Blob1 0) + 1 -> peek ptr >>= \a -> return (Blob1 a) + 2 -> peekArray 2 ptr >>= \[a,b] -> return (Blob2 a b) + 3 -> peekArray 3 ptr >>= \[a,b,c] -> return (Blob3 a b c) + 4 -> peekArray 4 ptr >>= \[a,b,c,d] -> return (Blob4 a b c d) + 5 -> peekArray 5 ptr >>= \[a,b,c,d,e] -> return (Blob5 a b c d e) + 6 -> peekArray 6 ptr >>= \[a,b,c,d,e,f] -> return (Blob6 a b c d e f) + _ -> do + mut <- newByteArray (shiftL n 3) + copyPtrToByteArray ptr mut 0 (shiftL n 3) + ba <- unsafeFreezeByteArray mut + return (BlobN ba) + +-------------------------------------------------------------------------------- +-- * Wrappers for C implementations +-- +-- $wrapper +-- +-- As above, these involve copying of the data (both inputs and outputs); +-- so they first allocate temporary buffers, copy the data into them +-- call the C function, and copy the result to a new 'Blob'. +-- +-- Naming conventions: For example @CFun21@ means 2 Blob inputs and 1 Blob output. +-- +type CFun10 a = CInt -> Ptr Word64 -> IO a +type CFun20 a = CInt -> Ptr Word64 -> CInt -> Ptr Word64 -> IO a +type CFun11 a = CInt -> Ptr Word64 -> Ptr CInt -> Ptr Word64 -> IO a +type CFun21 a = CInt -> Ptr Word64 -> + CInt -> Ptr Word64 -> Ptr CInt -> Ptr Word64 -> IO a + +type CFun11_ = CFun11 () +type CFun21_ = CFun21 () + +-- | Allocate a temporary buffer, copy the content of the Blob there, +-- and call the C function +wrapCFun10_IO :: CFun10 a -> Blob -> IO a +wrapCFun10_IO action blob = do + let !n = blobSizeInWords blob + allocaArray n $ \ptr1 -> do + pokeBlob ptr1 blob + action (fromIntegral n) ptr1 + +-- | Allocate two temporary buffers, copy the content of the two Blobs there, +-- and call the C function +wrapCFun20_IO :: CFun20 a -> Blob -> Blob -> IO a +wrapCFun20_IO action blob1 blob2 = do + let !n1 = blobSizeInWords blob1 + let !n2 = blobSizeInWords blob2 + allocaArray n1 $ \ptr1 -> do + pokeBlob ptr1 blob1 + allocaArray n2 $ \ptr2 -> do + pokeBlob ptr2 blob2 + action (fromIntegral n1) ptr1 (fromIntegral n2) ptr2 + +-- | Allocate a temporary buffer, copy the content of the Blob there (unfortunately +-- we have to do this, because the GHC runtime does not allow direct manipulation of the heap, +-- even though we /know/ the heap layout...); then allocate another temporary buffer of +-- the given length (measured in words), call the C function which can fill this second +-- buffer, finally create a new Blob from the content of the second buffer +-- (another copying happens here). +-- +wrapCFun11_IO :: CFun11 a -> Int -> Blob -> IO (a,Blob) +wrapCFun11_IO action m blob = do + let !n = blobSizeInWords blob + allocaArray n $ \ptr1 -> do + pokeBlob ptr1 blob + allocaArray m $ \ptr2 -> do + alloca $ \q -> do + y <- action (fromIntegral n) ptr1 q ptr2 + k <- peek q + new <- peekBlob (fromIntegral k) ptr2 + return (y,new) + +wrapCFun21_IO :: CFun21 a -> Int -> Blob -> Blob -> IO (a,Blob) +wrapCFun21_IO action m blob1 blob2 = do + let !n1 = blobSizeInWords blob1 + allocaArray n1 $ \ptr1 -> do + pokeBlob ptr1 blob1 + let !n2 = blobSizeInWords blob2 + allocaArray n2 $ \ptr2 -> do + pokeBlob ptr2 blob2 + allocaArray m $ \ptr3 -> do + alloca $ \q -> do + y <- action (fromIntegral n1) ptr1 (fromIntegral n2) ptr2 q ptr3 + k <- peek q + new <- peekBlob (fromIntegral k) ptr3 + return (y,new) + +{-# NOINLINE wrapCFun10 #-} +wrapCFun10 :: CFun10 a -> Blob -> a +wrapCFun10 action blob = Unsafe.unsafePerformIO $ wrapCFun10_IO action blob + +{-# NOINLINE wrapCFun20 #-} +wrapCFun20 :: CFun20 a -> Blob -> Blob -> a +wrapCFun20 action blob1 blob2 = Unsafe.unsafePerformIO $ wrapCFun20_IO action blob1 blob2 + +{-# NOINLINE wrapCFun11 #-} +wrapCFun11 :: CFun11 a -> (Int -> Int) -> Blob -> (a,Blob) +wrapCFun11 action f blob = Unsafe.unsafePerformIO $ do + let !n = blobSizeInWords blob + wrapCFun11_IO action (f n) blob + +{-# NOINLINE wrapCFun11_ #-} +wrapCFun11_ :: CFun11_ -> (Int -> Int) -> Blob -> Blob +wrapCFun11_ action f blob = Unsafe.unsafePerformIO $ do + let !n = blobSizeInWords blob + snd <$> wrapCFun11_IO action (f n) blob + +{-# NOINLINE wrapCFun21 #-} +wrapCFun21 :: CFun21 a -> (Int -> Int -> Int) -> Blob -> Blob -> (a,Blob) +wrapCFun21 action f blob1 blob2 = Unsafe.unsafePerformIO $ do + let !n1 = blobSizeInWords blob1 + let !n2 = blobSizeInWords blob2 + wrapCFun21_IO action (f n1 n2) blob1 blob2 + +{-# NOINLINE wrapCFun21_ #-} +wrapCFun21_ :: CFun21_ -> (Int -> Int -> Int) -> Blob -> Blob -> Blob +wrapCFun21_ action f blob1 blob2 = Unsafe.unsafePerformIO $ do + let !n1 = blobSizeInWords blob1 + let !n2 = blobSizeInWords blob2 + snd <$> wrapCFun21_IO action (f n1 n2) blob1 blob2 + +-------------------------------------------------------------------------------- + +foreign import ccall unsafe "identity" c_identity :: CFun11_ -- for testing + +foreign import ccall unsafe "tail" c_tail :: CFun11_ +foreign import ccall unsafe "cons" c_cons :: Word64 -> CFun11_ +foreign import ccall unsafe "snoc" c_snoc :: Word64 -> CFun11_ + +foreign import ccall unsafe "rotate_left" c_rotate_left :: CInt -> CFun11_ +foreign import ccall unsafe "rotate_right" c_rotate_right :: CInt -> CFun11_ + +foreign import ccall unsafe "shift_left_strict" c_shift_left_strict :: CInt -> CFun11_ +foreign import ccall unsafe "shift_left_nonstrict" c_shift_left_nonstrict :: CInt -> CFun11_ +foreign import ccall unsafe "shift_right" c_shift_right :: CInt -> CFun11_ + +-------------------------------------------------------------------------------- +-- * Resizing + extendToSize :: Int -> Blob -> Blob extendToSize tgt blob | n >= tgt = blob @@ -311,20 +575,19 @@ -- * map and zipWith mapBlob :: (Word64 -> Word64) -> Blob -> Blob -mapBlob f blob = case blob of +mapBlob f !blob = case blob of Blob1 a -> Blob1 (f a) Blob2 a b -> Blob2 (f a) (f b) Blob3 a b c -> Blob3 (f a) (f b) (f c) Blob4 a b c d -> Blob4 (f a) (f b) (f c) (f d) Blob5 a b c d e -> Blob5 (f a) (f b) (f c) (f d) (f e) Blob6 a b c d e y -> Blob6 (f a) (f b) (f c) (f d) (f e) (f y) - BlobN arr# -> runST $ do + BlobN ba -> runST $ do let !n = blobSizeInWords blob - let ba = ByteArray arr# mut <- newByteArray (shiftL n 3) forM_ [0..n-1] $ \i -> writeByteArray mut i $ f (indexByteArray ba i) - new@(ByteArray new#) <- unsafeFreezeByteArray mut - return (BlobN new#) + new <- unsafeFreezeByteArray mut + return (BlobN new) shortZipWith :: (Word64 -> Word64 -> Word64) -> Blob -> Blob -> Blob shortZipWith f !blob1 !blob2 @@ -335,6 +598,7 @@ n1 = blobSizeInWords blob1 n2 = blobSizeInWords blob2 +-- | Extend the shorter blob with zeros longZipWith :: (Word64 -> Word64 -> Word64) -> Blob -> Blob -> Blob longZipWith f !blob1 !blob2 | n1 == n2 = unsafeZipWith f blob1 blob2 @@ -344,6 +608,7 @@ n1 = blobSizeInWords blob1 n2 = blobSizeInWords blob2 +-- | We assume that the two blobs has the same size! unsafeZipWith :: (Word64 -> Word64 -> Word64) -> Blob -> Blob -> Blob unsafeZipWith f !blob1 !blob2 = case (blob1,blob2) of ( Blob1 a , Blob1 p ) -> Blob1 (f a p) @@ -352,29 +617,30 @@ ( Blob4 a b c d , Blob4 p q r s ) -> Blob4 (f a p) (f b q) (f c r) (f d s) ( Blob5 a b c d e , Blob5 p q r s t ) -> Blob5 (f a p) (f b q) (f c r) (f d s) (f e t) ( Blob6 a b c d e y , Blob6 p q r s t u ) -> Blob6 (f a p) (f b q) (f c r) (f d s) (f e t) (f y u) - ( BlobN one# , BlobN two# ) -> + ( BlobN ba1 , BlobN ba2 ) -> runST $ do let !n = blobSizeInWords blob1 - ba1 = ByteArray one# - ba2 = ByteArray two# mut <- newByteArray (shiftL n 3) forM_ [0..n-1] $ \i -> writeByteArray mut i $ f (indexByteArray ba1 i) (indexByteArray ba2 i) - new@(ByteArray new#) <- unsafeFreezeByteArray mut - return (BlobN new#) + new <- unsafeFreezeByteArray mut + return (BlobN new) _ -> error "FATAL ERROR: should not happen" -------------------------------------------------------------------------------- +-- | Implementation note: When necessary, the bitwise operations consider the blobs +-- extended to infinity with zero withs. This is especially important with 'shiftL', +-- which may /NOT/ extend the blob size if the new bits are all zero. instance Bits Blob where (.&.) = shortZipWith (.&.) (.|.) = longZipWith (.|.) xor = longZipWith xor complement = mapBlob complement - shiftL = error "shiftR" - shiftR = error "shiftR" - rotateL = error "rotateL" - rotateR = error "rotateR" + shiftL blob k = wrapCFun11_ (c_shift_left_nonstrict (fromIntegral k)) f blob where f n = n + shiftR (k+63) 6 + shiftR blob k = wrapCFun11_ (c_shift_right (fromIntegral k)) id blob + rotateL blob k = wrapCFun11_ (c_rotate_left (fromIntegral k)) id blob + rotateR blob k = wrapCFun11_ (c_rotate_right (fromIntegral k)) id blob #if MIN_VERSION_base(4,12,0) bitSizeMaybe = Just . blobSizeInBits @@ -385,9 +651,9 @@ zeroBits = Blob1 0 isSigned _ = False - popCount blob = foldl' (+) 0 (map popCount $ blobToWordList blob) + popCount blob = L.foldl' (+) 0 (map popCount $ blobToWordList blob) - testBit !blob !k = if q >= n then False else testBit (peekWord blob q) r where + testBit !blob !k = if q >= n then False else testBit (indexWord blob q) r where (q,r) = divMod k 64 n = blobSizeInWords blob @@ -399,5 +665,23 @@ finiteBitSize = blobSizeInBits #endif + -------------------------------------------------------------------------------- - +-- * ByteArray helpers + +baToList :: ByteArray -> [Word64] +baToList = foldrByteArray (:) [] + +baSizeInWords :: ByteArray -> Int +baSizeInWords ba = shiftR (sizeofByteArray ba) 3 + +-- copyByteArrayToAddr# :: ByteArray# -> Int# -> Addr# -> Int# -> State# s -> State# s +-- copyAddrToByteArray# :: Addr# -> MutableByteArray# s -> Int# -> Int# -> State# s -> State# s + +copyByteArrayToPtr :: ByteArray -> Int -> Ptr a -> Int -> IO () +copyByteArrayToPtr (ByteArray ba#) (I# ofs) (Ptr p) (I# n) = primitive_ $ copyByteArrayToAddr# ba# ofs p n + +copyPtrToByteArray :: Ptr a -> MutableByteArray (PrimState IO) -> Int -> Int -> IO () +copyPtrToByteArray (Ptr p) (MutableByteArray mut#) (I# ofs) (I# n) = primitive_ $ copyAddrToByteArray# p mut# ofs n + +-------------------------------------------------------------------------------- diff --git a/src/Data/Vector/Compact/IntVec.hs b/src/Data/Vector/Compact/IntVec.hs --- a/src/Data/Vector/Compact/IntVec.hs +++ b/src/Data/Vector/Compact/IntVec.hs @@ -2,14 +2,53 @@ -- | Signed integer version of dynamic word vectors. -- -- See "Data.Vector.Compact.WordVec" for more details. +-- +-- Note: for unsigned integers, you really should use 'WordVec' instead, +-- because that is significantly faster, and has much more specialized functions +-- implemented. +-- +-- This module should be imported qualified (to avoid name clashes with Prelude). +-- {-# LANGUAGE BangPatterns #-} -module Data.Vector.Compact.IntVec where +module Data.Vector.Compact.IntVec + ( -- * The dynamic Word vector type + IntVec(..) + , Shape(..) + , vecShape + , vecLen , vecBits + -- * Show instance + , showIntVec , showsPrecIntVec + -- * Empty vector, singleton + , null , empty + , singleton , isSingleton + -- * Conversion to\/from lists + , fromList , fromList' , fromList'' + , lenMinMax + , toList , toRevList + -- * Indexing + , unsafeIndex , safeIndex + -- * Head, tail, etc + , head , tail , cons , uncons + , last , snoc -- init, unsnoc + , concat + -- * Generic operations + , fold + , naiveMap , boundedMap + , naiveZipWith , boundedZipWith , listZipWith + -- * Number of bits needed + , bitsNeededForMinMax + , bitsNeededFor + , roundBits + ) + where -------------------------------------------------------------------------------- +import Prelude hiding ( head , tail , init , last , null , concat ) +import qualified Data.List as L + import Data.Bits -import Data.List import Data.Vector.Compact.WordVec ( Shape(..) ) import qualified Data.Vector.Compact.WordVec as Dyn @@ -20,14 +59,25 @@ -- | A dynamic int vector is a (small) vector of small signed integers stored compactly newtype IntVec = IntVec Dyn.WordVec - deriving (Eq,Ord) - + -- deriving (Eq,Ord) -- WARNING: deriving Eq and Ord result in INCORRECT instances!! + +instance Eq IntVec where + (==) x y = (vecLen x == vecLen y) && (toList x == toList y) + +instance Ord IntVec where + compare x y = case compare (vecLen x) (vecLen y) of + LT -> LT + GT -> GT + EQ -> compare (toList x) (toList y) + vecShape :: IntVec -> Shape vecShape (IntVec dyn) = Dyn.vecShape dyn +-- | The length of the vector vecLen :: IntVec -> Int vecLen (IntVec dyn) = Dyn.vecLen dyn +-- | The number of bits per element used to encode the vector vecBits :: IntVec -> Int vecBits (IntVec dyn) = Dyn.vecBits dyn @@ -48,7 +98,7 @@ . shows (toList intvec) -------------------------------------------------------------------------------- --- * Empty +-- * Empty, singleton empty :: IntVec empty = fromList [] @@ -56,6 +106,14 @@ null :: IntVec -> Bool null v = vecLen v == 0 +singleton :: Int -> IntVec +singleton i = fromList [i] + +isSingleton :: IntVec -> Maybe Int +isSingleton (IntVec dynvec) = case Dyn.isSingleton dynvec of + Nothing -> Nothing + Just w -> Just $ word2int (Dyn.vecBits dynvec) w + -------------------------------------------------------------------------------- -- * Conversion from\/to lists @@ -63,20 +121,36 @@ toList (IntVec dynvec) = map (word2int bits) $ Dyn.toList dynvec where !bits = Dyn.vecBits dynvec +-- | @toRevList == reverse . toList@ +toRevList :: IntVec -> [Int] +toRevList (IntVec dynvec) = map (word2int bits) $ Dyn.toRevList dynvec where + !bits = Dyn.vecBits dynvec + +-- | Note: @fromList xs = fromList' (lenMinMax xs)@ fromList :: [Int] -> IntVec fromList xs = IntVec $ Dyn.fromList' (Dyn.Shape len bits) $ map (int2word bits) xs where (!len,!minMax) = lenMinMax xs !bits = roundBits (bitsNeededForMinMax minMax) +-- | usage: @fromList' (len,(min,max)) xs@ where @min@ and @max@ are the minimum and +-- maximum (or just a lower and upper bound) appearing in the list. fromList' :: (Int,(Int,Int)) -> [Int] -> IntVec fromList' (!len,!minMax) xs = IntVec $ Dyn.fromList' (Dyn.Shape len bits) $ map (int2word bits) xs where !bits = roundBits (bitsNeededForMinMax minMax) +-- | Don't use this unless you really know what you are doing! +fromList'' :: Shape -> [Int] -> IntVec +fromList'' shape@(Shape len !bits) xs = IntVec $ Dyn.fromList' shape $ map (int2word bits) xs + +-- | Computes the length, minimum and maximum of a list, traversing it only +-- once (instead of 3 times). lenMinMax :: [Int] -> (Int,(Int,Int)) lenMinMax = go 0 0 0 where go !cnt !p !q (x:xs) = go (cnt+1) (min x p) (max x q) xs go !cnt !p !q [] = (cnt,(p,q)) +-------------------------------------------------------------------------------- + int2word :: Int -> (Int -> Word) int2word !bits = i2w where !mask = shiftL 1 bits - 1 :: Word @@ -98,6 +172,7 @@ -------------------------------------------------------------------------------- -- * Indexing +-- | Indexing starts from 0. No bound checks are done. unsafeIndex :: Int -> IntVec -> Int unsafeIndex idx (IntVec dynvec) = word2int bits (Dyn.unsafeIndex idx dynvec) where !bits = Dyn.vecBits dynvec @@ -105,28 +180,98 @@ safeIndex :: Int -> IntVec -> Maybe Int safeIndex idx (IntVec dynvec) = (word2int bits) <$> (Dyn.safeIndex idx dynvec) where !bits = Dyn.vecBits dynvec + +-------------------------------------------------------------------------------- +-- * Head, tail, etc head :: IntVec -> Int head (IntVec dynvec) = word2int bits (Dyn.head dynvec) where !bits = Dyn.vecBits dynvec --------------------------------------------------------------------------------- +last :: IntVec -> Int +last (IntVec dynvec) = word2int bits (Dyn.last dynvec) where + !bits = Dyn.vecBits dynvec +tail :: IntVec -> IntVec +tail (IntVec dynvec) = IntVec (Dyn.tail dynvec) + +uncons :: IntVec -> Maybe (Int,IntVec) +uncons (IntVec dynvec) = case Dyn.uncons dynvec of + Nothing -> Nothing + Just (w,tl) -> Just (word2int bits w , IntVec tl) + where + bits = Dyn.vecBits dynvec + {- +-- | For testing purposes only +uncons_naive :: IntVec -> Maybe (Int,IntVec) +uncons_naive vec = if null vec + then Nothing + else Just (head vec, tail vec) +-} + +-- | Prepends an element +cons :: Int -> IntVec -> IntVec +cons k ivec@(IntVec vec) = IntVec $ Dyn.fromList' shape' $ map (int2word bits') (k : toList ivec) where + (Shape len bits) = Dyn.vecShape vec + bits' = roundBits $ max bits (bitsNeededFor k) + shape' = Shape (len+1) bits' + +-- | Appends an element +snoc :: IntVec -> Int -> IntVec +snoc ivec@(IntVec vec) k = IntVec $ Dyn.fromList' shape' $ map (int2word bits') (toList ivec ++ [k]) where + (Shape len bits) = Dyn.vecShape vec + bits' = roundBits $ max bits (bitsNeededFor k) + shape' = Shape (len+1) bits' + concat :: IntVec -> IntVec -> IntVec -concat u v = fromList' (Shape (lu+lv) (max bu bv)) (toList u ++ toList v) where +concat u v = fromList'' (Shape (lu+lv) (max bu bv)) (toList u ++ toList v) where Shape lu bu = vecShape u Shape lv bv = vecShape v --} + +-------------------------------------------------------------------------------- +-- * Generic operations + +-- | Left fold +fold :: (a -> Int -> a) -> a -> IntVec -> a +fold f x v = L.foldl' f x (toList v) + +naiveMap :: (Int -> Int) -> IntVec -> IntVec +naiveMap f u = fromList (map f $ toList u) + +-- | If you have (nearly sharp) lower and upper bounds for the result of your of function +-- on your vector, mapping can be more efficient +boundedMap :: (Int,Int) -> (Int -> Int) -> IntVec -> IntVec +boundedMap minMax f vec = fromList'' (Shape l bits) (toList vec) where + l = vecLen vec + bits = roundBits $ bitsNeededForMinMax minMax + +naiveZipWith :: (Int -> Int -> Int) -> IntVec -> IntVec -> IntVec +naiveZipWith f u v = fromList $ L.zipWith f (toList u) (toList v) + +-- | If you have (nearly sharp) lower and upper bounds for the result of your of function +-- on your vector, zipping can be more efficient +boundedZipWith :: (Int,Int) -> (Int -> Int -> Int) -> IntVec -> IntVec -> IntVec +boundedZipWith minMax f vec1 vec2 = fromList'' (Shape l bits) $ L.zipWith f (toList vec1) (toList vec2) where + l = min (vecLen vec1) (vecLen vec2) + bits = roundBits $ bitsNeededForMinMax minMax + +listZipWith :: (Int -> Int -> a) -> IntVec -> IntVec -> [a] +listZipWith f u v = L.zipWith f (toList u) (toList v) -------------------------------------------------------------------------------- -- * helpers for counting the necessary number of bits +-- | usage: @bitsNeededForMinMax (min,max)@ bitsNeededForMinMax :: (Int,Int) -> Int bitsNeededForMinMax (p,q) = max (bitsNeededFor p) (bitsNeededFor q) +-- | Note: this automatically rounds up to multiples of 4 bitsNeededFor :: Int -> Int -bitsNeededFor bound +bitsNeededFor = roundBits . bitsNeededFor' + +bitsNeededFor' :: Int -> Int +bitsNeededFor' bound | bound >= 0 = ceilingLog2 ( bound + 1) + 1 -- +8 needs 5 bits (-16..+15) | bound < 0 = ceilingLog2 (abs bound ) + 1 -- -8 needs 4 bits (-8 ..+7 ) where diff --git a/src/Data/Vector/Compact/WordVec.hs b/src/Data/Vector/Compact/WordVec.hs --- a/src/Data/Vector/Compact/WordVec.hs +++ b/src/Data/Vector/Compact/WordVec.hs @@ -5,7 +5,7 @@ -- This is data structure engineered to store large amount of -- small vectors of small elements compactly on memory. -- --- For example the list @[1..14] :: [Int]@ consumes 576 bytes (72 words) on +-- For example the list @[1..14] :: [Int]@ consumes 560 bytes (14x5=70 words) on -- a 64 bit machine, while the corresponding 'WordVec' takes only -- 16 bytes (2 words), and the one corresponding to @[101..115]@ still only -- 24 bytes (3 words). @@ -21,21 +21,89 @@ -- 20x improvement in memory usage). In any case the primary goal -- here is optimized memory usage. -- --- TODO: ability to add user-defined (fixed-length) header, it can be useful --- for some applications +-- This module should be imported qualified (to avoid name clashes with Prelude). -- +-- TODO: ability to add user-defined (fixed-length) header, it can be +-- potentially useful for some applications +-- -{-# LANGUAGE BangPatterns #-} -module Data.Vector.Compact.WordVec where +{-# LANGUAGE CPP, BangPatterns, ForeignFunctionInterface #-} +module Data.Vector.Compact.WordVec + ( -- * The dynamic Word vector type + WordVec(..) + , Shape(..) + , vecShape , vecShape' + , vecLen , vecBits , vecIsSmall + -- * Show instance + , showWordVec , showsPrecWordVec + -- * Empty vector, singleton + , null , empty + , singleton , isSingleton + -- * Conversion to\/from lists + , fromList , fromListN , fromList' + , toList , toRevList + -- * Indexing + , unsafeIndex , safeIndex + -- * Head, tail, etc + , head , tail , cons , uncons + , last , snoc -- init, unsnoc + , concat + -- * Specialized operations + -- + -- $spec + -- + -- ** Specialized folds + , sum , maximum + -- ** Specialized \"zipping folds\" + , eqStrict , eqExtZero + , cmpStrict , cmpExtZero + , lessOrEqual , partialSumsLessOrEqual + -- ** Specialized zips + , add , subtract + -- ** Specialized maps + , scale + -- ** Specialized scans + , partialSums + -- * Generic operations + , fold + , naiveMap , boundedMap + , naiveZipWith , boundedZipWith , listZipWith + -- * Number of bits needed + , bitsNeededFor , bitsNeededFor' + , roundBits + ) + where -------------------------------------------------------------------------------- +import Prelude hiding ( head , tail , init , last , null , concat , subtract , sum , maximum ) +import qualified Data.List as L + import Data.Bits import Data.Word -import Data.Vector.Compact.Blob +import Foreign.C +import Data.Vector.Compact.Blob hiding ( head , tail , last ) +import qualified Data.Vector.Compact.Blob as Blob + -------------------------------------------------------------------------------- + +-- ???? how to determine this properly... +-- why on earth isn't this stuff properly documented?!?!?!?!? +#ifdef x86_64_HOST_ARCH +#define MACHINE_WORD_BITS 64 +#elif i386_HOST_ARCH +#define MACHINE_WORD_BITS 32 +#elif i686_HOST_ARCH +#define MACHINE_WORD_BITS 32 +#elif aarch64_HOST_ARCH +#define MACHINE_WORD_BITS 64 +#else +#define MACHINE_WORD_BITS 32 +#endif + +-------------------------------------------------------------------------------- -- * The dynamic Word vector type -- | Dynamic word vectors are internally 'Blob'-s, which the first few bits @@ -50,8 +118,18 @@ -- We use the very first bit to decide which of these two encoding we use. -- (if we would make a sum type instead, it would take 2 extra words...) -- -newtype WordVec = WordVec Blob - -- deriving Show +-- About the instances: +-- +-- * the @Eq@ instance is strict: @x == y@ iff @toList x == toList y@. +-- For an equality which disregards trailing zeros, see 'eqExtZero' +-- +-- * the @Ord@ instance first compares the length, +-- then if the lengths are equal, compares the content lexicographically. +-- For a comparison which disregards the length, and lexicographically +-- compares the sequences extended with zeros, see 'cmpExtZero' +-- +newtype WordVec + = WordVec Blob -- | The \"shape\" of a dynamic word vector data Shape = Shape @@ -63,9 +141,10 @@ vecShape :: WordVec -> Shape vecShape = snd . vecShape' +-- | @vecShape' vec == (vecIsSmall vec , vecShape vec)@ vecShape' :: WordVec -> (Bool,Shape) vecShape' (WordVec blob) = (isSmall,shape) where - !h = blobHead blob + !h = Blob.head blob !h2 = shiftR h 1 !isSmall = (h .&. 1) == 0 shape = if isSmall @@ -74,11 +153,16 @@ mkShape :: Word64 -> Word64 -> Shape mkShape !x !y = Shape (fromIntegral x) (fromIntegral y) +-- | @True@ if the internal representation is the \"small\" one vecIsSmall :: WordVec -> Bool -vecIsSmall (WordVec blob) = (blobHead blob .&. 1) == 0 +vecIsSmall (WordVec !blob) = (Blob.head blob .&. 1) == 0 -vecLen, vecBits :: WordVec -> Int +-- | The length of the vector +vecLen :: WordVec -> Int vecLen = shapeLen . vecShape + +-- | The number of bits per element used to encode the vector +vecBits :: WordVec -> Int vecBits = shapeBits . vecShape -------------------------------------------------------------------------------- @@ -98,27 +182,54 @@ . showChar ' ' . shows (toList dynvec) +-- | The Eq instance is strict: @x == y@ iff @toList x == toList y@. +-- For an equality which disregards trailing zeros, see 'eqExtZero'. instance Eq WordVec where - (==) x y = (vecLen x == vecLen y) && (toList x == toList y) + (==) x y = eqStrict x y + -- (==) x y = (vecLen x == vecLen y) && (toList x == toList y) +-- | The Ord instance first compares the length, then if the lengths are equal, +-- compares the content lexicographically. For a different ordering, see 'cmpExtZero'. instance Ord WordVec where + compare x y = cmpStrict x y +{- compare x y = case compare (vecLen x) (vecLen y) of LT -> LT GT -> GT EQ -> compare (toList x) (toList y) +-} -------------------------------------------------------------------------------- --- * Empty vectors +-- * Empty vector, singleton empty :: WordVec empty = fromList [] null :: WordVec -> Bool -null v = (vecLen v == 0) +null (WordVec !blob) = + -- null v = (vecLen v == 0) + let !h = Blob.head blob + in (h .&. 0xf9 == 0) || (h .&. 0xffffffe1 == 1) + -- 0xf9 = 000 ... 00|11111001 + -- 0xffffffe1 = 111 ... 11|11100001 + +{- +null_naive :: WordVec -> Bool +null_naive v = (vecLen v == 0) +-} +singleton :: Word -> WordVec +singleton !x = fromListN 1 x [x] where + +isSingleton :: WordVec -> Maybe Word +isSingleton !v = case (vecLen v) of + 1 -> Just (head v) + _ -> Nothing + -------------------------------------------------------------------------------- -- * Indexing +-- | No boundary check is done. Indexing starts from 0. unsafeIndex :: Int -> WordVec -> Word unsafeIndex idx dynvec@(WordVec blob) = case isSmall of @@ -136,16 +247,86 @@ False -> extractSmallWord bits blob (32 + bits*idx) where (isSmall, Shape len bits) = vecShape' dynvec + +-------------------------------------------------------------------------------- +-- * Head, tail, etc +-- | Note: For the empty vector, @head@ returns 0 head :: WordVec -> Word -head dynvec@(WordVec blob) = - case vecIsSmall dynvec of - True -> extractSmallWord bits blob 8 - False -> extractSmallWord bits blob 32 +head dynvec@(WordVec blob) + | null dynvec = 0 + | otherwise = case vecIsSmall dynvec of + True -> extractSmallWord bits blob 8 + False -> extractSmallWord bits blob 32 where bits = vecBits dynvec +-- | Note: For the empty vector, @last@ returns 0 +last :: WordVec -> Word +last dynvec@(WordVec blob) + | len == 0 = 0 + | otherwise = case isSmall of + True -> extractSmallWord bits blob ( 8 + bits*(len-1)) + False -> extractSmallWord bits blob (32 + bits*(len-1)) + where + (isSmall, Shape len bits) = vecShape' dynvec + -------------------------------------------------------------------------------- + +-- | Note: For the empty vector, @tail@ returns (another) empty vector +tail :: WordVec -> WordVec +tail = tail_v2 + +-- | Prepends an element +cons :: Word -> WordVec -> WordVec +cons = cons_v2 + +-- | Appends an element +snoc :: WordVec -> Word -> WordVec +snoc = snoc_v2 + +uncons :: WordVec -> Maybe (Word, WordVec) +uncons = uncons_v2 + +concat :: WordVec -> WordVec -> WordVec +concat u v = fromList' (Shape (lu+lv) (max bu bv)) (toList u ++ toList v) where + Shape lu bu = vecShape u + Shape lv bv = vecShape v + +-------------------------------------------------------------------------------- + +foreign import ccall unsafe "vec_identity" c_vec_identity :: CFun11_ -- for testing +foreign import ccall unsafe "vec_tail" c_vec_tail :: CFun11_ +foreign import ccall unsafe "vec_head_tail" c_vec_head_tail :: CFun11 Word64 +foreign import ccall unsafe "vec_cons" c_vec_cons :: Word64 -> CFun11_ +foreign import ccall unsafe "vec_snoc" c_vec_snoc :: Word64 -> CFun11_ + +tail_v2 :: WordVec -> WordVec +tail_v2 (WordVec blob) = WordVec $ wrapCFun11_ c_vec_tail id blob + +cons_v2 :: Word -> WordVec -> WordVec +cons_v2 y vec@(WordVec blob) = WordVec $ wrapCFun11_ (c_vec_cons (fromIntegral y)) f blob where + f !n = max (n+2) worstcase + len = vecLen vec + worstcase = shiftR (32 + bitsNeededFor y * (len+1) + 63) 6 + -- it can happen that we cons (2^64-1) to a long vector of 4 bit numbers... + -- now it either fits in the old bits, in which case we need at most 1 new word + -- (maybe two, if we also switch from small header to big header at the same time???) + -- or does not, which is computed by @worstcase@ + +snoc_v2 :: WordVec -> Word -> WordVec +snoc_v2 vec@(WordVec blob) y = WordVec $ wrapCFun11_ (c_vec_snoc (fromIntegral y)) f blob where + f !n = max (n+2) worstcase + len = vecLen vec + worstcase = shiftR (32 + bitsNeededFor y * (len+1) + 63) 6 + +uncons_v2 :: WordVec -> Maybe (Word,WordVec) +uncons_v2 vec@(WordVec blob) = if null vec + then Nothing + else let (hd,tl) = wrapCFun11 c_vec_head_tail id blob + in Just (fromIntegral hd , WordVec tl) + +-------------------------------------------------------------------------------- -- * Conversion to\/from lists toList :: WordVec -> [Word] @@ -178,13 +359,13 @@ !elem = mask (this .|. shiftL that (64-bitOfs)) in elem : worker newOfs' (k-1) (shiftR that newOfs' : rest') [] -> error "WordVec/toList: FATAL ERROR! this should not happen" - --- | Another implementation of 'toList', for testing purposes only -toList_naive :: WordVec -> [Word] -toList_naive dynvec@(WordVec blob) = + +-- | @toRevList vec == reverse (toList vec)@, but should be faster (?) +toRevList :: WordVec -> [Word] +toRevList dynvec@(WordVec blob) = case isSmall of - True -> [ extractSmallWord bits blob ( 8 + bits*i) | i<-[0..len-1] ] - False -> [ extractSmallWord bits blob (32 + bits*i) | i<-[0..len-1] ] + True -> [ extractSmallWord bits blob ( 8 + bits*(len-i)) | i<-[1..len] ] + False -> [ extractSmallWord bits blob (32 + bits*(len-i)) | i<-[1..len] ] where (isSmall, Shape len bits) = vecShape' dynvec @@ -194,8 +375,17 @@ fromList [] = fromList' (Shape 0 4) [] fromList xs = fromList' (Shape l b) xs where l = length xs - b = bitsNeededFor (maximum xs) - + b = bitsNeededFor (L.maximum xs) + +-- | This is faster than 'fromList' +fromListN + :: Int -- ^ length + -> Word -- ^ maximum (or just an upper bound) + -> [Word] -- ^ elements + -> WordVec +fromListN len max = fromList' (Shape len (bitsNeededFor max)) + +-- | If you know the shape in advance, it\'s faster to use this function fromList' :: Shape -> [Word] -> WordVec fromList' (Shape len bits0) words | bits <= 16 && len <= 31 = WordVec $ mkBlob (mkHeader 0 2) 8 words @@ -227,8 +417,140 @@ in current' : worker (k-1) (shiftR this (64-bitOfs)) newOfs' rest -------------------------------------------------------------------------------- --- * Some more operations +-- * Specialized operations +-- +-- $spec +-- +-- These are are faster than the generic operations below, and should be preferred +-- to those. +-- +-------------------------------------------------------------------------------- +-- ** Specialized folds + +-- | Sum of the elements of the vector +sum :: WordVec -> Word +sum (WordVec blob) = fromIntegral $ wrapCFun10 c_vec_sum blob + +-- | Maximum of the elements of the vector +maximum :: WordVec -> Word +maximum (WordVec blob) = fromIntegral $ wrapCFun10 c_vec_max blob + +foreign import ccall unsafe "vec_sum" c_vec_sum :: CFun10 Word64 +foreign import ccall unsafe "vec_max" c_vec_max :: CFun10 Word64 + +-------------------------------------------------------------------------------- +-- ** Specialized \"zipping folds\" +-- + +foreign import ccall unsafe "vec_equal_strict" c_equal_strict :: CFun20 CInt +foreign import ccall unsafe "vec_equal_extzero" c_equal_extzero :: CFun20 CInt +foreign import ccall unsafe "vec_compare_strict" c_compare_strict :: CFun20 CInt +foreign import ccall unsafe "vec_compare_extzero" c_compare_extzero :: CFun20 CInt +foreign import ccall unsafe "vec_less_or_equal" c_less_or_equal :: CFun20 CInt +foreign import ccall unsafe "vec_partial_sums_less_or_equal" c_partial_sums_less_or_equal :: CFun20 CInt + +-- | Strict equality of vectors (same length, same content) +eqStrict :: WordVec -> WordVec -> Bool +eqStrict (WordVec blob1) (WordVec blob2) = (0 /= wrapCFun20 c_equal_strict blob1 blob2) + +-- | Equality of vectors extended with zeros to infinity +eqExtZero :: WordVec -> WordVec -> Bool +eqExtZero (WordVec blob1) (WordVec blob2) = (0 /= wrapCFun20 c_equal_extzero blob1 blob2) + +cintToOrdering :: CInt -> Ordering +cintToOrdering !k + | k < 0 = LT + | k > 0 = GT + | otherwise = EQ + +-- | Strict comparison of vectors (first compare the lengths; if the lengths are the same then compare lexicographically) +cmpStrict :: WordVec -> WordVec -> Ordering +cmpStrict (WordVec blob1) (WordVec blob2) = cintToOrdering $ wrapCFun20 c_compare_strict blob1 blob2 + +-- | Lexicographic ordering of vectors extended with zeros to infinity +cmpExtZero :: WordVec -> WordVec -> Ordering +cmpExtZero (WordVec blob1) (WordVec blob2) = cintToOrdering $ wrapCFun20 c_compare_extzero blob1 blob2 + +-- | Pointwise comparison of vectors extended with zeros to infinity +lessOrEqual :: WordVec -> WordVec -> Bool +lessOrEqual (WordVec blob1) (WordVec blob2) = (0 /= wrapCFun20 c_less_or_equal blob1 blob2) + +-- | Pointwise comparison of partial sums of vectors extended with zeros to infinity +-- +-- For example @[x1,x2,x3] <= [y1,y2,y3]@ iff (@x1 <=y1 && x1+x2 <= y1+y2 && x1+x2+x3 <= y1+y2+y3@). +-- +partialSumsLessOrEqual :: WordVec -> WordVec -> Bool +partialSumsLessOrEqual (WordVec blob1) (WordVec blob2) = + (0 /= wrapCFun20 c_partial_sums_less_or_equal blob1 blob2) + +-------------------------------------------------------------------------------- +-- ** Specialized zips +-- + +foreign import ccall unsafe "vec_add" c_vec_add :: CFun21_ +foreign import ccall unsafe "vec_sub_overflow" c_vec_sub_overflow :: CFun21 CInt + +-- | Pointwise addition of vectors. The shorter one is extended by zeros. +add :: WordVec -> WordVec -> WordVec +add vec1@(WordVec blob1) vec2@(WordVec blob2) = WordVec $ wrapCFun21_ c_vec_add f blob1 blob2 where + -- WARNING! memory allocation is _very_ tricky here! + -- worst case: we have a very long vector with 4 bits/elem, + -- and a very short vector with 64 bits/elem! + -- even @max b1 b2@ is not enough, because it can overflow... + f _ _ = 1 + shiftR ( (max b1 b2 + 4)*(max l1 l2) + 63 ) 6 + Shape !l1 !b1 = vecShape vec1 + Shape !l2 !b2 = vecShape vec2 + +-- | Pointwise subtraction of vectors. The shorter one is extended by zeros. +-- If any element would become negative, we return Nothing +subtract :: WordVec -> WordVec -> Maybe WordVec +subtract vec1@(WordVec blob1) vec2@(WordVec blob2) = + case (wrapCFun21 c_vec_sub_overflow f blob1 blob2) of + (0 , blob3) -> Just (WordVec blob3) + (_ , _ ) -> Nothing + where + f _ _ = 1 + shiftR ( (max b1 b2 + 4)*(max l1 l2) + 63 ) 6 + Shape !l1 !b1 = vecShape vec1 + Shape !l2 !b2 = vecShape vec2 + +-------------------------------------------------------------------------------- +-- ** Specialized maps + +foreign import ccall unsafe "vec_scale" c_vec_scale :: Word64 -> CFun11_ + +-- | Pointwise multiplication by a constant. +scale :: Word -> WordVec -> WordVec +scale s vec@(WordVec blob) = WordVec $ wrapCFun11_ (c_vec_scale (fromIntegral s)) f blob where + f _ = shiftR (32 + len*newbits + 63) 6 + Shape !len !bits = vecShape vec + bound = if s <= shiftL 1 (64-bits) + then (2^bits - 1) * s + else (2^64 - 1) + newbits = bitsNeededFor bound + +-------------------------------------------------------------------------------- +-- ** Specialized scans + +foreign import ccall unsafe "vec_partial_sums" c_vec_partial_sums :: CFun11 Word64 + +-- | @toList (partialSums vec) == tail (scanl (+) 0 $ toList vec)@ +partialSums :: WordVec -> WordVec +partialSums vec@(WordVec blob) = WordVec $ snd $ wrapCFun11 c_vec_partial_sums f blob where + f _ = shiftR (32 + len*newbits + 63) 6 + Shape !len !bits = vecShape vec + bound = if len <= shiftL 1 (64-bits) + then (2^bits - 1) * (fromIntegral len :: Word) -- worst case: @replicate N (2^bits-1)@ + else (2^64 - 1) + newbits = bitsNeededFor bound + +-------------------------------------------------------------------------------- +-- * Some generic operations + +-- | Left fold +fold :: (a -> Word -> a) -> a -> WordVec -> a +fold f x v = L.foldl' f x (toList v) + naiveMap :: (Word -> Word) -> WordVec -> WordVec naiveMap f u = fromList (map f $ toList u) @@ -239,33 +561,64 @@ l = vecLen vec bits = bitsNeededFor bound -concat :: WordVec -> WordVec -> WordVec -concat u v = fromList' (Shape (lu+lv) (max bu bv)) (toList u ++ toList v) where - Shape lu bu = vecShape u - Shape lv bv = vecShape v - naiveZipWith :: (Word -> Word -> Word) -> WordVec -> WordVec -> WordVec -naiveZipWith f u v = fromList $ zipWith f (toList u) (toList v) +naiveZipWith f u v = fromList $ L.zipWith f (toList u) (toList v) -- | If you have a (nearly sharp) upper bound to the result of your of function -- on your vector, zipping can be more efficient boundedZipWith :: Word -> (Word -> Word -> Word) -> WordVec -> WordVec -> WordVec -boundedZipWith bound f vec1 vec2 = fromList' (Shape l bits) $ zipWith f (toList vec1) (toList vec2) where +boundedZipWith bound f vec1 vec2 = fromList' (Shape l bits) $ L.zipWith f (toList vec1) (toList vec2) where l = min (vecLen vec1) (vecLen vec2) bits = bitsNeededFor bound + +listZipWith :: (Word -> Word -> a) -> WordVec -> WordVec -> [a] +listZipWith f u v = L.zipWith f (toList u) (toList v) -------------------------------------------------------------------------------- -- * Misc helpers +-- | Number of bits needed to encode a given number, rounded up to multiples of four bitsNeededFor :: Word -> Int -bitsNeededFor bound = ceilingLog2 (bound + 1) where -- for example, if maximum is 16, log2 = 4 but we need 5 bits +bitsNeededFor = bitsNeededForHs - -- | Smallest integer @k@ such that @2^k@ is larger or equal to @n@ - ceilingLog2 :: Word -> Int - ceilingLog2 0 = 0 - ceilingLog2 n = 1 + go (n-1) where - go 0 = -1 - go k = 1 + go (shiftR k 1) +-- | Number of bits needed to encode a given number +bitsNeededFor' :: Word -> Int +bitsNeededFor' = bitsNeededForHs' + +bitsNeededForHs :: Word -> Int +bitsNeededForHs = roundBits . bitsNeededForHs' + +bitsNeededForHs' :: Word -> Int +bitsNeededForHs' bound + | bound == 0 = 1 -- this is handled incorrectly by the formula below + | bound+1 == 0 = MACHINE_WORD_BITS -- and this handled incorrectly because of overflow + | otherwise = ceilingLog2 (bound + 1) -- for example, if maximum is 16, log2 = 4 but we need 5 bits + where + -- | Smallest integer @k@ such that @2^k@ is larger or equal to @n@ + ceilingLog2 :: Word -> Int + ceilingLog2 0 = 0 + ceilingLog2 n = 1 + go (n-1) where + go 0 = -1 + go k = 1 + go (shiftR k 1) + +{- + +-- apparently, the C implementation is _not_ faster... + +foreign import ccall unsafe "export_required_bits_not_rounded" export_required_bits_not_rounded :: Word64 -> CInt +foreign import ccall unsafe "export_required_bits" export_required_bits :: Word64 -> CInt + +bitsNeededForC :: Word -> Int +bitsNeededForC = fromIntegral . export_required_bits . fromIntegral + +bitsNeededForC' :: Word -> Int +bitsNeededForC' = fromIntegral . export_required_bits_not_rounded . fromIntegral +-} + +-- | We only allow multiples of 4. +roundBits :: Int -> Int +roundBits 0 = 4 +roundBits k = shiftL (shiftR (k+3) 2) 2 -------------------------------------------------------------------------------- diff --git a/test/TestSuite.hs b/test/TestSuite.hs --- a/test/TestSuite.hs +++ b/test/TestSuite.hs @@ -11,9 +11,18 @@ import qualified Tests.WordVec import qualified Tests.IntVec +-- import System.Random + -------------------------------------------------------------------------------- -main = defaultMain tests +{- +speedtest_main = do + setStdGen (mkStdGen 12345) -- for reproducible timing! + ... +-} + +main = do + defaultMain tests tests :: TestTree tests = testGroup "Tests" diff --git a/test/Tests/Blob.hs b/test/Tests/Blob.hs --- a/test/Tests/Blob.hs +++ b/test/Tests/Blob.hs @@ -6,28 +6,50 @@ -------------------------------------------------------------------------------- +import Data.Bits import Data.Word import Data.List as L import Data.Primitive.ByteArray -import Data.Vector.Compact.Blob +import Data.Vector.Compact.Blob as B import Test.Tasty import Test.Tasty.HUnit -------------------------------------------------------------------------------- +-- | An alternate implementation using 'testBit', for testing purposes only +extractSmallWord64_naive :: Int -> Blob -> Int -> Word64 +extractSmallWord64_naive n blob ofs = sum [ shiftL 1 i | i<-[0..n-1] , testBit blob (ofs+i) ] + +-------------------------------------------------------------------------------- + all_tests = tests_blobs tests_blobs = testGroup "unit tests for Blobs" - [ testCase "toList . fromList == id" $ forall_ w64_lists prop_from_to_list - , testCase "fromList . toList == id" $ forall_ blobs prop_to_from_blob - , testCase "toList vs. peekWord" $ forall_ blobs prop_tolist_vs_peek - , testCase "blobHead vs. list head" $ forall_ w64_lists prop_head_of_list - , testCase "eqBlob vs. naive" $ forall_ blob_pairs prop_eq_vs_naive - , testCase "fromBA . toBA == id" $ forall_ blobs prop_to_from_bytearray - , testCase "toBA . fromBA == pad" $ forall_ byte_lists prop_from_to_bytearray + [ testCase "toList . fromList == id" $ forall_ w64_lists prop_from_to_list + , testCase "fromList . toList == id" $ forall_ blobs prop_to_from_blob + , testCase "toList vs. peekWord" $ forall_ blobs prop_tolist_vs_peek + , testCase "blobHead vs. list head" $ forall_ w64_lists prop_head_of_list + , testCase "eqBlob vs. naive" $ forall_ blob_pairs prop_eq_vs_naive + , testCase "fromBA . toBA == id" $ forall_ blobs prop_to_from_bytearray + , testCase "toBA . fromBA == pad" $ forall_ byte_lists prop_from_to_bytearray + -- + , testCase "tail . cons == id" $ forall_ blobs prop_tail_cons + -- + , testCase "rotateL . rotateR == id" $ forall_ blobs prop_rotateL_rotateR + , testCase "rotateR . rotateL == id" $ forall_ blobs prop_rotateR_rotateL + , testCase "rotateL [1]" $ forall_ [0..130] prop_rotateL_one + , testCase "rotateR [1]" $ forall_ [0..130] prop_rotateR_one + , testCase "rotateL [1,0,0]" $ forall_ [0..530] prop_rotateL_one_blob3 + , testCase "rotateR [0,0,z]" $ forall_ [0..530] prop_rotateR_one_blob3 + -- + , testCase "shiftR_64 . shiftL_64 == id" $ forall_ blobs prop_shiftR64_shiftL64 + , testCase "shiftL_64 . shiftR_64 = ..." $ forall_ blobs prop_shiftL64_shiftR64 + , testCase "shiftR 64 == tail" $ forall_ blobs prop_shiftR64_is_tail + , testCase "shiftR . shiftL ~= id" $ forall_ blobs prop_shiftR_shiftL + , testCase "shiftL [1]" $ forall_ [0..530] prop_shiftL_one ] forall_ :: [a] -> (a -> Bool) -> Assertion @@ -36,6 +58,54 @@ -------------------------------------------------------------------------------- -- * inputs +w1,w2 :: Word64 +w1 = 0x1234567890abcdef +w2 = 0xefcdab9078563412 + +rndWords1, rndWords2 :: [Word64] + +rndWords1 = + [0x636bab7ce12a9f82,0xacb1651d835dea8a,0x53ee90f119788b9c,0xb3e5f34f1a7b14df,0x108afb5d5ef5420b + ,0xad8a4492183b5890,0x29338ff97c3c4b30,0xab38a94da16f8831,0x3eebedc58090704f,0xfb0de19c60305cfc + ,0xe1f9cc292fd9fd99,0x661cbeb165fc5369,0x45b9557954fa8197,0x5999769d181a0e58,0xbbbc1a04bf66f1a7 + ,0x8b14a985d9e28575,0xa59e5b55e5bd190f,0xc86dcfa51ec3f984,0x8ada7c8b2a057cc2,0x38085a65b744237d + ,0xcaf0303da21ecc16,0xc660be5561054445,0xeaf4360fc8aa9031,0xd60926fff096d2ab,0x1a03ac007c232799 + ,0xe584496fb0608fc0,0x09079c518d7206e3,0x4c5a70bebce35d84,0xf8200f7b4dae3d6a,0xb4a6c4c4b5d16a38 + ,0x5451bc4e7fb2cfeb,0x6525fd92075937b6,0x406eae19f78ec53e,0x6bf94e6694a523ec,0xfeb03e90e7faba04 + ,0xf10a0d43f2c28d4b,0xe44aa8f952a71ce4,0x6ef966e88c95e8cc,0x1f0a91de603823fe,0x540fc5e2212688c5 + ,0x98538a3191aa00c4,0x82a866243999920b,0x7c531773ea9a576e,0x683e1d213a11048e,0x001441c42aea5812 + ,0x278d5e7726567c30,0xe0af4e322898ba88,0x6693727353fd7cb0,0x87fc50f2604236e4,0x888c720e479ba9ba + ,0x9aad9422fa1cd5ce,0x77b7d3c79b10764a,0x693f8172f63598a7,0x61daefcdebed2f6e,0x999306f3808e557e + ,0x5d58c9725e45aed4,0x56aa23006a204e29,0x7b9f38d5f920cb7b,0xd1a79fbb76b795d0,0xe268fbe17c6672d3 + ,0xb8d8989b3d3dba70,0xa12c3645bb8cdf93,0x732cb315500e9ff6,0x45614a09f768ab04,0xd7a44a81f4a4626d + ] + +rndWords2 = + [0x888efd3a8e9d3d0f,0x3f7c8015929c4d0b,0x9dd1576dfee9144d,0x706e50d93edb988f,0x1bf0fbd22cb2fe42 + ,0xea3cb4e94f17df4a,0xf24bc3fb3eb11658,0xfdafd48ff4a6560e,0xf3fb721e97c24e28,0xe4189f3b29bab63c + ,0x9888e3eed7a10abe,0x3f655e3ac3f99dfe,0xbdaf179b0ab70e1b,0x933489e815343e8d,0xc552296865200a2b + ,0x0d7ee80f467eeb04,0xfd6b778e8babe925,0x1bb1255de3cd8786,0xc4a3b6f573ec6af1,0xd2dbf7d5b6a3be2e + ,0xa1e22b70a36c96f8,0xa617f27f72fe8ff2,0x1e5024e9159ff0ae,0x603f5a9c50986495,0xe2f17131fb59bd93 + ,0x8d01a921a59646b7,0xc85ae16a975b9e97,0x3ff90c245b9ad063,0x4889782b52890c9f,0xa79b24111b5558d9 + ,0x9bc09895641cdd04,0xda8a2359ad7d0335,0x4e925260435f69a0,0x28fd2829993ee9b7,0xb38af0536c569fee + ,0x78462beafef9e57c,0x2613166891ab270e,0x13775c45b6efb6c4,0x69921e49759bdda0,0x68bfc6f1643ff245 + ,0xf7dcc717f0549f2d,0x423e6b5bbfa2a8b6,0x5555fe6cc05d5519,0xc90f85f121c6adf0,0x68ec14dc41c1174d + ,0x231473484795b255,0xa13c37dfbf1bfc44,0x5e3831b09836a34b,0x93fb41b63eedc8f8,0x07b3f59a37165884 + ,0x0b760cb359a1b6af,0x478b77fb0f54786e,0x024b1a3ddc880e54,0x3d565cd1272a3d81,0x9984a2d9ab5a4741 + ,0x8234b289e02f7aec,0x24e65cdf39319ded,0x7d3fd5f3b4f6b5a2,0x9289ff5cc281dcbd,0xb189bc7fee3c6c7c + ,0x14f628aa402b4ab0,0x9dbe78e5f46e9e33,0xe2eeee94e9b874c3,0x6c236ec875b93341,0xf0429b07692c691b + ] + +{- +-- generate random words +import Data.List ; import Control.Monad ; import System.Random ; import Data.Word ; import Text.Printf +main = do + list <- replicateM 65 randomIO :: IO [Word64] + putStrLn $ "[" ++ intercalate "," (map (printf "0x%016x") list) ++ "]" +-} + +-------------------------------------------------------------------------------- + w64_lists_of_length :: Int -> [[Word64]] w64_lists_of_length ni = [ [ i | i<-[1..n] ] @@ -43,6 +113,10 @@ , [ i | i<-[2^32-1-n .. 2^32-1 ] ] , [ i | i<-[2^32-1 .. 2^32+n-2 ] ] , [ i | i<-[2^64-1-n .. 2^64-1 ] ] + , take ni rndWords1 + , take ni rndWords2 + , take ni $ reverse rndWords1 + , take ni $ reverse rndWords2 ] where n = fromIntegral ni :: Word64 @@ -58,30 +132,76 @@ eqBlob_naive b1 b2 = blobToWordList b1 == blobToWordList b2 -baToList :: ByteArray -> [Word8] -baToList = foldrByteArray (:) [] +baToByteList :: ByteArray -> [Word8] +baToByteList = foldrByteArray (:) [] pad :: [Word8] -> [Word8] pad xs = xs ++ replicate (len8-len) 0 where len = length xs len8 = 8 * (div (len+7) 8) +local_longZipWith :: Integral a => (a -> a -> b) -> [a] -> [a] -> [b] +local_longZipWith f (x:xs) (y:ys) = f x y : local_longZipWith f xs ys +local_longZipWith f xs [] = zipWith f xs (repeat 0) +local_longZipWith f [] ys = zipWith f (repeat 0) ys + +eqWithZeros :: Blob -> Blob -> Bool +eqWithZeros a b = and $ local_longZipWith (==) (blobToWordList a) (blobToWordList b) + -------------------------------------------------------------------------------- -- * properties prop_from_to_list list = blobToWordList (blobFromWordList list) == list prop_to_from_blob blob = blobFromWordList (blobToWordList blob) == blob -prop_tolist_vs_peek blob = [ peekWord blob i | i<-[0..n-1] ] == blobToWordList blob where +prop_tolist_vs_peek blob = [ B.indexWord blob i | i<-[0..n-1] ] == blobToWordList blob where n = blobSizeInWords blob -prop_head_of_list list = blobHead (blobFromWordList list) == L.head list +prop_head_of_list list = B.head (blobFromWordList list) == L.head list prop_eq_vs_naive (b1,b2) = eqBlob b1 b2 == eqBlob_naive b1 b2 prop_to_from_bytearray blob = blobFromByteArray (blobToByteArray blob) == blob -prop_from_to_bytearray list = baToList (blobToByteArray (blobFromByteArray ba)) == pad list where +prop_from_to_bytearray list = baToByteList (blobToByteArray (blobFromByteArray ba)) == pad list where ba = byteArrayFromList list +-------------------------------------------------------------------------------- + +prop_tail_cons blob = B.tail (B.consWord 0x1234567890abcdef blob) == blob +prop_shiftR64_is_tail blob = (shiftR blob 64) == B.tail blob +prop_shiftR64_shiftL64 blob = shiftR (shiftL blob 64) 64 == blob +prop_shiftL64_shiftR64 blob = shiftL (shiftR blob 64) 64 == B.consWord 0 (B.tail blob) + +prop_shiftR_shiftL blob = and [ shiftR (shiftL blob i) i `eqWithZeros` blob | i<-[0..201] ] + +prop_rotateL_rotateR blob = and [ rotateL (rotateR blob i) i == blob | i<-[0..201] ] +prop_rotateR_rotateL blob = and [ rotateR (rotateL blob i) i == blob | i<-[0..201] ] + +prop_rotateL_one i = rotateL (Blob1 1) i == Blob1 (rotateL 1 i) +prop_rotateR_one i = rotateR (Blob1 1) i == Blob1 (rotateR 1 i) + +prop_rotateL_one_blob3 i = x == y where + x = blobToWordList $ rotateL (blobFromWordList [1,0,0]) i + y = case divMod i 64 of + (q0,r) -> case (mod q0 3, r) of + (0,r) -> [ shiftL 1 r , 0 , 0 ] + (1,r) -> [ 0 , shiftL 1 r , 0 ] + (2,r) -> [ 0 , 0 , shiftL 1 r ] + +prop_rotateR_one_blob3 i = x == y where + x = blobToWordList $ rotateR (blobFromWordList [0,0,z]) i + z = shiftL 1 63 + y = case divMod i 64 of + (q0,r) -> case (mod q0 3, r) of + (0,r) -> [ 0 , 0 , shiftR z r ] + (1,r) -> [ 0 , shiftR z r , 0 ] + (2,r) -> [ shiftR z r , 0 , 0 ] + +prop_shiftL_one i = x == y where + x = blobToWordList $ shiftL (Blob1 1) i + y = case divMod i 64 of + (q,r) -> replicate q 0 ++ [shiftL 1 r] + + -------------------------------------------------------------------------------- diff --git a/test/Tests/IntVec.hs b/test/Tests/IntVec.hs --- a/test/Tests/IntVec.hs +++ b/test/Tests/IntVec.hs @@ -6,8 +6,10 @@ -------------------------------------------------------------------------------- +import Control.Monad ( liftM ) + import Data.Int -import Data.List as L +import qualified Data.List as L import Data.Vector.Compact.IntVec as V @@ -24,25 +26,55 @@ -------------------------------------------------------------------------------- +-- | For testing purposes only +uncons_naive :: IntVec -> Maybe (Int,IntVec) +uncons_naive vec = if V.null vec + then Nothing + else Just (V.head vec, V.tail vec) + +-------------------------------------------------------------------------------- + all_tests = testGroup "unit tests for IntVec-s" - [ tests_small + [ tests_unit + , tests_small , tests_bighead ] +tests_unit = testGroup "misc unit tests" + [ testCase "equality with different bit sizes" $ assertBool "failed" $ check_eq_bitsizes + , testCase "head of empty == 0" $ assertBool "failed" $ (V.head V.empty == 0) + , testCase "last of empty == 0" $ assertBool "failed" $ (V.last V.empty == 0) + , testCase "tail of empty == empty" $ assertBool "failed" $ (V.tail V.empty == V.empty) + ] + tests_small = testGroup "unit tests for small dynamic int vectors" - [ testCase "toList . fromList == id" $ forall_ small_Lists prop_from_to_list - , testCase "fromList . toList == id" $ forall_ small_Vecs prop_to_from_vec - , testCase "fromList vs. indexing" $ forall_ small_Lists prop_fromlist_vs_index - , testCase "vec head vs. list head" $ forall_ small_NELists prop_head_of_list - , testCase "head vs. indexing" $ forall_ small_NEVecs prop_head_vs_index + [ testCase "toList . fromList == id" $ forall_ small_Lists prop_from_to_list + , testCase "fromList . toList == id" $ forall_ small_Vecs prop_to_from_vec + , testCase "fromList vs. indexing" $ forall_ small_Lists prop_fromlist_vs_index + , testCase "toRevList == reverse . toList" $ forall_ small_Vecs prop_toRevList + , testCase "vec head vs. list head" $ forall_ small_NELists prop_head_of_list + , testCase "vec last vs. list last" $ forall_ small_NELists prop_last_of_list + , testCase "head vs. indexing" $ forall_ small_NEVecs prop_head_vs_index + , testCase "cons . uncons == id" $ forall_ small_NEVecs prop_cons_uncons + , testCase "uncons . cons == id" $ forall_ small_cons prop_uncons_cons + , testCase "uncons vs. naive" $ forall_ small_Vecs prop_uncons_vs_naive + , testCase "uncons vs. list" $ forall_ small_Vecs prop_uncons_vs_list + , testCase "cons vs. list" $ forall_ small_cons prop_cons_vs_list ] tests_bighead = testGroup "unit tests for small dynamic int vectors with big heads" - [ testCase "toList . fromList == id" $ forall_ bighead_Lists prop_from_to_list - , testCase "fromList . toList == id" $ forall_ bighead_Vecs prop_to_from_vec - , testCase "fromList vs. indexing" $ forall_ bighead_Lists prop_fromlist_vs_index - , testCase "vec head vs. list head" $ forall_ bighead_NELists prop_head_of_list - , testCase "head vs. indexing" $ forall_ bighead_NEVecs prop_head_vs_index + [ testCase "toList . fromList == id" $ forall_ bighead_Lists prop_from_to_list + , testCase "fromList . toList == id" $ forall_ bighead_Vecs prop_to_from_vec + , testCase "fromList vs. indexing" $ forall_ bighead_Lists prop_fromlist_vs_index + , testCase "toRevList == reverse . toList" $ forall_ bighead_Vecs prop_toRevList + , testCase "vec head vs. list head" $ forall_ bighead_NELists prop_head_of_list + , testCase "vec last vs. list last" $ forall_ bighead_NELists prop_last_of_list + , testCase "head vs. indexing" $ forall_ bighead_NEVecs prop_head_vs_index + , testCase "cons . uncons == id" $ forall_ bighead_NEVecs prop_cons_uncons + , testCase "uncons . cons == id" $ forall_ bighead_cons prop_uncons_cons + , testCase "uncons vs. naive" $ forall_ bighead_Vecs prop_uncons_vs_naive + , testCase "uncons vs. list" $ forall_ bighead_Vecs prop_uncons_vs_list + , testCase "cons vs. list" $ forall_ bighead_cons prop_cons_vs_list ] forall_ :: [a] -> (a -> Bool) -> Assertion @@ -68,6 +100,9 @@ small_NEVecs :: [NEVec] small_NEVecs = [ NEVec (V.fromList xs) | NEList xs <- small_NELists ] +small_cons :: [(Int,Vec)] +small_cons = [ (x, Vec (V.fromList xs)) | NEList (x:xs) <- small_NELists ] + -------------------------------------------------------------------------------- add_bighead :: List -> [List] @@ -84,7 +119,17 @@ bighead_NELists = [ NEList xs | List xs <- bighead_Lists ] :: [NEList] bighead_NEVecs = [ NEVec v | Vec v <- bighead_Vecs ] :: [NEVec] +bighead_cons = [ (x, Vec (V.fromList xs)) | NEList (x:xs) <- bighead_NELists ] + -------------------------------------------------------------------------------- + +check_eq_bitsizes = and + [ V.fromList [k] == V.fromList' (1,(-2^b,2^b)) [k] + | k<-[-16..15] + , b<-[4..31] + ] + +-------------------------------------------------------------------------------- -- * properties prop_from_to_list (List list) = V.toList (V.fromList list) == list @@ -94,7 +139,21 @@ vec = V.fromList list n = V.vecLen vec +prop_toRevList (Vec vec) = V.toRevList vec == reverse (V.toList vec) + prop_head_of_list (NEList list) = V.head (V.fromList list) == L.head list +prop_last_of_list (NEList list) = V.last (V.fromList list) == L.last list prop_head_vs_index (NEVec vec ) = V.head vec == unsafeIndex 0 vec + +prop_cons_uncons (NEVec vec) = liftM (uncurry V.cons) (V.uncons vec) == Just vec +prop_uncons_cons (w,Vec vec) = V.uncons (V.cons w vec) == Just (w,vec) +prop_uncons_vs_naive (Vec vec) = V.uncons vec == uncons_naive vec + +prop_uncons_vs_list (Vec vec) = unconsToList (V.uncons vec) == L.uncons (V.toList vec) +prop_cons_vs_list (w,Vec vec) = V.toList (V.cons w vec) == w : (V.toList vec) + +unconsToList mb = case mb of + Nothing -> Nothing + Just (i,vec) -> Just (i, V.toList vec) -------------------------------------------------------------------------------- diff --git a/test/Tests/WordVec.hs b/test/Tests/WordVec.hs --- a/test/Tests/WordVec.hs +++ b/test/Tests/WordVec.hs @@ -1,50 +1,322 @@ -- | Tests for dynamic word vectors -{-# LANGUAGE CPP,BangPatterns #-} +{-# LANGUAGE CPP, BangPatterns, ForeignFunctionInterface #-} module Tests.WordVec where -------------------------------------------------------------------------------- +import Control.Monad + import Data.Word +import Data.Bits +import Data.Maybe import Data.List as L import Data.Vector.Compact.WordVec as V +import Data.Vector.Compact.Blob as B +import Foreign.C.Types + import Test.Tasty import Test.Tasty.HUnit +import System.Random +import System.IO.Unsafe as Unsafe + +-------------------------------------------------------------------------------- + + +-- the fuck, we are almost writing 2020 here, why am i doing this?!.... #ifdef x86_64_HOST_ARCH arch_bits = 64 #elif i386_HOST_ARCH arch_bits = 32 +#elif i686_HOST_ARCH +arch_bits = 32 +#elif aarch64_HOST_ARCH +arch_bits = 64 #else arch_bits = 32 #endif -------------------------------------------------------------------------------- +-- helper functions -all_tests = testGroup "unit tests for WordVec-s" - [ tests_small +listMax :: [Word] -> Word +listMax [] = 0 +listMax xs = L.maximum xs + +listLongZipWith :: (Word -> Word -> a) -> [Word] -> [Word] -> [a] +listLongZipWith f = go where + go (x:xs) (y:ys) = f x y : go xs ys + go (x:xs) [] = f x 0 : go xs [] + go [] (y:ys) = f 0 y : go [] ys + go [] [] = [] + +-------------------------------------------------------------------------------- + +cmpStrict_naive :: WordVec -> WordVec -> Ordering +cmpStrict_naive x y = case compare (vecLen x) (vecLen y) of + LT -> LT + GT -> GT + EQ -> compare (toList x) (toList y) + +cmpExtZero_naive :: WordVec -> WordVec -> Ordering +cmpExtZero_naive x y = go (toList x) (toList y) where + go (x:xs) (y:ys) = case compare x y of + LT -> LT + GT -> GT + EQ -> go xs ys + go (x:xs) [] = go (x:xs) [0] + go [] (y:ys) = go [0] (y:ys) + go [] [] = EQ + +-------------------------------------------------------------------------------- + +bitsNeededForHs :: Word -> Int +bitsNeededForHs = roundBits . bitsNeededForHs' + +bitsNeededForHs' :: Word -> Int +bitsNeededForHs' bound + | bound == 0 = 1 -- this is handled incorrectly by the formula below + | bound+1 == 0 = arch_bits -- MACHINE_WORD_BITS -- and this handled incorrectly because of overflow + | otherwise = ceilingLog2 (bound + 1) -- for example, if maximum is 16, log2 = 4 but we need 5 bits + where + -- | Smallest integer @k@ such that @2^k@ is larger or equal to @n@ + ceilingLog2 :: Word -> Int + ceilingLog2 0 = 0 + ceilingLog2 n = 1 + go (n-1) where + go 0 = -1 + go k = 1 + go (shiftR k 1) + +-- apparently, the C implementation is _not_ faster... +bitsNeededForC :: Word -> Int +bitsNeededForC = fromIntegral . export_required_bits . fromIntegral + +bitsNeededForC' :: Word -> Int +bitsNeededForC' = fromIntegral . export_required_bits_not_rounded . fromIntegral + +foreign import ccall unsafe "export_required_bits_not_rounded" export_required_bits_not_rounded :: Word64 -> CInt +foreign import ccall unsafe "export_required_bits" export_required_bits :: Word64 -> CInt + +-------------------------------------------------------------------------------- +-- * Naive, reference implementations + +-- | Another implementation of 'toList', for testing purposes only +toList_extract :: WordVec -> [Word] +toList_extract dynvec@(WordVec blob) = + case isSmall of + True -> [ B.extractSmallWord bits blob ( 8 + bits*i) | i<-[0..len-1] ] + False -> [ B.extractSmallWord bits blob (32 + bits*i) | i<-[0..len-1] ] + where + (isSmall, Shape len bits) = vecShape' dynvec + +-------------------------------------------------------------------------------- + +tail_v1 :: WordVec -> WordVec +tail_v1 dynvec + | len == 0 = empty + | otherwise = fromList' (Shape (len-1) bits) (L.tail $ toList dynvec) + where + (Shape len bits) = vecShape dynvec + +uncons_v1 :: WordVec -> Maybe (Word,WordVec) +uncons_v1 vec + | len == 0 = Nothing + | otherwise = Just $ case toList vec of { (w:ws) -> (w , fromList' (Shape (len-1) bits) ws) } + where + (Shape len bits) = vecShape vec + +cons_v1 :: Word -> WordVec -> WordVec +cons_v1 w vec = fromList' shape' (w : toList vec) where + (Shape len bits) = vecShape vec + bits' = max bits (bitsNeededFor w) + shape' = Shape (len+1) bits' + +snoc_v1 :: WordVec -> Word -> WordVec +snoc_v1 vec w = fromList' shape' (toList vec ++ [w]) where + (Shape len bits) = vecShape vec + bits' = max bits (bitsNeededFor w) + shape' = Shape (len+1) bits' + +tail_v2 = V.tail +cons_v2 = V.cons +snoc_v2 = V.snoc +uncons_v2 = V.uncons + +-------------------------------------------------------------------------------- + +null_naive :: WordVec -> Bool +null_naive v = (vecLen v == 0) + +tail_naive :: WordVec -> WordVec +tail_naive vec = if V.null vec + then empty + else fromList $ L.tail $ V.toList vec + +-- | For testing purposes only +uncons_naive :: WordVec -> Maybe (Word,WordVec) +uncons_naive vec = if V.null vec + then Nothing + else Just (V.head vec, tail_naive vec) + +add_naive :: WordVec -> WordVec -> WordVec +add_naive vec1 vec2 = V.fromList $ listLongZipWith (+) (V.toList vec1) (V.toList vec2) + +sub_naive :: WordVec -> WordVec -> Maybe WordVec +sub_naive vec1 vec2 = case and (listLongZipWith (>=) (V.toList vec1) (V.toList vec2)) of + True -> Just $ V.fromList $ listLongZipWith (-) (V.toList vec1) (V.toList vec2) + False -> Nothing + +scale_naive :: Word -> WordVec -> WordVec +scale_naive s = V.fromList . L.map (*s) . V.toList + +partialSums_naive :: WordVec -> WordVec +partialSums_naive = V.fromList . L.tail . L.scanl' (+) 0 . V.toList + +-------------------------------------------------------------------------------- + +all_tests = testGroup "tests for WordVec-s" + [ tests_unit + , tests_small + , tests_rnd , tests_bighead ] +tests_unit = testGroup "misc unit tests" + [ testCase "equality with different bit sizes" $ assertBool "failed" $ check_eq_bitsizes + , testCase "head of empty == 0" $ assertBool "failed" $ (V.head V.empty == 0) + , testCase "last of empty == 0" $ assertBool "failed" $ (V.last V.empty == 0) + , testCase "tail of empty == empty" $ assertBool "failed" $ (V.tail V.empty == V.empty) + , testCase "bitsNeededFor C vs. ref" $ forall_ around_powers_of_two (\k -> bitsNeededForHs k == bitsNeededForC k) + , testCase "bitsNeededFor' C vs. ref" $ forall_ around_powers_of_two (\k -> bitsNeededForHs' k == bitsNeededForC' k) + , testCase "cons_v2 crash (left-shift) is fixed" $ assertBool "failed" $ cons_v2_crash + ] + +-- [ k | k<-around_powers_of_two, bitsNeededFor' k /= bitsNeededForReference' k ] + tests_small = testGroup "unit tests for small dynamic word vectors" - [ testCase "toList . fromList == id" $ forall_ small_Lists prop_from_to_list - , testCase "fromList . toList == id" $ forall_ small_Vecs prop_to_from_vec - , testCase "fromList vs. indexing" $ forall_ small_Lists prop_fromlist_vs_index - , testCase "toList vs. naive" $ forall_ small_Vecs prop_tolist_vs_naive - , testCase "vec head vs. list head" $ forall_ small_NELists prop_head_of_list - , testCase "head vs. indexing" $ forall_ small_NEVecs prop_head_vs_index + [ testGroup "conversion, basic operations (small)" + [ testCase "toList . fromList == id" $ forall_ small_Lists prop_from_to_list + , testCase "fromList . toList == id" $ forall_ small_Vecs prop_to_from_vec + , testCase "fromList vs. indexing" $ forall_ small_Lists prop_fromlist_vs_index + , testCase "toList vs. naive" $ forall_ small_Vecs prop_tolist_vs_naive + , testCase "toRevList == reverse . toList" $ forall_ small_Vecs prop_toRevList + , testCase "vec head vs. list head" $ forall_ small_NELists prop_head_of_list + , testCase "vec last vs. list last" $ forall_ small_NELists prop_last_of_list + , testCase "vec tail vs. list tail" $ forall_ small_NEVecs prop_tail_of_list + , testCase "tail_v1 vs. tail_v2" $ forall_ small_Vecs prop_tail_v1_vs_v2 + , testCase "cons_v1 vs. cons_v2" $ forall_ small_cons prop_cons_v1_vs_v2 + , testCase "snoc_v1 vs. snoc_v2" $ forall_ small_cons prop_snoc_v1_vs_v2 + , testCase "uncons_v1 vs. uncons_v2" $ forall_ small_Vecs prop_uncons_v1_vs_v2 + , testCase "head vs. indexing" $ forall_ small_NEVecs prop_head_vs_index + , testCase "cons . uncons == id" $ forall_ small_NEVecs prop_cons_uncons + , testCase "uncons . cons == id" $ forall_ small_cons prop_uncons_cons + , testCase "uncons vs. naive" $ forall_ small_Vecs prop_uncons_vs_naive + , testCase "uncons vs. list" $ forall_ small_Vecs prop_uncons_vs_list + , testCase "cons vs. list" $ forall_ small_cons prop_cons_vs_list + , testCase "snoc vs. list" $ forall_ small_cons prop_snoc_vs_list + ] + , testGroup "\"advanced\" operations (small)" + [ testCase "sum vs. list" $ forall_ small_Vecs prop_sum_vs_list + , testCase "max vs. list" $ forall_ small_Vecs prop_max_vs_list + , testCase "strict equality vs. list" $ forall_ small_pairs prop_strict_eq_vs_list + , testCase "strict comparison vs. list" $ forall_ small_pairs prop_strict_cmp_vs_list + , testCase "ext0 equality vs. list" $ forall_ small_pairs prop_ext0_eq_vs_list + , testCase "ext0 comparison vs. list" $ forall_ small_pairs prop_ext0_cmp_vs_list + , testCase "less or equal vs. list" $ forall_ small_pairs prop_less_or_equal_vs_list + , testCase "add vs. naive" $ forall_ small_pairs prop_add_vs_naive + , testCase "sub vs. naive" $ forall_ small_pairs prop_sub_vs_naive + , testCase "add is commutative" $ forall_ small_pairs prop_add_commutative + , testCase "partial sums vs. naive" $ forall_ small_Vecs prop_psums_vs_naive + , testCase "scale 3 vs. naive" $ forall_ small_Vecs (prop_scale_vs_naive 3) + , testCase "scale 14 vs. naive" $ forall_ small_Vecs (prop_scale_vs_naive 14) + , testCase "scale 254 vs. naive" $ forall_ small_Vecs (prop_scale_vs_naive 254) + , testCase "scale 65000 vs. naive" $ forall_ small_Vecs (prop_scale_vs_naive 65000) + ] ] +tests_rnd = testGroup "tests for random dynamic word vectors" + [ testGroup "conversion, basic operations (random)" + [ testCase "toList . fromList == id" $ forall_ rnd_Lists prop_from_to_list + , testCase "fromList . toList == id" $ forall_ rnd_Vecs prop_to_from_vec + , testCase "fromList vs. indexing" $ forall_ rnd_Lists prop_fromlist_vs_index + , testCase "toList vs. naive" $ forall_ rnd_Vecs prop_tolist_vs_naive + , testCase "toRevList == reverse . toList" $ forall_ rnd_Vecs prop_toRevList + , testCase "vec head vs. list head" $ forall_ rnd_NELists prop_head_of_list + , testCase "vec last vs. list last" $ forall_ rnd_NELists prop_last_of_list + , testCase "vec tail vs. list tail" $ forall_ rnd_NEVecs prop_tail_of_list + , testCase "tail_v1 vs. tail_v2" $ forall_ rnd_Vecs prop_tail_v1_vs_v2 + , testCase "cons_v1 vs. cons_v2" $ forall_ rnd_cons prop_cons_v1_vs_v2 + , testCase "snoc_v1 vs. snoc_v2" $ forall_ rnd_cons prop_snoc_v1_vs_v2 + , testCase "uncons_v1 vs. uncons_v2" $ forall_ rnd_Vecs prop_uncons_v1_vs_v2 + , testCase "head vs. indexing" $ forall_ rnd_NEVecs prop_head_vs_index + , testCase "cons . uncons == id" $ forall_ rnd_NEVecs prop_cons_uncons + , testCase "uncons . cons == id" $ forall_ rnd_cons prop_uncons_cons + , testCase "uncons vs. naive" $ forall_ rnd_Vecs prop_uncons_vs_naive + , testCase "uncons vs. list" $ forall_ rnd_Vecs prop_uncons_vs_list + , testCase "cons vs. list" $ forall_ rnd_cons prop_cons_vs_list + , testCase "snoc vs. list" $ forall_ rnd_cons prop_snoc_vs_list + ] + , testGroup "\"advanced\" operations (random)" + [ testCase "sum vs. list" $ forall_ rnd_Vecs prop_sum_vs_list + , testCase "max vs. list" $ forall_ rnd_Vecs prop_max_vs_list + , testCase "strict equality vs. list" $ forall_ rnd_pairs prop_strict_eq_vs_list + , testCase "strict comparison vs. list" $ forall_ rnd_pairs prop_strict_cmp_vs_list + , testCase "ext0 equality vs. list" $ forall_ rnd_pairs prop_ext0_eq_vs_list + , testCase "ext0 comparison vs. list " $ forall_ rnd_pairs prop_ext0_cmp_vs_list + , testCase "less or equal vs. list" $ forall_ rnd_pairs prop_less_or_equal_vs_list + , testCase "add vs. naive" $ forall_ rnd_pairs prop_add_vs_naive + , testCase "sub vs. naive" $ forall_ rnd_pairs prop_sub_vs_naive + , testCase "add is commutative" $ forall_ rnd_pairs prop_add_commutative + , testCase "partial sums vs. naive" $ forall_ rnd_Vecs prop_psums_vs_naive + , testCase "scale 3 vs. naive" $ forall_ rnd_Vecs (prop_scale_vs_naive 3) + , testCase "scale 14 vs. naive" $ forall_ rnd_Vecs (prop_scale_vs_naive 14) + , testCase "scale 254 vs. naive" $ forall_ rnd_Vecs (prop_scale_vs_naive 254) + , testCase "scale 65000 vs. naive" $ forall_ rnd_Vecs (prop_scale_vs_naive 65000) + ] + ] + tests_bighead = testGroup "unit tests for small dynamic word vectors with big heads" - [ testCase "toList . fromList == id" $ forall_ bighead_Lists prop_from_to_list - , testCase "fromList . toList == id" $ forall_ bighead_Vecs prop_to_from_vec - , testCase "fromList vs. indexing" $ forall_ bighead_Lists prop_fromlist_vs_index - , testCase "toList vs. naive" $ forall_ bighead_Vecs prop_tolist_vs_naive - , testCase "vec head vs. list head" $ forall_ bighead_NELists prop_head_of_list - , testCase "head vs. indexing" $ forall_ bighead_NEVecs prop_head_vs_index + [ testGroup "conversion, basic operations (big head)" + [ testCase "toList . fromList == id" $ forall_ bighead_Lists prop_from_to_list + , testCase "fromList . toList == id" $ forall_ bighead_Vecs prop_to_from_vec + , testCase "fromList vs. indexing" $ forall_ bighead_Lists prop_fromlist_vs_index + , testCase "toList vs. naive" $ forall_ bighead_Vecs prop_tolist_vs_naive + , testCase "toRevList == reverse . toList" $ forall_ bighead_Vecs prop_toRevList + , testCase "vec head vs. list head" $ forall_ bighead_NELists prop_head_of_list + , testCase "vec last vs. list last" $ forall_ bighead_NELists prop_last_of_list + , testCase "vec tail vs. list tail" $ forall_ bighead_NEVecs prop_tail_of_list + , testCase "tail_v1 vs. tail_v2" $ forall_ bighead_Vecs prop_tail_v1_vs_v2 + , testCase "cons_v1 vs. cons_v2" $ forall_ bighead_cons prop_cons_v1_vs_v2 + , testCase "snoc_v1 vs. snoc_v2" $ forall_ bighead_cons prop_snoc_v1_vs_v2 + , testCase "uncons_v1 vs. uncons_v2" $ forall_ bighead_Vecs prop_uncons_v1_vs_v2 + , testCase "head vs. indexing" $ forall_ bighead_NEVecs prop_head_vs_index + , testCase "cons . uncons == id" $ forall_ bighead_NEVecs prop_cons_uncons + , testCase "uncons . cons == id" $ forall_ bighead_cons prop_uncons_cons + , testCase "uncons vs. naive" $ forall_ bighead_Vecs prop_uncons_vs_naive + , testCase "uncons vs. list" $ forall_ bighead_Vecs prop_uncons_vs_list + , testCase "cons vs. list" $ forall_ bighead_cons prop_cons_vs_list + , testCase "snoc vs. list" $ forall_ bighead_cons prop_snoc_vs_list + ] + , testGroup "\"advanced\" operations (big head)" + [ testCase "sum vs. list" $ forall_ bighead_Vecs prop_sum_vs_list + , testCase "max vs. list" $ forall_ bighead_Vecs prop_max_vs_list + , testCase "strict equality vs. list" $ forall_ bighead_pairs prop_strict_eq_vs_list + , testCase "strict comparison vs. list" $ forall_ bighead_pairs prop_strict_cmp_vs_list + , testCase "ext0 equality vs. list" $ forall_ bighead_pairs prop_ext0_eq_vs_list + , testCase "ext0 comparison vs. list " $ forall_ bighead_pairs prop_ext0_cmp_vs_list + , testCase "less or equal vs. list" $ forall_ bighead_pairs prop_less_or_equal_vs_list + , testCase "add vs. naive" $ forall_ bighead_pairs prop_add_vs_naive + , testCase "sub vs. naive" $ forall_ bighead_pairs prop_sub_vs_naive + , testCase "add is commutative" $ forall_ bighead_pairs prop_add_commutative + , testCase "partial sums vs. naive" $ forall_ bighead_Vecs prop_psums_vs_naive + , testCase "scale 3 vs. naive" $ forall_ bighead_Vecs (prop_scale_vs_naive 3) + , testCase "scale 14 vs. naive" $ forall_ bighead_Vecs (prop_scale_vs_naive 14) + , testCase "scale 254 vs. naive" $ forall_ bighead_Vecs (prop_scale_vs_naive 254) + , testCase "scale 65000 vs. naive" $ forall_ bighead_Vecs (prop_scale_vs_naive 65000) + ] ] forall_ :: [a] -> (a -> Bool) -> Assertion @@ -52,6 +324,19 @@ -------------------------------------------------------------------------------- -- * inputs + +randomSublistIO :: Double -> [a] -> IO [a] +randomSublistIO prob xs = catMaybes <$> mapM f xs where + f !x = do + s <- randomRIO (0,1) + if (s < prob) + then return (Just x) + else return Nothing + +{-# NOINLINE randomSublist #-} +randomSublist :: Double -> [a] -> [a] +randomSublist prob xs = Unsafe.unsafePerformIO (randomSublistIO prob xs) + newtype List = List [Word] deriving Show newtype NEList = NEList [Word] deriving Show @@ -70,32 +355,171 @@ small_NEVecs :: [NEVec] small_NEVecs = [ NEVec (V.fromList xs) | NEList xs <- small_NELists ] +small_cons :: [(Word,Vec)] +small_cons = [ (x, Vec (V.fromList xs)) | NEList (x:xs) <- small_NELists ] + +smaller_Vecs :: [Vec] +smaller_Vecs = randomSublist 0.25 small_Vecs + +small_pairs :: [(Vec,Vec)] +small_pairs = [ (u,v) | u <- smaller_Vecs , v <- smaller_Vecs ] + -------------------------------------------------------------------------------- +randomVec :: Int -> Int -> IO Vec +randomVec maxlen maxbits = do + len <- randomRIO (0,maxlen ) + bits <- randomRIO (0,maxbits) + let bnd = 2^bits - 1 + xs <- replicateM len (randomRIO (0,bnd)) + -- print (bnd,xs) + return $ Vec $ V.fromListN len bnd xs + +{-# NOINLINE rnd_Vecs #-} +rnd_Vecs :: [Vec] +rnd_Vecs = L.concat $ L.concat $ Unsafe.unsafePerformIO $ do + forM [1,5..100] $ \len -> do + forM [4,8..64] $ \maxbits -> do + let k = maxbits -- for smaller bit depths, we need less test cases + replicateM k (randomVec len maxbits) + +rnd_Lists = [ List (V.toList vec) | Vec vec <- rnd_Vecs ] + +rnd_NELists = [ NEList l | List l <- rnd_Lists , not (L.null l) ] +rnd_NEVecs = [ NEVec v | Vec v <- rnd_Vecs , not (V.null v) ] + +rnd_cons = [ (x, Vec (V.fromList xs)) | NEList (x:xs) <- rnd_NELists ] + +smaller_rnd_Vecs :: [Vec] +smaller_rnd_Vecs = randomSublist 0.025 rnd_Vecs + +rnd_pairs :: [(Vec,Vec)] +rnd_pairs = [ (u,v) | u <- smaller_rnd_Vecs , v <- smaller_rnd_Vecs ] + +-------------------------------------------------------------------------------- + add_bighead :: List -> [List] add_bighead (List xs) = [ List (2^k-1 : xs) | k<-[1..arch_bits-1] ] ++ [ List (2^k : xs) | k<-[1..arch_bits-1] ] ++ [ List (2^k+1 : xs) | k<-[1..arch_bits-1] ] -bighead_Lists = concatMap add_bighead small_Lists :: [List] -bighead_Vecs = [ Vec (V.fromList xs) | List xs <- bighead_Lists ] :: [Vec] +bighead_Lists_orig = concatMap add_bighead small_Lists :: [List] +bighead_Lists = randomSublist 0.25 bighead_Lists_orig :: [List] + +bighead_Vecs = [ Vec (V.fromList xs) | List xs <- bighead_Lists ] :: [Vec] bighead_NELists = [ NEList xs | List xs <- bighead_Lists ] :: [NEList] bighead_NEVecs = [ NEVec v | Vec v <- bighead_Vecs ] :: [NEVec] +bighead_cons = [ (x, Vec (V.fromList xs)) | NEList (x:xs) <- bighead_NELists ] + +smaller_bighead_Vecs :: [Vec] +smaller_bighead_Vecs = randomSublist 0.01 bighead_Vecs + +bighead_pairs :: [(Vec,Vec)] +bighead_pairs = [ (u,v) | u <- smaller_bighead_Vecs , v <- smaller_bighead_Vecs ] + -------------------------------------------------------------------------------- +-- * misc tests + +-- cons fatal error (caused by non-strict left shift of blobs): +cons_v2_crash = V.toList (cons_v2 x xs) == [159407,244557,137175,96511,42979] where + (x, Vec xs) = (159407, Vec (fromList' (Shape {shapeLen = 4, shapeBits = 20}) [244557,137175,96511,42979])) + +check_eq_bitsizes = and + [ V.fromList [k] == V.fromList' (Shape 1 b) [k] + | k<-[0..15] + , b<-[4..31] + ] + +around_powers_of_two :: [Word] +around_powers_of_two = [0..7] ++ stuff ++ [2^nn-i | i<-[1..8] ] where + stuff = [ 2^i + fromIntegral ofs | i<-[2..nn-1] , ofs <- [-2..2::Int] ] + nn = {- 64 -} arch_bits + +-------------------------------------------------------------------------------- -- * properties prop_from_to_list (List list) = V.toList (V.fromList list) == list prop_to_from_vec (Vec vec ) = V.fromList (V.toList vec ) == vec -prop_tolist_vs_naive (Vec vec) = (V.toList vec == V.toList_naive vec) +prop_tolist_vs_naive (Vec vec) = (V.toList vec == toList_extract vec) +prop_toRevList (Vec vec) = V.toRevList vec == reverse (V.toList vec) prop_fromlist_vs_index (List list) = [ unsafeIndex i vec | i<-[0..n-1] ] == list where vec = V.fromList list n = V.vecLen vec prop_head_of_list (NEList list) = V.head (V.fromList list) == L.head list +prop_last_of_list (NEList list) = V.last (V.fromList list) == L.last list +prop_tail_of_list (NEVec vec ) = V.toList (V.tail vec) == L.tail (V.toList vec) prop_head_vs_index (NEVec vec ) = V.head vec == unsafeIndex 0 vec + +prop_cons_uncons (NEVec vec) = liftM (uncurry V.cons) (V.uncons vec) == Just vec +prop_uncons_cons (w,Vec vec) = V.uncons (V.cons w vec) == Just (w,vec) +prop_uncons_vs_naive (Vec vec) = V.uncons vec == uncons_naive vec + +prop_uncons_vs_list (Vec vec) = unconsToList (V.uncons vec) == L.uncons (V.toList vec) +prop_cons_vs_list (w,Vec vec) = V.toList (V.cons w vec) == w : (V.toList vec) +prop_snoc_vs_list (w,Vec vec) = V.toList (V.snoc vec w) == (V.toList vec) ++ [w] + +unconsToList mb = case mb of + Nothing -> Nothing + Just (w,vec) -> Just (w, V.toList vec) + +-------------------------------------------------------------------------------- + +prop_tail_v1_vs_v2 (Vec vec) = tail_v1 vec == tail_v2 vec +prop_cons_v1_vs_v2 (y,(Vec vec)) = cons_v1 y vec == cons_v2 y vec +prop_snoc_v1_vs_v2 (y,(Vec vec)) = snoc_v1 vec y == snoc_v2 vec y +prop_uncons_v1_vs_v2 (Vec vec ) = uncons_v1 vec == uncons_v2 vec + +{- +bad_tail = [ vec | Vec vec <- bighead_Vecs , tail_v1 vec /= tail_v2 vec ] +bad_cons = [ (y,vec) | (y,Vec vec) <- bighead_cons , cons_v1 y vec /= cons_v2 y vec ] + +bad_cons_uncons = [ vec | NEVec vec <- bighead_NEVecs , liftM (uncurry V.cons) (V.uncons vec) /= Just vec ] +-} + +bad_snoc = [ (y,vec) | (y,Vec vec) <- small_cons , snoc_v1 vec y /= snoc_v2 vec y ] + +-------------------------------------------------------------------------------- + +prop_max_vs_list (Vec vec) = V.maximum vec == listMax (V.toList vec) +prop_sum_vs_list (Vec vec) = V.sum vec == L.sum (V.toList vec) + +{- +bad_max = [ vec | v@(Vec vec) <- bighead_Vecs , not (prop_max_vs_list v) ] +bad_sum = [ vec | v@(Vec vec) <- bighead_Vecs , not (prop_sum_vs_list v) ] +-} + +-------------------------------------------------------------------------------- +-- TODO: randomized tests for these + +eqListExt0 :: [Word] -> [Word] -> Bool +eqListExt0 xs ys = and (listLongZipWith (==) xs ys) + +leListExt0 :: [Word] -> [Word] -> Bool +leListExt0 xs ys = and (listLongZipWith (<=) xs ys) + +prop_strict_eq_vs_list (Vec vec1 , Vec vec2) = eqStrict vec1 vec2 == (V.toList vec1 == V.toList vec2) +prop_strict_cmp_vs_list (Vec vec1 , Vec vec2) = cmpStrict vec1 vec2 == (vec1 `cmpStrict_naive` vec2) +prop_ext0_eq_vs_list (Vec vec1 , Vec vec2) = eqExtZero vec1 vec2 == eqListExt0 (V.toList vec1) (V.toList vec2) +prop_ext0_cmp_vs_list (Vec vec1 , Vec vec2) = cmpExtZero vec1 vec2 == (vec1 `cmpExtZero_naive` vec2) +prop_less_or_equal_vs_list (Vec vec1 , Vec vec2) = lessOrEqual vec1 vec2 == leListExt0 (V.toList vec1) (V.toList vec2) + +-------------------------------------------------------------------------------- + +prop_add_commutative (Vec vec1 , Vec vec2) = V.add vec1 vec2 == V.add vec2 vec1 +prop_add_vs_naive (Vec vec1 , Vec vec2) = V.add vec1 vec2 == add_naive vec1 vec2 +prop_sub_vs_naive (Vec vec1 , Vec vec2) = V.subtract vec1 vec2 == sub_naive vec1 vec2 + +-- bad_add = [ (bad1,bad2) | (Vec bad1, Vec bad2) <- bighead_pairs , not (prop_add_vs_naive (Vec bad1 , Vec bad2)) ] + +prop_psums_vs_naive (Vec vec) = V.partialSums vec == partialSums_naive vec +prop_scale_vs_naive s (Vec vec) = V.scale s vec == scale_naive s vec + +bad_psums = [ bad | Vec bad <- bighead_Vecs , not (prop_psums_vs_naive (Vec bad)) ] +bad_scale3 = [ bad | Vec bad <- bighead_Vecs , not (prop_scale_vs_naive 3 (Vec bad)) ] --------------------------------------------------------------------------------