source: branches/items/ISLAND_HOPPING/trnsprob.cxx

Last change on this file was 10148, checked in by westram, 12 years ago
  • makedepend changed behavior between 1.0.2 and 1.0.4
    • 1.0.4 does generate impossible dependencies (see [10143] for an example)
    • force correct behavior by defining SIMPLE_ARB_ASSERT on CLI (when invoking compiler and makedepend)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.9 KB
Line 
1// =============================================================
2/*                                                               */
3//   File      : trnsprob.c
4//   Purpose   :
5/*                                                               */
6//   Institute of Microbiology (Technical University Munich)
7//   http://www.arb-home.de/
8/*                                                               */
9// =============================================================
10
11#include "trnsprob.h"
12#include "defs.h"
13#include "mem.h"
14
15#include <arb_simple_assert.h>
16#include <stdio.h>
17
18#define ih_assert(bed) arb_assert(bed)
19#define EPS 0.00001
20
21//============================================================================
22
23static void identity(double **i) {
24    i[T][T]=1; i[T][C]=0; i[T][A]=0; i[T][G]=0;
25    i[C][T]=0; i[C][C]=1; i[C][A]=0; i[C][G]=0;
26    i[A][T]=0; i[A][C]=0; i[A][A]=1; i[A][G]=0;
27    i[G][T]=0; i[G][C]=0; i[G][A]=0; i[G][G]=1;
28}
29
30//............................................................................
31
32static void copy(double **i,double **j) {
33    i[T][T]=j[T][T]; i[T][C]=j[T][C]; i[T][A]=j[T][A]; i[T][G]=j[T][G];
34    i[C][T]=j[C][T]; i[C][C]=j[C][C]; i[C][A]=j[C][A]; i[C][G]=j[C][G];
35    i[A][T]=j[A][T]; i[A][C]=j[A][C]; i[A][A]=j[A][A]; i[A][G]=j[A][G];
36    i[G][T]=j[G][T]; i[G][C]=j[G][C]; i[G][A]=j[G][A]; i[G][G]=j[G][G];
37}
38
39//............................................................................
40
41static void ipol(double **i,double **j,double **k,double f) { double g; g=1.0-f;
42    i[T][T]=g*j[T][T]+f*k[T][T]; i[T][C]=g*j[T][C]+f*k[T][C]; i[T][A]=g*j[T][A]+f*k[T][A]; i[T][G]=g*j[T][G]+f*k[T][G];
43    i[C][T]=g*j[C][T]+f*k[C][T]; i[C][C]=g*j[C][C]+f*k[C][C]; i[C][A]=g*j[C][A]+f*k[C][A]; i[C][G]=g*j[C][G]+f*k[C][G];
44    i[A][T]=g*j[A][T]+f*k[A][T]; i[A][C]=g*j[A][C]+f*k[A][C]; i[A][A]=g*j[A][A]+f*k[A][A]; i[A][G]=g*j[A][G]+f*k[A][G];
45    i[G][T]=g*j[G][T]+f*k[G][T]; i[G][C]=g*j[G][C]+f*k[G][C]; i[G][A]=g*j[G][A]+f*k[G][A]; i[G][G]=g*j[G][G]+f*k[G][G];
46}
47
48//............................................................................
49
50static void addmul(double **i,double **j,double f) {
51    i[T][T]+=j[T][T]*f; i[T][C]+=j[T][C]*f; i[T][A]+=j[T][A]*f; i[T][G]+=j[T][G]*f;
52    i[C][T]+=j[C][T]*f; i[C][C]+=j[C][C]*f; i[C][A]+=j[C][A]*f; i[C][G]+=j[C][G]*f;
53    i[A][T]+=j[A][T]*f; i[A][C]+=j[A][C]*f; i[A][A]+=j[A][A]*f; i[A][G]+=j[A][G]*f;
54    i[G][T]+=j[G][T]*f; i[G][C]+=j[G][C]*f; i[G][A]+=j[G][A]*f; i[G][G]+=j[G][G]*f;
55}
56
57//............................................................................
58
59static void dot(double **i,double **j,double **k) {
60    i[T][T]=j[T][T]*k[T][T]+j[T][C]*k[C][T]+j[T][A]*k[A][T]+j[T][G]*k[G][T];
61    i[T][C]=j[T][T]*k[T][C]+j[T][C]*k[C][C]+j[T][A]*k[A][C]+j[T][G]*k[G][C];
62    i[T][A]=j[T][T]*k[T][A]+j[T][C]*k[C][A]+j[T][A]*k[A][A]+j[T][G]*k[G][A];
63    i[T][G]=j[T][T]*k[T][G]+j[T][C]*k[C][G]+j[T][A]*k[A][G]+j[T][G]*k[G][G];
64    i[C][T]=j[C][T]*k[T][T]+j[C][C]*k[C][T]+j[C][A]*k[A][T]+j[C][G]*k[G][T];
65    i[C][C]=j[C][T]*k[T][C]+j[C][C]*k[C][C]+j[C][A]*k[A][C]+j[C][G]*k[G][C];
66    i[C][A]=j[C][T]*k[T][A]+j[C][C]*k[C][A]+j[C][A]*k[A][A]+j[C][G]*k[G][A];
67    i[C][G]=j[C][T]*k[T][G]+j[C][C]*k[C][G]+j[C][A]*k[A][G]+j[C][G]*k[G][G];
68    i[A][T]=j[A][T]*k[T][T]+j[A][C]*k[C][T]+j[A][A]*k[A][T]+j[A][G]*k[G][T];
69    i[A][C]=j[A][T]*k[T][C]+j[A][C]*k[C][C]+j[A][A]*k[A][C]+j[A][G]*k[G][C];
70    i[A][A]=j[A][T]*k[T][A]+j[A][C]*k[C][A]+j[A][A]*k[A][A]+j[A][G]*k[G][A];
71    i[A][G]=j[A][T]*k[T][G]+j[A][C]*k[C][G]+j[A][A]*k[A][G]+j[A][G]*k[G][G];
72    i[G][T]=j[G][T]*k[T][T]+j[G][C]*k[C][T]+j[G][A]*k[A][T]+j[G][G]*k[G][T];
73    i[G][C]=j[G][T]*k[T][C]+j[G][C]*k[C][C]+j[G][A]*k[A][C]+j[G][G]*k[G][C];
74    i[G][A]=j[G][T]*k[T][A]+j[G][C]*k[C][A]+j[G][A]*k[A][A]+j[G][G]*k[G][A];
75    i[G][G]=j[G][T]*k[T][G]+j[G][C]*k[C][G]+j[G][A]*k[A][G]+j[G][G]*k[G][G];
76}
77
78//============================================================================
79
80char encodeBase(char c) {
81
82    switch(c) {
83        case 'U': return T;
84        case 'T': return T;
85        case 'C': return C;
86        case 'A': return A;
87        case 'G': return G;
88        default : return N;
89    }
90
91}
92
93//............................................................................
94
95char decodeBase(char c) {
96
97    switch(c) {
98        case T: return 'T';
99        case C: return 'C';
100        case A: return 'A';
101        case G: return 'G';
102        case N: return '?';
103        default : return '-';
104    }
105
106}
107
108//============================================================================
109
110void normalizeBaseFreqs(
111                        double *F,double fT,double fC,double fA,double fG
112                        ) {
113    double s;
114
115    s=fT+fC+fA+fG;
116
117    if(s<REAL_MIN) {fT=1.; fC=1.; fA=1.; fG=1.; s=4.;}
118
119    fT/=s; fC/=s; fA/=s; fG/=s;
120
121    if(fT<ATOLPARAM) fT=ATOLPARAM;
122    if(fC<ATOLPARAM) fC=ATOLPARAM;
123    if(fA<ATOLPARAM) fA=ATOLPARAM;
124    if(fG<ATOLPARAM) fG=ATOLPARAM;
125
126    F[T]=fT; F[C]=fC; F[A]=fA; F[G]=fG;
127
128}
129
130//............................................................................
131
132void normalizeRateParams(
133                         double *X,double a,double b,double c,double d,double e,double f
134                         ) { // TC TA TG CA CG AG
135    double s;
136
137    s=a+b+c+d+e+f;
138
139    if(s<REAL_MIN) {a=1.; b=1.; c=1.; d=1.; e=1.; s=6.;}
140
141    a/=s; b/=s; c/=s; d/=s; e/=s;
142
143    if(a<ATOLPARAM) a=ATOLPARAM;
144    if(b<ATOLPARAM) b=ATOLPARAM;
145    if(c<ATOLPARAM) c=ATOLPARAM;
146    if(d<ATOLPARAM) d=ATOLPARAM;
147    if(e<ATOLPARAM) e=ATOLPARAM;
148
149    X[0]=a; X[1]=b; X[2]=c; X[3]=d; X[4]=e;
150
151}
152
153//============================================================================
154
155void initTrnsprob(double ****PP) {
156    double ***P; short j;
157
158    P=(double***)newVector(128,sizeof(double **));
159    for(j=0;j<128;j++) P[j]=(double **)newMatrix(N,N,sizeof(double));
160
161    *PP=P;
162
163}
164
165//............................................................................
166
167void uninitTrnsprob(double ****PP) {
168    double ***P; short j;
169
170    P=*PP;
171
172    for(j=0;j<128;j++) freeMatrix(&P[j]);
173    freeBlock(PP);
174
175}
176
177//============================================================================
178
179void getTrnsprob(double **P,double ***PP,double d) {
180    long I,J;
181    double *Y[N],YT[N],YC[N],YA[N],YG[N];
182    double *Z[N],ZT[N],ZC[N],ZA[N],ZG[N];
183    double *X[N],XT[N],XC[N],XA[N],XG[N];
184
185    d/=EPS; I=(long)d;
186
187    if( I < 0        ) {copy(P,PP[0]); return;}
188
189    X[T]=XT; X[C]=XC; X[A]=XA; X[G]=XG;
190
191    if( I < 32       ) {ipol(X,PP[I],PP[I==31?33:I+1],d-I); copy(P,X); return;}
192
193    if( I < 1024     ) {
194        J=I%32;     ipol(X,PP[J],PP[J==31?33:J+1],d-I);
195        dot(P,PP[ I/32    + 32 ],X);
196        return;
197    }
198
199    Y[T]=YT; Y[C]=YC; Y[A]=YA; Y[G]=YG;
200
201    if( I < 32768   ) {
202        J=I%32;     ipol(X,PP[J],PP[J==31?33:J+1],d-I);
203        J=I%1024;    dot(Y,PP[ J/32    + 32 ],X);
204        dot(P,PP[ I/1024  + 64 ],Y);
205        return;
206    }
207
208    Z[T]=ZT; Z[C]=ZC; Z[A]=ZA; Z[G]=ZG;
209
210    if( I < 1048576 ) {
211        J=I%32;     ipol(X,PP[J],PP[J==31?33:J+1],d-I);
212        J=I%1024;    dot(Y,PP[ J/32    + 32 ],X);
213        J=I%32768;   dot(Z,PP[ J/1024  + 64 ],Y);
214        dot(P,PP[ I/32768 + 96 ],Z);
215        return;
216    }
217
218    dot(Y,PP[ 31      + 32 ],PP[31]);
219    dot(Z,PP[ 31      + 64 ],Y);
220    dot(P,PP[ 31      + 96 ],Z);
221
222}
223
224//============================================================================
225
226void updateTrnsprob(double ***PP,double *F,double *X,short M) {
227    double s,a,b,c,d,e,f;
228    double *Q[N],QT[N],QC[N],QA[N],QG[N];
229    double *R[N],RT[N],RC[N],RA[N],RG[N];
230    double *S[N],ST[N],SC[N],SA[N],SG[N];
231
232    Q[T]=QT; Q[C]=QC; Q[A]=QA; Q[G]=QG;
233    R[T]=RT; R[C]=RC; R[A]=RA; R[G]=RG;
234    S[T]=ST; S[C]=SC; S[A]=SA; S[G]=SG;
235
236    switch(M) {
237        case REV:   a = X[0]; b=X[1]; c=X[2]; d=X[3]; e=X[4]; f=1.0-(a+b+c+d+e); break;
238        case TN93:  a = X[0]; b=c=d=e=X[1]; f=1.0-(a+b+c+d+e); break;
239        case HKY85: a = f=X[0]; b=c=d=e=0.25*(1.0-(a+f)); break;
240        default : ih_assert(0); break;
241    }
242
243    s=0.5/(a*F[T]*F[C]+b*F[T]*F[A]+c*F[T]*F[G]+d*F[C]*F[A]+e*F[C]*F[G]+f*F[A]*F[G]);
244
245    Q[T][C]=a*F[C]*s; Q[T][A]=b*F[A]*s; Q[T][G]=c*F[G]*s;
246    Q[C][T]=a*F[T]*s;                   Q[C][A]=d*F[A]*s; Q[C][G]=e*F[G]*s;
247    Q[A][T]=b*F[T]*s; Q[A][C]=d*F[C]*s;                   Q[A][G]=f*F[G]*s;
248    Q[G][T]=c*F[T]*s; Q[G][C]=e*F[C]*s; Q[G][A]=f*F[A]*s;
249
250    Q[T][T] = - (           Q[T][C] + Q[T][A] + Q[T][G] );
251    Q[C][C] = - ( Q[C][T] +           Q[C][A] + Q[C][G] );
252    Q[A][A] = - ( Q[A][T] + Q[A][C] +           Q[A][G] );
253    Q[G][G] = - ( Q[G][T] + Q[G][C] + Q[G][A]           );
254
255    identity(PP[1]);
256    addmul(PP[1],Q,EPS);
257    dot(R,Q,Q);
258    addmul(PP[1],R,EPS*EPS/2.0);
259    dot(S,R,Q);
260    addmul(PP[1],S,EPS*EPS*EPS/6.0);
261    dot(S,R,R);
262    addmul(PP[1],S,EPS*EPS*EPS*EPS/24.0);
263
264    identity(PP[ 0]);
265    identity(PP[32]);
266    identity(PP[64]);
267    identity(PP[96]);
268
269    dot(PP[  2      ],PP[  1      ],PP[  1      ]);
270    dot(PP[  4      ],PP[  2      ],PP[  2      ]);
271    dot(PP[  8      ],PP[  4      ],PP[  4      ]);
272    dot(PP[ 16      ],PP[  8      ],PP[  8      ]);
273    dot(PP[  1 + 32 ],PP[ 16      ],PP[ 16      ]);
274    dot(PP[  2 + 32 ],PP[  1 + 32 ],PP[  1 + 32 ]);
275    dot(PP[  4 + 32 ],PP[  2 + 32 ],PP[  2 + 32 ]);
276    dot(PP[  8 + 32 ],PP[  4 + 32 ],PP[  4 + 32 ]);
277    dot(PP[ 16 + 32 ],PP[  8 + 32 ],PP[  8 + 32 ]);
278    dot(PP[  1 + 64 ],PP[ 16 + 32 ],PP[ 16 + 32 ]);
279    dot(PP[  2 + 64 ],PP[  1 + 64 ],PP[  1 + 64 ]);
280    dot(PP[  4 + 64 ],PP[  2 + 64 ],PP[  2 + 64 ]);
281    dot(PP[  8 + 64 ],PP[  4 + 64 ],PP[  4 + 64 ]);
282    dot(PP[ 16 + 64 ],PP[  8 + 64 ],PP[  8 + 64 ]);
283    dot(PP[  1 + 96 ],PP[ 16 + 64 ],PP[ 16 + 64 ]);
284    dot(PP[  2 + 96 ],PP[  1 + 96 ],PP[  1 + 96 ]);
285    dot(PP[  4 + 96 ],PP[  2 + 96 ],PP[  2 + 96 ]);
286    dot(PP[  8 + 96 ],PP[  4 + 96 ],PP[  4 + 96 ]);
287    dot(PP[ 16 + 96 ],PP[  8 + 96 ],PP[  8 + 96 ]);
288
289    dot(PP[  3      ],PP[  2      ],PP[  1      ]);
290    dot(PP[  6      ],PP[  3      ],PP[  3      ]);
291    dot(PP[ 12      ],PP[  6      ],PP[  6      ]);
292    dot(PP[ 24      ],PP[ 12      ],PP[ 12      ]);
293    dot(PP[  3 + 32 ],PP[  2 + 32 ],PP[  1 + 32 ]);
294    dot(PP[  6 + 32 ],PP[  3 + 32 ],PP[  3 + 32 ]);
295    dot(PP[ 12 + 32 ],PP[  6 + 32 ],PP[  6 + 32 ]);
296    dot(PP[ 24 + 32 ],PP[ 12 + 32 ],PP[ 12 + 32 ]);
297    dot(PP[  3 + 64 ],PP[  2 + 64 ],PP[  1 + 64 ]);
298    dot(PP[  6 + 64 ],PP[  3 + 64 ],PP[  3 + 64 ]);
299    dot(PP[ 12 + 64 ],PP[  6 + 64 ],PP[  6 + 64 ]);
300    dot(PP[ 24 + 64 ],PP[ 12 + 64 ],PP[ 12 + 64 ]);
301    dot(PP[  3 + 96 ],PP[  2 + 96 ],PP[  1 + 96 ]);
302    dot(PP[  6 + 96 ],PP[  3 + 96 ],PP[  3 + 96 ]);
303    dot(PP[ 12 + 96 ],PP[  6 + 96 ],PP[  6 + 96 ]);
304    dot(PP[ 24 + 96 ],PP[ 12 + 96 ],PP[ 12 + 96 ]);
305
306    dot(PP[  5      ],PP[  3      ],PP[  2      ]);
307    dot(PP[ 10      ],PP[  5      ],PP[  5      ]);
308    dot(PP[ 20      ],PP[ 10      ],PP[ 10      ]);
309    dot(PP[  5 + 32 ],PP[  3 + 32 ],PP[  2 + 32 ]);
310    dot(PP[ 10 + 32 ],PP[  5 + 32 ],PP[  5 + 32 ]);
311    dot(PP[ 20 + 32 ],PP[ 10 + 32 ],PP[ 10 + 32 ]);
312    dot(PP[  5 + 64 ],PP[  3 + 64 ],PP[  2 + 64 ]);
313    dot(PP[ 10 + 64 ],PP[  5 + 64 ],PP[  5 + 64 ]);
314    dot(PP[ 20 + 64 ],PP[ 10 + 64 ],PP[ 10 + 64 ]);
315    dot(PP[  5 + 96 ],PP[  3 + 96 ],PP[  2 + 96 ]);
316    dot(PP[ 10 + 96 ],PP[  5 + 96 ],PP[  5 + 96 ]);
317    dot(PP[ 20 + 96 ],PP[ 10 + 96 ],PP[ 10 + 96 ]);
318
319    dot(PP[  7      ],PP[  4      ],PP[  3      ]);
320    dot(PP[ 14      ],PP[  7      ],PP[  7      ]);
321    dot(PP[ 28      ],PP[ 14      ],PP[ 14      ]);
322    dot(PP[  7 + 32 ],PP[  4 + 32 ],PP[  3 + 32 ]);
323    dot(PP[ 14 + 32 ],PP[  7 + 32 ],PP[  7 + 32 ]);
324    dot(PP[ 28 + 32 ],PP[ 14 + 32 ],PP[ 14 + 32 ]);
325    dot(PP[  7 + 64 ],PP[  4 + 64 ],PP[  3 + 64 ]);
326    dot(PP[ 14 + 64 ],PP[  7 + 64 ],PP[  7 + 64 ]);
327    dot(PP[ 28 + 64 ],PP[ 14 + 64 ],PP[ 14 + 64 ]);
328    dot(PP[  7 + 96 ],PP[  4 + 96 ],PP[  3 + 96 ]);
329    dot(PP[ 14 + 96 ],PP[  7 + 96 ],PP[  7 + 96 ]);
330    dot(PP[ 28 + 96 ],PP[ 14 + 96 ],PP[ 14 + 96 ]);
331
332    dot(PP[  9      ],PP[  5      ],PP[  4      ]);
333    dot(PP[ 18      ],PP[  9      ],PP[  9      ]);
334    dot(PP[  9 + 32 ],PP[  5 + 32 ],PP[  4 + 32 ]);
335    dot(PP[ 18 + 32 ],PP[  9 + 32 ],PP[  9 + 32 ]);
336    dot(PP[  9 + 64 ],PP[  5 + 64 ],PP[  4 + 64 ]);
337    dot(PP[ 18 + 64 ],PP[  9 + 64 ],PP[  9 + 64 ]);
338    dot(PP[  9 + 96 ],PP[  5 + 96 ],PP[  4 + 96 ]);
339    dot(PP[ 18 + 96 ],PP[  9 + 96 ],PP[  9 + 96 ]);
340
341    dot(PP[ 11      ],PP[  6      ],PP[  5      ]);
342    dot(PP[ 22      ],PP[ 11      ],PP[ 11      ]);
343    dot(PP[ 11 + 32 ],PP[  6 + 32 ],PP[  5 + 32 ]);
344    dot(PP[ 22 + 32 ],PP[ 11 + 32 ],PP[ 11 + 32 ]);
345    dot(PP[ 11 + 64 ],PP[  6 + 64 ],PP[  5 + 64 ]);
346    dot(PP[ 22 + 64 ],PP[ 11 + 64 ],PP[ 11 + 64 ]);
347    dot(PP[ 11 + 96 ],PP[  6 + 96 ],PP[  5 + 96 ]);
348    dot(PP[ 22 + 96 ],PP[ 11 + 96 ],PP[ 11 + 96 ]);
349
350    dot(PP[ 13      ],PP[  7      ],PP[  6      ]);
351    dot(PP[ 26      ],PP[ 13      ],PP[ 13      ]);
352    dot(PP[ 13 + 32 ],PP[  7 + 32 ],PP[  6 + 32 ]);
353    dot(PP[ 26 + 32 ],PP[ 13 + 32 ],PP[ 13 + 32 ]);
354    dot(PP[ 13 + 64 ],PP[  7 + 64 ],PP[  6 + 64 ]);
355    dot(PP[ 26 + 64 ],PP[ 13 + 64 ],PP[ 13 + 64 ]);
356    dot(PP[ 13 + 96 ],PP[  7 + 96 ],PP[  6 + 96 ]);
357    dot(PP[ 26 + 96 ],PP[ 13 + 96 ],PP[ 13 + 96 ]);
358
359    dot(PP[ 15      ],PP[  8      ],PP[  7      ]);
360    dot(PP[ 30      ],PP[ 15      ],PP[ 15      ]);
361    dot(PP[ 15 + 32 ],PP[  8 + 32 ],PP[  7 + 32 ]);
362    dot(PP[ 30 + 32 ],PP[ 15 + 32 ],PP[ 15 + 32 ]);
363    dot(PP[ 15 + 64 ],PP[  8 + 64 ],PP[  7 + 64 ]);
364    dot(PP[ 30 + 64 ],PP[ 15 + 64 ],PP[ 15 + 64 ]);
365    dot(PP[ 15 + 96 ],PP[  8 + 96 ],PP[  7 + 96 ]);
366    dot(PP[ 30 + 96 ],PP[ 15 + 96 ],PP[ 15 + 96 ]);
367
368    dot(PP[ 17      ],PP[  9      ],PP[  8      ]);
369    dot(PP[ 17 + 32 ],PP[  9 + 32 ],PP[  8 + 32 ]);
370    dot(PP[ 17 + 64 ],PP[  9 + 64 ],PP[  8 + 64 ]);
371    dot(PP[ 17 + 96 ],PP[  9 + 96 ],PP[  8 + 96 ]);
372
373    dot(PP[ 19      ],PP[ 10      ],PP[  9      ]);
374    dot(PP[ 19 + 32 ],PP[ 10 + 32 ],PP[  9 + 32 ]);
375    dot(PP[ 19 + 64 ],PP[ 10 + 64 ],PP[  9 + 64 ]);
376    dot(PP[ 19 + 96 ],PP[ 10 + 96 ],PP[  9 + 96 ]);
377
378    dot(PP[ 21      ],PP[ 11      ],PP[ 10      ]);
379    dot(PP[ 21 + 32 ],PP[ 11 + 32 ],PP[ 10 + 32 ]);
380    dot(PP[ 21 + 64 ],PP[ 11 + 64 ],PP[ 10 + 64 ]);
381    dot(PP[ 21 + 96 ],PP[ 11 + 96 ],PP[ 10 + 96 ]);
382
383    dot(PP[ 23      ],PP[ 12      ],PP[ 11      ]);
384    dot(PP[ 23 + 32 ],PP[ 12 + 32 ],PP[ 11 + 32 ]);
385    dot(PP[ 23 + 64 ],PP[ 12 + 64 ],PP[ 11 + 64 ]);
386    dot(PP[ 23 + 96 ],PP[ 12 + 96 ],PP[ 11 + 96 ]);
387
388    dot(PP[ 25      ],PP[ 13      ],PP[ 12      ]);
389    dot(PP[ 25 + 32 ],PP[ 13 + 32 ],PP[ 12 + 32 ]);
390    dot(PP[ 25 + 64 ],PP[ 13 + 64 ],PP[ 12 + 64 ]);
391    dot(PP[ 25 + 96 ],PP[ 13 + 96 ],PP[ 12 + 96 ]);
392
393    dot(PP[ 27      ],PP[ 14      ],PP[ 13      ]);
394    dot(PP[ 27 + 32 ],PP[ 14 + 32 ],PP[ 13 + 32 ]);
395    dot(PP[ 27 + 64 ],PP[ 14 + 64 ],PP[ 13 + 64 ]);
396    dot(PP[ 27 + 96 ],PP[ 14 + 96 ],PP[ 13 + 96 ]);
397
398    dot(PP[ 29      ],PP[ 15      ],PP[ 14      ]);
399    dot(PP[ 29 + 32 ],PP[ 15 + 32 ],PP[ 14 + 32 ]);
400    dot(PP[ 29 + 64 ],PP[ 15 + 64 ],PP[ 14 + 64 ]);
401    dot(PP[ 29 + 96 ],PP[ 15 + 96 ],PP[ 14 + 96 ]);
402
403    dot(PP[ 31      ],PP[ 16      ],PP[ 15      ]);
404    dot(PP[ 31 + 32 ],PP[ 16 + 32 ],PP[ 15 + 32 ]);
405    dot(PP[ 31 + 64 ],PP[ 16 + 64 ],PP[ 15 + 64 ]);
406    dot(PP[ 31 + 96 ],PP[ 16 + 96 ],PP[ 15 + 96 ]);
407
408}
409
Note: See TracBrowser for help on using the repository browser.