shingles.c

   1 /*
   2   this is an attempted implementation of Broders shingles algorithm for similarity testing of documents
   3
   4   see http://gatekeeper.dec.com/pub/DEC/SRC/publications/broder/positano-final-wpnums.pdf
   5
   6   tridge@samba.org
   7 */
   8 #include <stdio.h>
   9 #include <stdlib.h>
  10 #include <string.h>
  11 #include <fcntl.h>
  12 #include <errno.h>
  13 #include <sys/mman.h>
  14 #include <sys/stat.h>
  15 #include <unistd.h>
  16 #include <ctype.h>
  17 #include <time.h>
  18
  19 #define MAX_WORD_SIZE 8
  20 #define W_SHINGLES 8
  21 #define SKETCH_SIZE 32
  22
  23 #define HASH_PRIME 0x01000193
  24 #define HASH_INIT 0x28021967
  25
  26 #define FLAG_IGNORE_HEADERS 1
  27
  28 #define ROLLING_WINDOW 20
  29
  30 typedef unsigned u32;
  31 typedef unsigned char uchar;
  32
  33 #define QSORT_CAST (int (*)(const void *, const void *))
  34
  35
  36 struct u32_set {
  37         int count;
  38         u32 v[SKETCH_SIZE*2];
  39 };
  40
  41
  42 static void parse_sketch(const uchar *str, struct u32_set *v)
  43 {
  44         int i, j;
  45         v->count = 0;
  46         for (i=0;i<SKETCH_SIZE;i++) {
  47                 v->v[i] = 0;
  48                 if (! *str) break;
  49                 for (j=0;j<4;j++) {
  50                         v->v[i] = (v->v[i] << 7) | *str++;
  51                 }
  52                 v->count++;
  53         }
  54 }
  55
  56 static int u32_comp(u32 *v1, u32 *v2)
  57 {
  58         if (*v1 > *v2) {
  59                 return 1;
  60         }
  61         if (*v1 < *v2) {
  62                 return -1;
  63         }
  64         return 0;
  65 }
  66
  67 static void u32_union(struct u32_set *v1, struct u32_set *v2, struct u32_set *u)
  68 {
  69         int i, n;
  70         u32 vv[SKETCH_SIZE*2];
  71
  72         n = v1->count + v2->count;
  73         memcpy(vv, v1->v, v1->count * sizeof(u32));
  74         memcpy(&vv[v1->count], v2->v, v2->count * sizeof(u32));
  75
  76         qsort(vv, n, sizeof(u32), QSORT_CAST u32_comp);
  77
  78         u->count=1;
  79         u->v[0] = vv[0];
  80
  81         for (i=1;i<n;i++) {
  82                 if (vv[i] != u->v[u->count-1]) {
  83                         u->v[u->count++] = vv[i];
  84                 }
  85         }
  86 }
  87
  88 static void u32_intersect(struct u32_set *v1, struct u32_set *v2, struct u32_set *u)
  89 {
  90         u32 vv[SKETCH_SIZE*2];
  91         int i;
  92
  93         memcpy(vv, v1->v, v1->count*sizeof(u32));
  94         memcpy(&vv[v1->count], v2->v, v2->count*sizeof(u32));
  95         qsort(vv, v1->count + v2->count, sizeof(u32), QSORT_CAST u32_comp);
  96
  97         u->count=0;
  98
  99         for (i=1;i<v1->count + v2->count; i++) {
 100                 if (vv[i] == vv[i-1]) {
 101                         u->v[u->count++] = vv[i];
 102                 }
 103         }
 104 }
 105
 106 static void u32_min(struct u32_set *v)
 107 {
 108         qsort(v->v, v->count, sizeof(u32), QSORT_CAST u32_comp);
 109         if (v->count > SKETCH_SIZE) {
 110                 v->count = SKETCH_SIZE;
 111         }
 112 }
 113
 114 /*
 115   given two sketch strings return a value indicating the degree to which they match.
 116 */
 117 static u32 sketch_match(const char *str1, const char *str2)
 118 {
 119         u32 score = 0;
 120         struct u32_set v1, v2;
 121         struct u32_set u1;
 122         struct u32_set v3, v4;
 123
 124         parse_sketch(str1, &v1);
 125         parse_sketch(str2, &v2);
 126
 127         u32_union(&v1, &v2, &u1);
 128
 129         u32_min(&u1);
 130         u32_intersect(&u1, &v1, &v3);
 131         u32_intersect(&v3, &v2, &v4);
 132
 133         score = (100 * v4.count) / u1.count;
 134
 135         return score;
 136 }
 137
 138 /*
 139   return the maximum match for a file containing a list of sketchs
 140 */
 141 static u32 sketch_match_db(const char *fname, const char *sum, u32 threshold)
 142 {
 143         FILE *f;
 144         char line[200];
 145         u32 best = 0;
 146
 147         f = fopen(fname, "r");
 148         if (!f) return 0;
 149
 150         /* on each line of the database we compute the sketch match
 151            score. We then pick the best score */
 152         while (fgets(line, sizeof(line)-1, f)) {
 153                 u32 score;
 154                 int len;
 155                 len = strlen(line);
 156                 if (line[len-1] == '\n') line[len-1] = 0;
 157
 158                 score = sketch_match(sum, line);
 159
 160                 if (score > best) {
 161                         best = score;
 162                         if (best >= threshold) break;
 163                 }
 164         }
 165
 166         fclose(f);
 167
 168         return best;
 169 }
 170
 171
 172 /* a simple non-rolling hash, based on the FNV hash */
 173 static inline u32 sum_hash(uchar c, u32 h)
 174 {
 175         h *= HASH_PRIME;
 176         h ^= c;
 177         return h;
 178 }
 179
 180 static struct {
 181         uchar window[ROLLING_WINDOW];
 182         u32 h1, h2, h3;
 183         u32 n;
 184 } roll_state;
 185
 186 /*
 187   a rolling hash, based on the Adler checksum. By using a rolling hash
 188   we can perform auto resynchronisation after inserts/deletes
 189
 190   internally, h1 is the sum of the bytes in the window and h2
 191   is the sum of the bytes times the index
 192
 193   h3 is a shift/xor based rolling hash, and is mostly needed to ensure that
 194   we can cope with large blocksize values
 195 */
 196 static inline u32 roll_hash(uchar c)
 197 {
 198         roll_state.h2 -= roll_state.h1;
 199         roll_state.h2 += ROLLING_WINDOW * c;
 200
 201         roll_state.h1 += c;
 202         roll_state.h1 -= roll_state.window[roll_state.n % ROLLING_WINDOW];
 203
 204         roll_state.window[roll_state.n % ROLLING_WINDOW] = c;
 205         roll_state.n++;
 206
 207         roll_state.h3 = (roll_state.h3 << 5) & 0xFFFFFFFF;
 208         roll_state.h3 ^= c;
 209
 210         return roll_state.h1 + roll_state.h2 + roll_state.h3;
 211 }
 212
 213 /*
 214   reset the state of the rolling hash and return the initial rolling hash value
 215 */
 216 static u32 roll_reset(void)
 217 {
 218         memset(&roll_state, 0, sizeof(roll_state));
 219         return 0;
 220 }
 221
 222 static uchar *advance(uchar *in, unsigned w, uchar *limit)
 223 {
 224         in++;
 225         while (w--) {
 226                 int i;
 227                 for (i=0;i<MAX_WORD_SIZE;i++) {
 228                         roll_hash(in[i]);
 229                         if (in >= limit || !isalnum(in[i])) break;
 230                 }
 231                 in += i;
 232                 while (in < limit && !isalnum(in[0])) {
 233                         in++;
 234                         roll_hash(in[i]);
 235                 }
 236         }
 237         return in;
 238 }
 239
 240 static u32 hash_emit(uchar *start, size_t len)
 241 {
 242         u32 v = HASH_INIT;
 243
 244         if (len == 0) return 0;
 245
 246         while (len--) {
 247                 if (isalnum(*start)) {
 248                         v = sum_hash(*start, v);
 249                 }
 250                 start++;
 251         }
 252         return v;
 253 }
 254
 255 struct el {
 256         u32 hash;
 257         int i;
 258 };
 259
 260 static int el_comp(struct el *v1, struct el *v2)
 261 {
 262         if (v2->hash > v1->hash) {
 263                 return 1;
 264         }
 265         if (v2->hash < v1->hash) {
 266                 return -1;
 267         }
 268         return 0;
 269 }
 270
 271 static int el_comp2(struct el *v1, struct el *v2)
 272 {
 273         return v1->i - v2->i;
 274 }
 275
 276 static uchar *shingle_sketch(uchar *in_0, size_t size, u32 flags)
 277 {
 278         uchar *in = in_0;
 279         uchar *p;
 280         struct el *sketch, *sketch2;
 281         uchar *ret;
 282         u32 sketch_size = 0;
 283         int i, n;
 284         const char *b64 = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/";
 285
 286         if (size == 0) {
 287                 return strdup("NULL");
 288         }
 289
 290         /* if we are ignoring email headers then skip past them now */
 291         if (flags & FLAG_IGNORE_HEADERS) {
 292                 uchar *s = strstr(in, "\n\n");
 293                 if (s) {
 294                         size -= (s+2 - in);
 295                         in = in_0 = s+2;
 296                 }
 297         }
 298
 299         /* allocate worst case size */
 300         sketch  = calloc(sizeof(struct el), SKETCH_SIZE + size+1);
 301         sketch2 = calloc(sizeof(struct el), SKETCH_SIZE + size+1);
 302         ret = calloc(SKETCH_SIZE+1, 4);
 303
 304         roll_reset();
 305
 306         p = advance(in, W_SHINGLES, in_0 + size);
 307
 308         sketch[sketch_size].hash = hash_emit(in, p-in);
 309         sketch[sketch_size].i = sketch_size++;
 310
 311         while (p < in_0 + size) {
 312                 in = advance(in, 1, in_0 + size);
 313                 p = advance(p, 1, in_0 + size);
 314                 sketch[sketch_size].hash = hash_emit(in, p-in);
 315                 sketch[sketch_size].i = sketch_size++;
 316         }
 317
 318         qsort(sketch, sketch_size, sizeof(struct el), QSORT_CAST el_comp);
 319
 320         sketch2[0] = sketch[0];
 321         n = 1;
 322         for (i=1;i<sketch_size;i++) {
 323                 if (sketch[i].hash != sketch2[n-1].hash) {
 324                         sketch2[n++] = sketch[i];
 325                 }
 326         }
 327         free(sketch);
 328         sketch = sketch2;
 329         sketch_size = n;
 330
 331         if (sketch_size > SKETCH_SIZE) {
 332                 qsort(sketch, SKETCH_SIZE, sizeof(struct el), QSORT_CAST el_comp2);
 333         }
 334
 335         for (i=0;i<SKETCH_SIZE && i < sketch_size;i++) {
 336                 u32 v = sketch[i].hash;
 337                 int j;
 338                 for (j=0;j<4;j++) {
 339                         ret[4*i + j] = b64[v % 64];
 340                         v >>= 6;
 341                 }
 342         }
 343         ret[4*i] = 0;
 344
 345         free(sketch);
 346
 347         return ret;
 348 }
 349
 350 /*
 351   return the shingles sketch on stdin
 352 */
 353 static char *sketch_stdin(u32 flags)
 354 {
 355         uchar buf[10*1024];
 356         uchar *msg;
 357         size_t length = 0;
 358         int n;
 359         char *sum;
 360
 361         msg = malloc(sizeof(buf));
 362         if (!msg) return NULL;
 363
 364         /* load the file, expanding the allocation as needed. */
 365         while (1) {
 366                 n = read(0, buf, sizeof(buf));
 367                 if (n == -1 && errno == EINTR) continue;
 368                 if (n <= 0) break;
 369
 370                 msg = realloc(msg, length + n);
 371                 if (!msg) return NULL;
 372
 373                 memcpy(msg+length, buf, n);
 374                 length += n;
 375         }
 376
 377         sum = shingle_sketch(msg, length, flags);
 378
 379         free(msg);
 380
 381         return sum;
 382 }
 383
 384
 385 /*
 386   return the sketch on a file
 387 */
 388 char *sketch_file(const char *fname, u32 flags)
 389 {
 390         int fd;
 391         char *sum;
 392         struct stat st;
 393         uchar *msg;
 394
 395         if (strcmp(fname, "-") == 0) {
 396                 return sketch_stdin(flags);
 397         }
 398
 399         fd = open(fname, O_RDONLY);
 400         if (fd == -1) {
 401                 perror(fname);
 402                 return NULL;
 403         }
 404
 405         if (fstat(fd, &st) == -1) {
 406                 perror("fstat");
 407                 return NULL;
 408         }
 409
 410         msg = mmap(NULL, st.st_size, PROT_READ, MAP_FILE|MAP_PRIVATE, fd, 0);
 411         if (msg == (uchar *)-1) {
 412                 perror("mmap");
 413                 return NULL;
 414         }
 415         close(fd);
 416
 417         sum = shingle_sketch(msg, st.st_size, flags);
 418
 419         munmap(msg, st.st_size);
 420
 421         return sum;
 422 }
 423
 424
 425 static void show_help(void)
 426 {
 427         printf("
 428 shingles v1.0 written by Andrew Tridgell <tridge@samba.org>
 429
 430 Syntax:
 431    shingles [options] <files>
 432 or
 433    shingles [options] -d sigs.txt -c SIG
 434 or
 435    shingles [options] -d sigs.txt -C file
 436
 437 When called with a list of filenames shingles will write out the
 438 signatures of each file on a separate line. You can specify the
 439 filename '-' for standard input.
 440
 441 When called with the second form, shingles will print the best score
 442 for the given signature with the signatures in the given database. A
 443 score of 100 means a perfect match, and a score of 0 means a complete
 444 mismatch.
 445
 446 When checking, shingles returns 0 (success) when the message *is* spam,
 447 1 for internal errors, and 2 for messages whose signature is not
 448 found.
 449
 450 The 3rd form is just like the second form, but you pass a file
 451 containing a message instead of a pre-computed signature.
 452
 453 Options:
 454    -W              ignore whitespace
 455    -H              skip past mail headers
 456    -T <threshold>  set the threshold above which shingles will stop
 457                    looking (default 90)
 458 ");
 459 }
 460
 461
 462 int main(int argc, char *argv[])
 463 {
 464         char *sum;
 465         extern char *optarg;
 466         extern int optind;
 467         int c;
 468         char *dbname = NULL;
 469         u32 score;
 470         int i;
 471         u32 flags = 0;
 472         u32 threshold = 90;
 473
 474         while ((c = getopt(argc, argv, "B:WHd:c:C:hT:")) != -1) {
 475                 switch (c) {
 476                 case 'H':
 477                         flags |= FLAG_IGNORE_HEADERS;
 478                         break;
 479
 480                 case 'd':
 481                         dbname = optarg;
 482                         break;
 483
 484                 case 'T':
 485                         threshold = atoi(optarg);
 486                         break;
 487
 488                 case 'c':
 489                         if (!dbname) {
 490                                 show_help();
 491                                 exit(1);
 492                         }
 493                         score = sketch_match_db(dbname, optarg,
 494                                                  threshold);
 495                         printf("%u\n", score);
 496                         exit(score >= threshold ? 0 : 2);
 497
 498                 case 'C':
 499                         if (!dbname) {
 500                                 show_help();
 501                                 exit(1);
 502                         }
 503                         score = sketch_match_db(dbname,
 504                                                  sketch_file(optarg, flags),
 505                                                  threshold);
 506                         printf("%u\n", score);
 507                         exit(score >= threshold ? 0 : 2);
 508
 509                 case 'h':
 510                 default:
 511                         show_help();
 512                         exit(0);
 513                 }
 514         }
 515
 516         argc -= optind;
 517         argv += optind;
 518
 519         if (argc == 0) {
 520                 show_help();
 521                 return 0;
 522         }
 523
 524         /* compute the sketch on a list of files */
 525         for (i=0;i<argc;i++) {
 526                 sum = sketch_file(argv[i], flags);
 527                 printf("%s\n", sum);
 528                 free(sum);
 529         }
 530
 531         return 0;
 532 }