+ target = 0;
+ for (i1 = 1; i1 < *n; i1++) {
+ for (i2 = i1 + 1; i2 <= *n; i2++) {
+ N = 0;
+
+ for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n)
+ if ((x[s1] ^ x[s2]) & 4) N++;
+
+ d[target] = ((double) N);
+ target++;
+ }
+ }
+}
+
+#define X(i, j) i - 1 + *n * (j - 1)
+
+#define DINDEX(i, j) n * (i - 1) - i*(i - 1)/2 + j - i - 1
+
+int give_index(int i, int j, int n)
+{
+ if (i > j) return(DINDEX(j, i));
+ else return(DINDEX(i, j));
+}
+
+void distDNA_indelblock(unsigned char *x, int *n, int *s, double *d)
+{
+ int i1, i2, s1, s2, target, N, start_block, end_block;
+
+ for (i1 = 1; i1 <= *n; i1++) {
+
+/* find a block of one or more '-' */
+
+ s1 = 1;
+ while (s1 < *s) {
+ if (x[X(i1, s1)] == 4) {
+ start_block = s1;
+ while (x[X(i1, s1)] == 4) s1++;
+ end_block = s1 - 1;
+
+/* then look if the same block is present in all the other sequences */
+
+ for (i2 = 1; i2 <= *n; i2++) {
+ if (i2 == i1) continue;
+
+ target = give_index(i1, i2, *n);
+
+ if (start_block > 1) {
+ if (x[X(i2, start_block - 1)] == 4) {
+ d[target]++;
+ continue;
+ }
+ }
+ if (end_block < *s) {
+ if (x[X(i2, end_block + 1)] == 4) {
+ d[target]++;
+ continue;
+ }
+ }
+ for (s2 = start_block; s2 <= end_block; s2++) {
+ if (x[X(i2, s2)] != 4) {
+ d[target]++;
+ continue;
+ }
+ }
+ }
+ s1 = end_block + 1;
+ }
+ s1++;
+ }
+ }
+}
+
+#undef X
+
+void distDNA_TsTv(unsigned char *x, int *n, int *s, double *d, int Ts, int pairdel)
+{
+ int i1, i2, s1, s2, target, Nd, Ns;
+
+ target = 0;
+ for (i1 = 1; i1 < *n; i1++) {
+ for (i2 = i1 + 1; i2 <= *n; i2++) {
+ Nd = Ns = 0;
+ for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
+ if (pairdel && !(KnownBase(x[s1]) && KnownBase(x[s2]))) continue;
+ COUNT_TS_TV
+ }
+ if (Ts) d[target] = ((double) Ns); /* output number of transitions */
+ else d[target] = ((double) Nd - Ns); /* output number of transversions */
+ target++;
+ }
+ }
+}
+
+void distDNA_raw(unsigned char *x, int *n, int *s, double *d, int scaled)
+{
+ int i1, i2, s1, s2, target, Nd;
+
+ target = 0;
+ for (i1 = 1; i1 < *n; i1++) {
+ for (i2 = i1 + 1; i2 <= *n; i2++) {
+ Nd = 0;
+ for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n)
+ if (DifferentBase(x[s1], x[s2])) Nd++;
+ if (scaled) d[target] = ((double) Nd / *s);
+ else d[target] = ((double) Nd);
+ target++;
+ }