X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=src%2Fdist_dna.c;h=7174af6d013d961c80dba33b453c78ad5e3c8bfe;hb=1ad48c7a70983375138a6500372db588c8a3a134;hp=62f8f64c9b5b375a702063f7d11825527a5d6d8b;hpb=bfaeca35ec130110327a3bb7a1f0fe3b66076a95;p=ape.git

diff --git a/src/dist_dna.c b/src/dist_dna.c
index 62f8f64..7174af6 100644
--- a/src/dist_dna.c
+++ b/src/dist_dna.c
@@ -1,6 +1,6 @@
-/* dist_dna.c       2011-02-18 */
+/* dist_dna.c       2012-01-10 */
 
-/* Copyright 2005-2011 Emmanuel Paradis
+/* Copyright 2005-2012 Emmanuel Paradis
 
 /* This file is part of the R-package `ape'. */
 /* See the file ../COPYING for licensing issues. */
@@ -535,8 +535,8 @@ void distDNA_TN93(unsigned char *x, int *n, int *s, double *d,
 		  double *BF, int *variance, double *var,
 		  int *gamma, double *alpha)
 {
-    int i1, i2, k, Nd, Ns1, Ns2, L, target, s1, s2;
-    double P1, P2, Q, A, B, C, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b;
+    int i1, i2, Nd, Ns1, Ns2, L, target, s1, s2;
+    double P1, P2, Q, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b;
 
     L = *s;
 
@@ -559,8 +559,8 @@ void distDNA_TN93_pairdel(unsigned char *x, int *n, int *s, double *d,
 			  double *BF, int *variance, double *var,
 			  int *gamma, double *alpha)
 {
-    int i1, i2, k, Nd, Ns1, Ns2, L, target, s1, s2;
-    double P1, P2, Q, A, B, C, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b;
+    int i1, i2, Nd, Ns1, Ns2, L, target, s1, s2;
+    double P1, P2, Q, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b;
 
     PREPARE_BF_TN93
 
@@ -812,8 +812,8 @@ void distDNA_BH87(unsigned char *x, int *n, int *s, double *d,
    currently not available.
    </FIXME> */
 {
-    int i1, i2, k, kb, s1, s2, m, Ntab[16], ROWsums[4], ndim = 4, info, ipiv[16];
-    double P12[16], P21[16], U[16];
+    int i1, i2, k, kb, s1, s2, m, Ntab[16], ROWsums[4];
+    double P12[16], P21[16];
 
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
@@ -987,33 +987,62 @@ void distDNA_ParaLin_pairdel(unsigned char *x, int *n, int *s, double *d,
 /* 	} */
 /*     } */
 /*     if (! *freq) for (i = 0; i < 4; i++) BF[i] /= m; */
-}
+/* } */
+
+/* void BaseProportion(unsigned char *x, int *n, double *BF) */
+/* { */
+/*     int i; */
+
+/*     for (i = 0; i < *n; i++) { */
+/* 	    switch (x[i]) { */
+/* 	    case 136 : BF[0]++; break; */
+/* 	    case 40 : BF[1]++; break; */
+/* 	    case 72 : BF[2]++; break; */
+/* 	    case 24 : BF[3]++; break; */
+/* 	    case 192 : BF[4]++; break; */
+/* 	    case 160 : BF[5]++; break; */
+/* 	    case 144 : BF[6]++; break; */
+/* 	    case 96 : BF[7]++; break; */
+/* 	    case 80 : BF[8]++; break; */
+/* 	    case 48 : BF[9]++; break; */
+/* 	    case 224 : BF[10]++; break; */
+/* 	    case 176 : BF[11]++; break; */
+/* 	    case 208 : BF[12]++; break; */
+/* 	    case 112 : BF[13]++; break; */
+/* 	    case 240 : BF[14]++; break; */
+/* 	    case 4 : BF[15]++; break; */
+/* 	    case 2 : BF[16]++; break; */
+/* 	    } */
+/*     } */
+/* } */
 
+/* a hash table is much faster than switch (2012-01-10) */
 void BaseProportion(unsigned char *x, int *n, double *BF)
 {
-    int i;
-
-    for (i = 0; i < *n; i++) {
-	    switch (x[i]) {
-	    case 136 : BF[0]++; break;
-	    case 40 : BF[1]++; break;
-	    case 72 : BF[2]++; break;
-	    case 24 : BF[3]++; break;
-	    case 192 : BF[4]++; break;
-	    case 160 : BF[5]++; break;
-	    case 144 : BF[6]++; break;
-	    case 96 : BF[7]++; break;
-	    case 80 : BF[8]++; break;
-	    case 48 : BF[9]++; break;
-	    case 224 : BF[10]++; break;
-	    case 176 : BF[11]++; break;
-	    case 208 : BF[12]++; break;
-	    case 112 : BF[13]++; break;
-	    case 240 : BF[14]++; break;
-	    case 4 : BF[15]++; break;
-	    case 2 : BF[16]++; break;
-	    }
-    }
+	int i;
+	double count[256];
+
+	memset(count, 0, 256*sizeof(double));
+
+	for (i = 0; i < *n; i++) count[x[i]]++;
+
+	BF[0] = count[136];
+	BF[1] = count[40];
+	BF[2] = count[72];
+	BF[3] = count[24];
+	BF[4] = count[192];
+	BF[5] = count[160];
+	BF[6] = count[144];
+	BF[7] = count[96];
+	BF[8] = count[80];
+	BF[9] = count[48];
+	BF[10] = count[224];
+	BF[11] = count[176];
+	BF[12] = count[208];
+	BF[13] = count[112];
+	BF[14] = count[240];
+	BF[15] = count[4];
+	BF[16] = count[2];
 }
 
 void SegSites(unsigned char *x, int *n, int *s, int *seg)