build_bdr() { char comm[100], line[100], *strsave(), *ckalloc(); int i, j; for (i = 0; i < n_seq-1; i++) { sprintf(line, "temp%d", i); sequence[i].sim = strsave(line); j = sequence[n_seq-1].hi_addr - sequence[n_seq-1].lo_addr; upper_bound[i] = (int *) ckalloc(sizeof(int)*j); lower_bound[i] = (int *) ckalloc(sizeof(int)*j); sprintf(comm, "%s 1 %s %s>%s", sim, sequence[n_seq-1].name, sequence[i].name, line); if (!sim_ready) system(comm); Read_SIM(sequence[i].sim, width, upper_bound[i], lower_bound[i]); } } int change_p(x, i, a) int x, i; int a; { int t, y,j; j = abs(i); t = ALLONE - mask[i]*3; y = x & t; y += a * mask[i]; return y; } int find_pos(x, i) int x, i; { int y; y = (x / mask[i]) & 3; return y; } next_p(current, avoid) int current, avoid; { int a; a = current+1; if (a != avoid) { if (a < 4) return a; if (avoid != 0) return -2; return -1; } a++; if (a <= 3) return a; return -2; } /* build blocks*/ int a[MAXSEQ], bk[MAXSEQ], c[MAXSEQ], n[MAXSEQ], m[MAXSEQ], bkt, bbkt; Hash_ptr aq[MAXSEQ]; /* init(k, i) give the first candidate of sequence k that matches position i of sequence n_seq-1, the base sequence. next_q(k, i) give the consequence of the init. There are two position of interests, indicated by lop and loq. lop indicate the first sequence that is different by two positions from the base sequence. loq indicate the position that we can be clear of the basis. */ int mm, aa, bb1, aa2, bb2, bktbkt, lop = MAXSEQ, loq = MAXSEQ; init(k, i) int k, i; { int j, v, t, a1, m, b1; /*if (i == 76) printf("init(%d, %d)", k, i);*/ if (k <= lop) lop = MAXSEQ; if (k >= n_seq-1) return 1; for (j = k; j < n_seq - 1; j++) { if (j > lop) { first_q(j); } else { a[j] = 0; bk[j] = bkt; n[j] = 0; c[j] = 0; v = find_pos(bkt, a[j]); if (v == 0) c[j] = 1; } while (a[j] < hash_size) { aq[j] = get_head(bk[j], j); while (aq[j] && aq[j]->pos <= low[j]) aq[j] = aq[j]->next; if (aq[j] == NULL || aq[j]->pos >= up[j]) { if (j > lop) if (!next_qq(j)) return next_q(j-1, i); else continue; bk[j] = change_p(bkt, a[j], c[j]); if (n[j]) t = next_p(3, v); else t = next_p(c[j], v); if (t < 0) { t += 2; a[j]++; c[j] = t; n[j] = 0; for (t = 0; t < j; t++) if (a[t] == a[j]) { c[j] = c[t]; n[j] = 1; break; } } else c[j] = t; } else { break;} } if (a[j] >= hash_size) { if (k > lop) return next_q(j-1, i); for (m = 0, t = 0; t < j; t++) { if (a[t] == 0) continue; if (!m) {a1 = a[t]; b1 = c[t]; m = 1; continue;} if (a1 != a[t]) return next_q(j-1, i); } for (bk[j] = first_p(m, a1, b1, bkt); bk[j] >= 0;) { aq[j] = get_head(bk[j], j); while (aq[j] && aq[j]->pos <= low[j]) aq[j] = aq[j]->next; if (aq[j] == NULL || aq[j]->pos >= up[j]) bk[j] = next_pp(); else break; } if (bk[j] < 0) return next_q(j-1, i); lop = j; loq = MAXSEQ; } } return 1; } next_q(k, i) int k, i; { int j, v, t; /*if (i == 76) printf("next_q(%d, %d)\n", k, i);*/ if (k < 0) return 0; j = k; for (;!(aq[j] = aq[j]->next) || (aq[j]->pos >= up[j]);) { while (a[j] < hash_size) { if (j >= lop) if (j == lop) { for (bk[j] = next_pp(); bk[j] >= 0; bk[j] = next_pp()) { aq[j] = get_head(bk[j], j); while (aq[j] && aq[j]->pos <= low[j]) aq[j] = aq[j]->next; if (!aq[j] || aq[j]->pos >= up[j]) continue; else break; } /*for*/ if (bk[j] < 0) break; /*while*/ } else { if (!next_qq(j)) return next_q(j-1, i); aq[j] = get_head(bk[j], j); } else { bk[j] = change_p(bkt, a[j], c[j]); aq[j] = get_head(bk[j], j); v = find_pos(bkt, a[j]); if (n[j]) t = next_p(3, v); else t = next_p(c[j], v); if (t < 0) { t += 2; a[j]++; c[j] = t; n[j] = 0; for (t = 0; t < j; t++) if (a[t] == a[j]) { c[j] = c[t]; n[j] = 1; break; } } else c[j] = t; } while (aq[j] && aq[j]->pos <= low[j]) aq[j] = aq[j]->next; if (aq[j] && aq[j]->pos < up[j]) return init(j+1, i); } if (j == 0) return 0; j--; } /*for*/ return init(j+1, i); } int first_p(m, a1, b1, bkt) int m, a1, b1,bkt; { int v1, v2; if (m == 1) { mm = 1; aa = a1; bb1 = b1; bbkt = bktbkt = change_p(bkt, aa, bb1); if (aa != 0) aa2 = 0; else aa2 = 1; v1 = find_pos(bkt, aa2); if (v1 == 0) bb2 = 1; else bb2 = 0; return(change_p(bktbkt, aa2, bb2)); } /*m == 0*/ mm = 0; bktbkt = -1; aa = 0; aa2 = 1; v1 = find_pos(bkt, aa); v2 = find_pos(bkt, aa2); if (!v1) bb1 = 1; else bb1 = 0; if (v2) bb1 = 0; else bb1 = 1; return change_p(change_p(bkt, aa, bb1), aa2, bb2); } int next_pp() { int v1, v2; if (mm == 1) { bb2 = next_p(bb2, v1); if (bb2 < 0) { bb2 += 2; aa2++; if (aa2 == aa) aa2++; if (aa2 >= hash_size) return -1; } return change_p(bktbkt, aa2, bb2); } bb2 = next_p(bb2, v2); if (bb2 < 0) { bb2 += 2; aa2++; if (aa2 >= hash_size) { bb1 = next_p(bb1, v1); if (bb1 < 0) { bb1+=2; aa++; if (aa >= hash_size-1) return -1; } aa2 = aa+1; bb2 = 0; } } return change_p(change_p(bkt, aa, bb1), aa2, bb2); } void check_c(i1, i2, j) int i1, i2, j; { int i; n[j] = 0; for (i = i1; i < i2; i++) if (a[i] == a[j]) { c[j] = c[i]; n[j] = 1; break; } } /* first_q give the first candidate of sequence j if j is larger lop. next_qq gives later ones. */ first_q(j) int j; { int v, v1, v2; if (mm == 1) { bk[j] = bktbkt; a[j] = 0; v = find_pos(bktbkt, 0); if (v) c[j] = 0; else c[j] = 1; check_c(lop+1, j, j); return; } if (j > loq) { bk[j] = bktbkt; a[j] = 0; n[j] = 0; if (0 == aa) { if ((v1=find_pos(bkt, 0)) != (v2=find_pos(bktbkt, 0))) { c[j] = v1; n[j] = 1; return; } } check_c(loq+1, j, j); return; } bk[j] = bk[lop]; a[j] = 0; v = find_pos(bktbkt, 0); if (v) c[j] = 0; else c[j] = 1; } int fs; next_qq(j) int j; { int v, v1, v2, t; if (j == lop) fatal("P"); if (mm == 1) { if (a[j] >= hash_size) return 0; bk[j] = change_p(bktbkt, a[j], c[j]); if (!n[j]) { v = find_pos(bktbkt, a[j]); v = next_p(c[j], v); if (v >= 0) { c[j] = v; return 1; } } /* all other case need to add 1 to a[j]*/ a[j]++; if (a[j] >= hash_size) return 1; v = find_pos(bktbkt, a[j]); if (v) c[j] = 0; else c[j] = 1; check_c(lop+1, j, j); return 1; } if (j >= loq) { if (a[j] >= hash_size) { if (j == loq && fs == 1) { bbkt = bktbkt = bk[j] = change_p(bkt, aa2, bb2); a[j] = 0; v = find_pos(bktbkt, 0); if (v) c[j] = 0; else c[j] = 1; fs++; return 1; } else return 0; } bk[j] = change_p(bktbkt, a[j], c[j]); if (!n[j]) { v = find_pos(bktbkt, a[j]); v = next_p(c[j], v); if (v >= 0) { c[j] = v; return 1; } } /* otherwise all need to add 1 to a[j];*/ a[j]++; if (j == loq) { if (a[j] == aa || a[j] == aa2) { a[j]++; if (a[j] == aa2) a[j]++; } } if (a[j] >= hash_size) return 1; if (j > loq) { if (a[j] == aa || a[j] == aa2) { if ((v1 = find_pos(bkt, 0)) == (v2=find_pos(bktbkt, 0))) if (a[j] == aa) c[j] = bb1; else c[j] = bb2; else {c[j] = v1; n[j] = 1;} return 1; } } v = find_pos(bktbkt, a[j]); if (v) c[j] = 0; else c[j] = 1; check_c(loq+1, j, j); return 1; } if (j < loq) { if (bk[j] == bk[lop]) bk[j] = bkt; else { loq = j; fs = 1; bktbkt = bk[j] = change_p(bkt, aa, bb1); a[j] = 0; v = find_pos(bktbkt, 0); if (v) c[j] = 0; else c[j] = 1; } return 1; } }