1 #include "simint/boys/boys.h"
2 #include "simint/ostei/gen/ostei_generated.h"
3 #include "simint/vectorization/vectorization.h"
4 #include <math.h>
5 #include <string.h>
6 
7 
ostei_d_d_d_p(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__d_d_d_p)8 int ostei_d_d_d_p(struct simint_multi_shellpair const P,
9                   struct simint_multi_shellpair const Q,
10                   double screen_tol,
11                   double * const restrict work,
12                   double * const restrict INT__d_d_d_p)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__d_d_d_p);
17     int ab, cd, abcd;
18     int istart, jstart;
19     int iprimcd, nprim_icd, icd;
20     const int check_screen = (screen_tol > 0.0);
21     int i, j;
22     int n;
23     int not_screened;
24     int real_abcd;
25     int iket;
26     int ibra;
27 
28     // partition workspace
29     double * const INT__d_s_d_s = work + (SIMINT_NSHELL_SIMD * 0);
30     double * const INT__d_s_f_s = work + (SIMINT_NSHELL_SIMD * 36);
31     double * const INT__f_s_d_s = work + (SIMINT_NSHELL_SIMD * 96);
32     double * const INT__f_s_f_s = work + (SIMINT_NSHELL_SIMD * 156);
33     double * const INT__g_s_d_s = work + (SIMINT_NSHELL_SIMD * 256);
34     double * const INT__g_s_f_s = work + (SIMINT_NSHELL_SIMD * 346);
35     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*496);
36     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
37     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_p_s = primwork + 8;
38     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_s = primwork + 17;
39     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_p_s = primwork + 38;
40     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_d_s = primwork + 65;
41     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 101;
42     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_p_s = primwork + 137;
43     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_d_s = primwork + 191;
44     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_f_s = primwork + 263;
45     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_s = primwork + 323;
46     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_p_s = primwork + 373;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_d_s = primwork + 463;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_f_s = primwork + 583;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_s_s = primwork + 683;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_p_s = primwork + 743;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_d_s = primwork + 878;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_f_s = primwork + 1058;
53     double * const hrrwork = (double *)(primwork + 1208);
54     double * const HRR_INT__d_p_d_s = hrrwork + 0;
55     double * const HRR_INT__d_p_f_s = hrrwork + 108;
56     double * const HRR_INT__d_d_d_s = hrrwork + 288;
57     double * const HRR_INT__d_d_f_s = hrrwork + 504;
58     double * const HRR_INT__f_p_d_s = hrrwork + 864;
59     double * const HRR_INT__f_p_f_s = hrrwork + 1044;
60 
61 
62     // Create constants
63     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
64     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
65     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
66     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
67     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
68 
69 
70     ////////////////////////////////////////
71     // Loop over shells and primitives
72     ////////////////////////////////////////
73 
74     real_abcd = 0;
75     istart = 0;
76     for(ab = 0; ab < P.nshell12_clip; ++ab)
77     {
78         const int iend = istart + P.nprim12[ab];
79 
80         cd = 0;
81         jstart = 0;
82 
83         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
84         {
85             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
86             int jend = jstart;
87             for(i = 0; i < nshellbatch; i++)
88                 jend += Q.nprim12[cd+i];
89 
90             // Clear the beginning of the workspace (where we are accumulating integrals)
91             memset(work, 0, SIMINT_NSHELL_SIMD * 496 * sizeof(double));
92             abcd = 0;
93 
94 
95             for(i = istart; i < iend; ++i)
96             {
97                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
98 
99                 if(check_screen)
100                 {
101                     // Skip this whole thing if always insignificant
102                     if((P.screen[i] * Q.screen_max) < screen_tol)
103                         continue;
104                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
105                 }
106 
107                 icd = 0;
108                 iprimcd = 0;
109                 nprim_icd = Q.nprim12[cd];
110                 double * restrict PRIM_PTR_INT__d_s_d_s = INT__d_s_d_s + abcd * 36;
111                 double * restrict PRIM_PTR_INT__d_s_f_s = INT__d_s_f_s + abcd * 60;
112                 double * restrict PRIM_PTR_INT__f_s_d_s = INT__f_s_d_s + abcd * 60;
113                 double * restrict PRIM_PTR_INT__f_s_f_s = INT__f_s_f_s + abcd * 100;
114                 double * restrict PRIM_PTR_INT__g_s_d_s = INT__g_s_d_s + abcd * 90;
115                 double * restrict PRIM_PTR_INT__g_s_f_s = INT__g_s_f_s + abcd * 150;
116 
117 
118 
119                 // Load these one per loop over i
120                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
121                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
122                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
123 
124                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
125 
126                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
127                 {
128                     // calculate the shell offsets
129                     // these are the offset from the shell pointed to by cd
130                     // for each element
131                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
132                     int lastoffset = 0;
133                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
134 
135                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
136                     {
137                         // Handle if the first element of the vector is a new shell
138                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
139                         {
140                             nprim_icd += Q.nprim12[cd + (++icd)];
141                             PRIM_PTR_INT__d_s_d_s += 36;
142                             PRIM_PTR_INT__d_s_f_s += 60;
143                             PRIM_PTR_INT__f_s_d_s += 60;
144                             PRIM_PTR_INT__f_s_f_s += 100;
145                             PRIM_PTR_INT__g_s_d_s += 90;
146                             PRIM_PTR_INT__g_s_f_s += 150;
147                         }
148                         iprimcd++;
149                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
150                         {
151                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
152                             {
153                                 shelloffsets[n] = shelloffsets[n-1] + 1;
154                                 lastoffset++;
155                                 nprim_icd += Q.nprim12[cd + (++icd)];
156                             }
157                             else
158                                 shelloffsets[n] = shelloffsets[n-1];
159                             iprimcd++;
160                         }
161                     }
162                     else
163                         iprimcd += SIMINT_SIMD_LEN;
164 
165                     // Do we have to compute this vector (or has it been screened out)?
166                     // (not_screened != 0 means we have to do this vector)
167                     if(check_screen)
168                     {
169                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
170                         if(vmax < screen_tol)
171                         {
172                             PRIM_PTR_INT__d_s_d_s += lastoffset*36;
173                             PRIM_PTR_INT__d_s_f_s += lastoffset*60;
174                             PRIM_PTR_INT__f_s_d_s += lastoffset*60;
175                             PRIM_PTR_INT__f_s_f_s += lastoffset*100;
176                             PRIM_PTR_INT__g_s_d_s += lastoffset*90;
177                             PRIM_PTR_INT__g_s_f_s += lastoffset*150;
178                             continue;
179                         }
180                     }
181 
182                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
183                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
184                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
185                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
186 
187 
188                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
189                     SIMINT_DBLTYPE PQ[3];
190                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
191                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
192                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
193                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
194                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
195                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
196 
197                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
198                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
199                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
200                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
201                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
202                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
203                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
204 
205                     // NOTE: Minus sign!
206                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
207                     SIMINT_DBLTYPE aop_PQ[3];
208                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
209                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
210                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
211 
212                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
213                     SIMINT_DBLTYPE aoq_PQ[3];
214                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
215                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
216                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
217                     // Put a minus sign here so we don't have to in RR routines
218                     a_over_q = SIMINT_NEG(a_over_q);
219 
220 
221                     //////////////////////////////////////////////
222                     // Fjt function section
223                     // Maximum v value: 7
224                     //////////////////////////////////////////////
225                     // The parameter to the Fjt function
226                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
227 
228 
229                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
230 
231 
232                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 7);
233                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
234                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
235                     for(n = 0; n <= 7; n++)
236                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
237 
238                     //////////////////////////////////////////////
239                     // Primitive integrals: Vertical recurrance
240                     //////////////////////////////////////////////
241 
242                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
243                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
244                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
245                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
246                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
247                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
248                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
249                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
250                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
251 
252 
253 
254                     // Forming PRIM_INT__p_s_s_s[7 * 3];
255                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
256                     {
257 
258                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
259                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]);
260 
261                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
262                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_s[n * 3 + 1]);
263 
264                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
265                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_s[n * 3 + 2]);
266 
267                     }
268 
269 
270 
271                     // Forming PRIM_INT__d_s_s_s[6 * 6];
272                     for(n = 0; n < 6; ++n)  // loop over orders of auxiliary function
273                     {
274 
275                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
276                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_s_s[n * 6 + 0]);
277                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__d_s_s_s[n * 6 + 0]);
278 
279                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
280                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_s_s[n * 6 + 1]);
281 
282                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
283                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_s_s[n * 6 + 2]);
284 
285                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
286                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_s_s[n * 6 + 3]);
287                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__d_s_s_s[n * 6 + 3]);
288 
289                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
290                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_s_s[n * 6 + 4]);
291 
292                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
293                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_s_s[n * 6 + 5]);
294                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__d_s_s_s[n * 6 + 5]);
295 
296                     }
297 
298 
299 
300                     // Forming PRIM_INT__d_s_p_s[3 * 18];
301                     for(n = 0; n < 3; ++n)  // loop over orders of auxiliary function
302                     {
303 
304                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
305                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__d_s_p_s[n * 18 + 0]);
306                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_p_s[n * 18 + 0]);
307 
308                         PRIM_INT__d_s_p_s[n * 18 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
309                         PRIM_INT__d_s_p_s[n * 18 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__d_s_p_s[n * 18 + 1]);
310 
311                         PRIM_INT__d_s_p_s[n * 18 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
312                         PRIM_INT__d_s_p_s[n * 18 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__d_s_p_s[n * 18 + 2]);
313 
314                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 1]);
315                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__d_s_p_s[n * 18 + 3]);
316                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_p_s[n * 18 + 3]);
317 
318                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 1]);
319                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__d_s_p_s[n * 18 + 4]);
320                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_p_s[n * 18 + 4]);
321 
322                         PRIM_INT__d_s_p_s[n * 18 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
323                         PRIM_INT__d_s_p_s[n * 18 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__d_s_p_s[n * 18 + 5]);
324 
325                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 2]);
326                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 2], PRIM_INT__d_s_p_s[n * 18 + 6]);
327                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_p_s[n * 18 + 6]);
328 
329                         PRIM_INT__d_s_p_s[n * 18 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 2]);
330                         PRIM_INT__d_s_p_s[n * 18 + 7] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 2], PRIM_INT__d_s_p_s[n * 18 + 7]);
331 
332                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 2]);
333                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 2], PRIM_INT__d_s_p_s[n * 18 + 8]);
334                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_p_s[n * 18 + 8]);
335 
336                         PRIM_INT__d_s_p_s[n * 18 + 9] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
337                         PRIM_INT__d_s_p_s[n * 18 + 9] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__d_s_p_s[n * 18 + 9]);
338 
339                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
340                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__d_s_p_s[n * 18 + 10]);
341                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_p_s[n * 18 + 10]);
342 
343                         PRIM_INT__d_s_p_s[n * 18 + 11] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
344                         PRIM_INT__d_s_p_s[n * 18 + 11] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__d_s_p_s[n * 18 + 11]);
345 
346                         PRIM_INT__d_s_p_s[n * 18 + 12] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 4]);
347                         PRIM_INT__d_s_p_s[n * 18 + 12] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 4], PRIM_INT__d_s_p_s[n * 18 + 12]);
348 
349                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 4]);
350                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 4], PRIM_INT__d_s_p_s[n * 18 + 13]);
351                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_p_s[n * 18 + 13]);
352 
353                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 4]);
354                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 4], PRIM_INT__d_s_p_s[n * 18 + 14]);
355                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_p_s[n * 18 + 14]);
356 
357                         PRIM_INT__d_s_p_s[n * 18 + 15] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
358                         PRIM_INT__d_s_p_s[n * 18 + 15] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__d_s_p_s[n * 18 + 15]);
359 
360                         PRIM_INT__d_s_p_s[n * 18 + 16] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
361                         PRIM_INT__d_s_p_s[n * 18 + 16] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__d_s_p_s[n * 18 + 16]);
362 
363                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
364                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__d_s_p_s[n * 18 + 17]);
365                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_p_s[n * 18 + 17]);
366 
367                     }
368 
369 
370 
371                     // Forming PRIM_INT__p_s_p_s[3 * 9];
372                     for(n = 0; n < 3; ++n)  // loop over orders of auxiliary function
373                     {
374 
375                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
376                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 0]);
377                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_p_s[n * 9 + 0]);
378 
379                         PRIM_INT__p_s_p_s[n * 9 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
380                         PRIM_INT__p_s_p_s[n * 9 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 1]);
381 
382                         PRIM_INT__p_s_p_s[n * 9 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
383                         PRIM_INT__p_s_p_s[n * 9 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 2]);
384 
385                         PRIM_INT__p_s_p_s[n * 9 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 1]);
386                         PRIM_INT__p_s_p_s[n * 9 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 3]);
387 
388                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
389                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 4]);
390                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_p_s[n * 9 + 4]);
391 
392                         PRIM_INT__p_s_p_s[n * 9 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
393                         PRIM_INT__p_s_p_s[n * 9 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 5]);
394 
395                         PRIM_INT__p_s_p_s[n * 9 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 2]);
396                         PRIM_INT__p_s_p_s[n * 9 + 6] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 6]);
397 
398                         PRIM_INT__p_s_p_s[n * 9 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 2]);
399                         PRIM_INT__p_s_p_s[n * 9 + 7] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 7]);
400 
401                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
402                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 8]);
403                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_p_s[n * 9 + 8]);
404 
405                     }
406 
407 
408                     VRR_K_d_s_d_s(
409                             PRIM_INT__d_s_d_s,
410                             PRIM_INT__d_s_p_s,
411                             PRIM_INT__d_s_s_s,
412                             PRIM_INT__p_s_p_s,
413                             Q_PA,
414                             a_over_q,
415                             aoq_PQ,
416                             one_over_2pq,
417                             one_over_2q,
418                             2);
419 
420 
421 
422                     // Forming PRIM_INT__s_s_p_s[3 * 3];
423                     for(n = 0; n < 3; ++n)  // loop over orders of auxiliary function
424                     {
425 
426                         PRIM_INT__s_s_p_s[n * 3 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
427                         PRIM_INT__s_s_p_s[n * 3 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_p_s[n * 3 + 0]);
428 
429                         PRIM_INT__s_s_p_s[n * 3 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
430                         PRIM_INT__s_s_p_s[n * 3 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_p_s[n * 3 + 1]);
431 
432                         PRIM_INT__s_s_p_s[n * 3 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
433                         PRIM_INT__s_s_p_s[n * 3 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_p_s[n * 3 + 2]);
434 
435                     }
436 
437 
438 
439                     // Forming PRIM_INT__p_s_d_s[2 * 18];
440                     for(n = 0; n < 2; ++n)  // loop over orders of auxiliary function
441                     {
442 
443                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_p_s[n * 9 + 0]);
444                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_p_s[(n+1) * 9 + 0], PRIM_INT__p_s_d_s[n * 18 + 0]);
445                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]), PRIM_INT__p_s_d_s[n * 18 + 0]);
446                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_d_s[n * 18 + 0]);
447 
448                         PRIM_INT__p_s_d_s[n * 18 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_p_s[n * 9 + 0]);
449                         PRIM_INT__p_s_d_s[n * 18 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_p_s[(n+1) * 9 + 0], PRIM_INT__p_s_d_s[n * 18 + 1]);
450 
451                         PRIM_INT__p_s_d_s[n * 18 + 3] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_p_s[n * 9 + 1]);
452                         PRIM_INT__p_s_d_s[n * 18 + 3] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_p_s[(n+1) * 9 + 1], PRIM_INT__p_s_d_s[n * 18 + 3]);
453                         PRIM_INT__p_s_d_s[n * 18 + 3] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]), PRIM_INT__p_s_d_s[n * 18 + 3]);
454 
455                         PRIM_INT__p_s_d_s[n * 18 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_p_s[n * 9 + 2]);
456                         PRIM_INT__p_s_d_s[n * 18 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_p_s[(n+1) * 9 + 2], PRIM_INT__p_s_d_s[n * 18 + 5]);
457                         PRIM_INT__p_s_d_s[n * 18 + 5] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]), PRIM_INT__p_s_d_s[n * 18 + 5]);
458 
459                         PRIM_INT__p_s_d_s[n * 18 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_p_s[n * 9 + 3]);
460                         PRIM_INT__p_s_d_s[n * 18 + 6] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_p_s[(n+1) * 9 + 3], PRIM_INT__p_s_d_s[n * 18 + 6]);
461                         PRIM_INT__p_s_d_s[n * 18 + 6] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_s[n * 3 + 1]), PRIM_INT__p_s_d_s[n * 18 + 6]);
462 
463                         PRIM_INT__p_s_d_s[n * 18 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_p_s[n * 9 + 3]);
464                         PRIM_INT__p_s_d_s[n * 18 + 7] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_p_s[(n+1) * 9 + 3], PRIM_INT__p_s_d_s[n * 18 + 7]);
465                         PRIM_INT__p_s_d_s[n * 18 + 7] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_d_s[n * 18 + 7]);
466 
467                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_p_s[n * 9 + 4]);
468                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_p_s[(n+1) * 9 + 4], PRIM_INT__p_s_d_s[n * 18 + 9]);
469                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_s[n * 3 + 1]), PRIM_INT__p_s_d_s[n * 18 + 9]);
470                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_d_s[n * 18 + 9]);
471 
472                         PRIM_INT__p_s_d_s[n * 18 + 11] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_p_s[n * 9 + 5]);
473                         PRIM_INT__p_s_d_s[n * 18 + 11] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_p_s[(n+1) * 9 + 5], PRIM_INT__p_s_d_s[n * 18 + 11]);
474                         PRIM_INT__p_s_d_s[n * 18 + 11] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_s[n * 3 + 1]), PRIM_INT__p_s_d_s[n * 18 + 11]);
475 
476                         PRIM_INT__p_s_d_s[n * 18 + 12] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_p_s[n * 9 + 6]);
477                         PRIM_INT__p_s_d_s[n * 18 + 12] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_p_s[(n+1) * 9 + 6], PRIM_INT__p_s_d_s[n * 18 + 12]);
478                         PRIM_INT__p_s_d_s[n * 18 + 12] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_s[n * 3 + 2]), PRIM_INT__p_s_d_s[n * 18 + 12]);
479 
480                         PRIM_INT__p_s_d_s[n * 18 + 13] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_p_s[n * 9 + 6]);
481                         PRIM_INT__p_s_d_s[n * 18 + 13] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_p_s[(n+1) * 9 + 6], PRIM_INT__p_s_d_s[n * 18 + 13]);
482 
483                         PRIM_INT__p_s_d_s[n * 18 + 15] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_p_s[n * 9 + 7]);
484                         PRIM_INT__p_s_d_s[n * 18 + 15] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_p_s[(n+1) * 9 + 7], PRIM_INT__p_s_d_s[n * 18 + 15]);
485                         PRIM_INT__p_s_d_s[n * 18 + 15] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_s[n * 3 + 2]), PRIM_INT__p_s_d_s[n * 18 + 15]);
486 
487                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_p_s[n * 9 + 8]);
488                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_p_s[(n+1) * 9 + 8], PRIM_INT__p_s_d_s[n * 18 + 17]);
489                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_s[n * 3 + 2]), PRIM_INT__p_s_d_s[n * 18 + 17]);
490                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_d_s[n * 18 + 17]);
491 
492                     }
493 
494 
495                     VRR_K_d_s_f_s(
496                             PRIM_INT__d_s_f_s,
497                             PRIM_INT__d_s_d_s,
498                             PRIM_INT__d_s_p_s,
499                             PRIM_INT__p_s_d_s,
500                             Q_PA,
501                             a_over_q,
502                             aoq_PQ,
503                             one_over_2pq,
504                             one_over_2q,
505                             1);
506 
507 
508 
509                     // Forming PRIM_INT__f_s_s_s[5 * 10];
510                     for(n = 0; n < 5; ++n)  // loop over orders of auxiliary function
511                     {
512 
513                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
514                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__f_s_s_s[n * 10 + 0]);
515                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]), PRIM_INT__f_s_s_s[n * 10 + 0]);
516 
517                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
518                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__f_s_s_s[n * 10 + 1]);
519 
520                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
521                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__f_s_s_s[n * 10 + 2]);
522 
523                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
524                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__f_s_s_s[n * 10 + 3]);
525 
526                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
527                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__f_s_s_s[n * 10 + 4]);
528 
529                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
530                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__f_s_s_s[n * 10 + 5]);
531 
532                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
533                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__f_s_s_s[n * 10 + 6]);
534                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_s[n * 3 + 1]), PRIM_INT__f_s_s_s[n * 10 + 6]);
535 
536                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
537                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__f_s_s_s[n * 10 + 7]);
538 
539                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
540                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__f_s_s_s[n * 10 + 8]);
541 
542                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
543                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__f_s_s_s[n * 10 + 9]);
544                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_s[n * 3 + 2]), PRIM_INT__f_s_s_s[n * 10 + 9]);
545 
546                     }
547 
548 
549                     VRR_K_f_s_p_s(
550                             PRIM_INT__f_s_p_s,
551                             PRIM_INT__f_s_s_s,
552                             PRIM_INT__d_s_s_s,
553                             Q_PA,
554                             aoq_PQ,
555                             one_over_2pq,
556                             3);
557 
558 
559                     VRR_K_f_s_d_s(
560                             PRIM_INT__f_s_d_s,
561                             PRIM_INT__f_s_p_s,
562                             PRIM_INT__f_s_s_s,
563                             PRIM_INT__d_s_p_s,
564                             Q_PA,
565                             a_over_q,
566                             aoq_PQ,
567                             one_over_2pq,
568                             one_over_2q,
569                             2);
570 
571 
572                     ostei_general_vrr_K(3, 0, 3, 0, 1,
573                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
574                             PRIM_INT__f_s_d_s, PRIM_INT__f_s_p_s, NULL, PRIM_INT__d_s_d_s, NULL, PRIM_INT__f_s_f_s);
575 
576 
577                     VRR_I_g_s_s_s(
578                             PRIM_INT__g_s_s_s,
579                             PRIM_INT__f_s_s_s,
580                             PRIM_INT__d_s_s_s,
581                             P_PA,
582                             a_over_p,
583                             aop_PQ,
584                             one_over_2p,
585                             4);
586 
587 
588                     VRR_K_g_s_p_s(
589                             PRIM_INT__g_s_p_s,
590                             PRIM_INT__g_s_s_s,
591                             PRIM_INT__f_s_s_s,
592                             Q_PA,
593                             aoq_PQ,
594                             one_over_2pq,
595                             3);
596 
597 
598                     ostei_general_vrr_K(4, 0, 2, 0, 2,
599                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
600                             PRIM_INT__g_s_p_s, PRIM_INT__g_s_s_s, NULL, PRIM_INT__f_s_p_s, NULL, PRIM_INT__g_s_d_s);
601 
602 
603                     ostei_general_vrr_K(4, 0, 3, 0, 1,
604                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
605                             PRIM_INT__g_s_d_s, PRIM_INT__g_s_p_s, NULL, PRIM_INT__f_s_d_s, NULL, PRIM_INT__g_s_f_s);
606 
607 
608 
609 
610                     ////////////////////////////////////
611                     // Accumulate contracted integrals
612                     ////////////////////////////////////
613                     if(lastoffset == 0)
614                     {
615                         contract_all(36, PRIM_INT__d_s_d_s, PRIM_PTR_INT__d_s_d_s);
616                         contract_all(60, PRIM_INT__d_s_f_s, PRIM_PTR_INT__d_s_f_s);
617                         contract_all(60, PRIM_INT__f_s_d_s, PRIM_PTR_INT__f_s_d_s);
618                         contract_all(100, PRIM_INT__f_s_f_s, PRIM_PTR_INT__f_s_f_s);
619                         contract_all(90, PRIM_INT__g_s_d_s, PRIM_PTR_INT__g_s_d_s);
620                         contract_all(150, PRIM_INT__g_s_f_s, PRIM_PTR_INT__g_s_f_s);
621                     }
622                     else
623                     {
624                         contract(36, shelloffsets, PRIM_INT__d_s_d_s, PRIM_PTR_INT__d_s_d_s);
625                         contract(60, shelloffsets, PRIM_INT__d_s_f_s, PRIM_PTR_INT__d_s_f_s);
626                         contract(60, shelloffsets, PRIM_INT__f_s_d_s, PRIM_PTR_INT__f_s_d_s);
627                         contract(100, shelloffsets, PRIM_INT__f_s_f_s, PRIM_PTR_INT__f_s_f_s);
628                         contract(90, shelloffsets, PRIM_INT__g_s_d_s, PRIM_PTR_INT__g_s_d_s);
629                         contract(150, shelloffsets, PRIM_INT__g_s_f_s, PRIM_PTR_INT__g_s_f_s);
630                         PRIM_PTR_INT__d_s_d_s += lastoffset*36;
631                         PRIM_PTR_INT__d_s_f_s += lastoffset*60;
632                         PRIM_PTR_INT__f_s_d_s += lastoffset*60;
633                         PRIM_PTR_INT__f_s_f_s += lastoffset*100;
634                         PRIM_PTR_INT__g_s_d_s += lastoffset*90;
635                         PRIM_PTR_INT__g_s_f_s += lastoffset*150;
636                     }
637 
638                 }  // close loop over j
639             }  // close loop over i
640 
641             //Advance to the next batch
642             jstart = SIMINT_SIMD_ROUND(jend);
643 
644             //////////////////////////////////////////////
645             // Contracted integrals: Horizontal recurrance
646             //////////////////////////////////////////////
647 
648 
649             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
650 
651 
652             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
653             {
654                 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
655 
656                 // set up HRR pointers
657                 double const * restrict HRR_INT__d_s_d_s = INT__d_s_d_s + abcd * 36;
658                 double const * restrict HRR_INT__d_s_f_s = INT__d_s_f_s + abcd * 60;
659                 double const * restrict HRR_INT__f_s_d_s = INT__f_s_d_s + abcd * 60;
660                 double const * restrict HRR_INT__f_s_f_s = INT__f_s_f_s + abcd * 100;
661                 double const * restrict HRR_INT__g_s_d_s = INT__g_s_d_s + abcd * 90;
662                 double const * restrict HRR_INT__g_s_f_s = INT__g_s_f_s + abcd * 150;
663                 double * restrict HRR_INT__d_d_d_p = INT__d_d_d_p + real_abcd * 648;
664 
665                 // form INT__d_p_d_s
666                 for(iket = 0; iket < 6; ++iket)
667                 {
668                     HRR_INT__d_p_d_s[0 * 6 + iket] = HRR_INT__f_s_d_s[0 * 6 + iket] + ( hAB[0] * HRR_INT__d_s_d_s[0 * 6 + iket] );
669 
670                     HRR_INT__d_p_d_s[1 * 6 + iket] = HRR_INT__f_s_d_s[1 * 6 + iket] + ( hAB[1] * HRR_INT__d_s_d_s[0 * 6 + iket] );
671 
672                     HRR_INT__d_p_d_s[2 * 6 + iket] = HRR_INT__f_s_d_s[2 * 6 + iket] + ( hAB[2] * HRR_INT__d_s_d_s[0 * 6 + iket] );
673 
674                     HRR_INT__d_p_d_s[3 * 6 + iket] = HRR_INT__f_s_d_s[1 * 6 + iket] + ( hAB[0] * HRR_INT__d_s_d_s[1 * 6 + iket] );
675 
676                     HRR_INT__d_p_d_s[4 * 6 + iket] = HRR_INT__f_s_d_s[3 * 6 + iket] + ( hAB[1] * HRR_INT__d_s_d_s[1 * 6 + iket] );
677 
678                     HRR_INT__d_p_d_s[5 * 6 + iket] = HRR_INT__f_s_d_s[4 * 6 + iket] + ( hAB[2] * HRR_INT__d_s_d_s[1 * 6 + iket] );
679 
680                     HRR_INT__d_p_d_s[6 * 6 + iket] = HRR_INT__f_s_d_s[2 * 6 + iket] + ( hAB[0] * HRR_INT__d_s_d_s[2 * 6 + iket] );
681 
682                     HRR_INT__d_p_d_s[7 * 6 + iket] = HRR_INT__f_s_d_s[4 * 6 + iket] + ( hAB[1] * HRR_INT__d_s_d_s[2 * 6 + iket] );
683 
684                     HRR_INT__d_p_d_s[8 * 6 + iket] = HRR_INT__f_s_d_s[5 * 6 + iket] + ( hAB[2] * HRR_INT__d_s_d_s[2 * 6 + iket] );
685 
686                     HRR_INT__d_p_d_s[9 * 6 + iket] = HRR_INT__f_s_d_s[3 * 6 + iket] + ( hAB[0] * HRR_INT__d_s_d_s[3 * 6 + iket] );
687 
688                     HRR_INT__d_p_d_s[10 * 6 + iket] = HRR_INT__f_s_d_s[6 * 6 + iket] + ( hAB[1] * HRR_INT__d_s_d_s[3 * 6 + iket] );
689 
690                     HRR_INT__d_p_d_s[11 * 6 + iket] = HRR_INT__f_s_d_s[7 * 6 + iket] + ( hAB[2] * HRR_INT__d_s_d_s[3 * 6 + iket] );
691 
692                     HRR_INT__d_p_d_s[12 * 6 + iket] = HRR_INT__f_s_d_s[4 * 6 + iket] + ( hAB[0] * HRR_INT__d_s_d_s[4 * 6 + iket] );
693 
694                     HRR_INT__d_p_d_s[13 * 6 + iket] = HRR_INT__f_s_d_s[7 * 6 + iket] + ( hAB[1] * HRR_INT__d_s_d_s[4 * 6 + iket] );
695 
696                     HRR_INT__d_p_d_s[14 * 6 + iket] = HRR_INT__f_s_d_s[8 * 6 + iket] + ( hAB[2] * HRR_INT__d_s_d_s[4 * 6 + iket] );
697 
698                     HRR_INT__d_p_d_s[15 * 6 + iket] = HRR_INT__f_s_d_s[5 * 6 + iket] + ( hAB[0] * HRR_INT__d_s_d_s[5 * 6 + iket] );
699 
700                     HRR_INT__d_p_d_s[16 * 6 + iket] = HRR_INT__f_s_d_s[8 * 6 + iket] + ( hAB[1] * HRR_INT__d_s_d_s[5 * 6 + iket] );
701 
702                     HRR_INT__d_p_d_s[17 * 6 + iket] = HRR_INT__f_s_d_s[9 * 6 + iket] + ( hAB[2] * HRR_INT__d_s_d_s[5 * 6 + iket] );
703 
704                 }
705 
706 
707                 // form INT__d_p_f_s
708                 for(iket = 0; iket < 10; ++iket)
709                 {
710                     HRR_INT__d_p_f_s[0 * 10 + iket] = HRR_INT__f_s_f_s[0 * 10 + iket] + ( hAB[0] * HRR_INT__d_s_f_s[0 * 10 + iket] );
711 
712                     HRR_INT__d_p_f_s[1 * 10 + iket] = HRR_INT__f_s_f_s[1 * 10 + iket] + ( hAB[1] * HRR_INT__d_s_f_s[0 * 10 + iket] );
713 
714                     HRR_INT__d_p_f_s[2 * 10 + iket] = HRR_INT__f_s_f_s[2 * 10 + iket] + ( hAB[2] * HRR_INT__d_s_f_s[0 * 10 + iket] );
715 
716                     HRR_INT__d_p_f_s[3 * 10 + iket] = HRR_INT__f_s_f_s[1 * 10 + iket] + ( hAB[0] * HRR_INT__d_s_f_s[1 * 10 + iket] );
717 
718                     HRR_INT__d_p_f_s[4 * 10 + iket] = HRR_INT__f_s_f_s[3 * 10 + iket] + ( hAB[1] * HRR_INT__d_s_f_s[1 * 10 + iket] );
719 
720                     HRR_INT__d_p_f_s[5 * 10 + iket] = HRR_INT__f_s_f_s[4 * 10 + iket] + ( hAB[2] * HRR_INT__d_s_f_s[1 * 10 + iket] );
721 
722                     HRR_INT__d_p_f_s[6 * 10 + iket] = HRR_INT__f_s_f_s[2 * 10 + iket] + ( hAB[0] * HRR_INT__d_s_f_s[2 * 10 + iket] );
723 
724                     HRR_INT__d_p_f_s[7 * 10 + iket] = HRR_INT__f_s_f_s[4 * 10 + iket] + ( hAB[1] * HRR_INT__d_s_f_s[2 * 10 + iket] );
725 
726                     HRR_INT__d_p_f_s[8 * 10 + iket] = HRR_INT__f_s_f_s[5 * 10 + iket] + ( hAB[2] * HRR_INT__d_s_f_s[2 * 10 + iket] );
727 
728                     HRR_INT__d_p_f_s[9 * 10 + iket] = HRR_INT__f_s_f_s[3 * 10 + iket] + ( hAB[0] * HRR_INT__d_s_f_s[3 * 10 + iket] );
729 
730                     HRR_INT__d_p_f_s[10 * 10 + iket] = HRR_INT__f_s_f_s[6 * 10 + iket] + ( hAB[1] * HRR_INT__d_s_f_s[3 * 10 + iket] );
731 
732                     HRR_INT__d_p_f_s[11 * 10 + iket] = HRR_INT__f_s_f_s[7 * 10 + iket] + ( hAB[2] * HRR_INT__d_s_f_s[3 * 10 + iket] );
733 
734                     HRR_INT__d_p_f_s[12 * 10 + iket] = HRR_INT__f_s_f_s[4 * 10 + iket] + ( hAB[0] * HRR_INT__d_s_f_s[4 * 10 + iket] );
735 
736                     HRR_INT__d_p_f_s[13 * 10 + iket] = HRR_INT__f_s_f_s[7 * 10 + iket] + ( hAB[1] * HRR_INT__d_s_f_s[4 * 10 + iket] );
737 
738                     HRR_INT__d_p_f_s[14 * 10 + iket] = HRR_INT__f_s_f_s[8 * 10 + iket] + ( hAB[2] * HRR_INT__d_s_f_s[4 * 10 + iket] );
739 
740                     HRR_INT__d_p_f_s[15 * 10 + iket] = HRR_INT__f_s_f_s[5 * 10 + iket] + ( hAB[0] * HRR_INT__d_s_f_s[5 * 10 + iket] );
741 
742                     HRR_INT__d_p_f_s[16 * 10 + iket] = HRR_INT__f_s_f_s[8 * 10 + iket] + ( hAB[1] * HRR_INT__d_s_f_s[5 * 10 + iket] );
743 
744                     HRR_INT__d_p_f_s[17 * 10 + iket] = HRR_INT__f_s_f_s[9 * 10 + iket] + ( hAB[2] * HRR_INT__d_s_f_s[5 * 10 + iket] );
745 
746                 }
747 
748 
749                 // form INT__f_p_d_s
750                 HRR_J_f_p(
751                     HRR_INT__f_p_d_s,
752                     HRR_INT__f_s_d_s,
753                     HRR_INT__g_s_d_s,
754                     hAB, 6);
755 
756                 // form INT__f_p_f_s
757                 HRR_J_f_p(
758                     HRR_INT__f_p_f_s,
759                     HRR_INT__f_s_f_s,
760                     HRR_INT__g_s_f_s,
761                     hAB, 10);
762 
763                 // form INT__d_d_d_s
764                 HRR_J_d_d(
765                     HRR_INT__d_d_d_s,
766                     HRR_INT__d_p_d_s,
767                     HRR_INT__f_p_d_s,
768                     hAB, 6);
769 
770                 // form INT__d_d_f_s
771                 HRR_J_d_d(
772                     HRR_INT__d_d_f_s,
773                     HRR_INT__d_p_f_s,
774                     HRR_INT__f_p_f_s,
775                     hAB, 10);
776 
777                 // form INT__d_d_d_p
778                 for(ibra = 0; ibra < 36; ++ibra)
779                 {
780                     HRR_INT__d_d_d_p[ibra * 18 + 0] = HRR_INT__d_d_f_s[ibra * 10 + 0] + ( hCD[0] * HRR_INT__d_d_d_s[ibra * 6 + 0] );
781 
782                     HRR_INT__d_d_d_p[ibra * 18 + 1] = HRR_INT__d_d_f_s[ibra * 10 + 1] + ( hCD[1] * HRR_INT__d_d_d_s[ibra * 6 + 0] );
783 
784                     HRR_INT__d_d_d_p[ibra * 18 + 2] = HRR_INT__d_d_f_s[ibra * 10 + 2] + ( hCD[2] * HRR_INT__d_d_d_s[ibra * 6 + 0] );
785 
786                     HRR_INT__d_d_d_p[ibra * 18 + 3] = HRR_INT__d_d_f_s[ibra * 10 + 1] + ( hCD[0] * HRR_INT__d_d_d_s[ibra * 6 + 1] );
787 
788                     HRR_INT__d_d_d_p[ibra * 18 + 4] = HRR_INT__d_d_f_s[ibra * 10 + 3] + ( hCD[1] * HRR_INT__d_d_d_s[ibra * 6 + 1] );
789 
790                     HRR_INT__d_d_d_p[ibra * 18 + 5] = HRR_INT__d_d_f_s[ibra * 10 + 4] + ( hCD[2] * HRR_INT__d_d_d_s[ibra * 6 + 1] );
791 
792                     HRR_INT__d_d_d_p[ibra * 18 + 6] = HRR_INT__d_d_f_s[ibra * 10 + 2] + ( hCD[0] * HRR_INT__d_d_d_s[ibra * 6 + 2] );
793 
794                     HRR_INT__d_d_d_p[ibra * 18 + 7] = HRR_INT__d_d_f_s[ibra * 10 + 4] + ( hCD[1] * HRR_INT__d_d_d_s[ibra * 6 + 2] );
795 
796                     HRR_INT__d_d_d_p[ibra * 18 + 8] = HRR_INT__d_d_f_s[ibra * 10 + 5] + ( hCD[2] * HRR_INT__d_d_d_s[ibra * 6 + 2] );
797 
798                     HRR_INT__d_d_d_p[ibra * 18 + 9] = HRR_INT__d_d_f_s[ibra * 10 + 3] + ( hCD[0] * HRR_INT__d_d_d_s[ibra * 6 + 3] );
799 
800                     HRR_INT__d_d_d_p[ibra * 18 + 10] = HRR_INT__d_d_f_s[ibra * 10 + 6] + ( hCD[1] * HRR_INT__d_d_d_s[ibra * 6 + 3] );
801 
802                     HRR_INT__d_d_d_p[ibra * 18 + 11] = HRR_INT__d_d_f_s[ibra * 10 + 7] + ( hCD[2] * HRR_INT__d_d_d_s[ibra * 6 + 3] );
803 
804                     HRR_INT__d_d_d_p[ibra * 18 + 12] = HRR_INT__d_d_f_s[ibra * 10 + 4] + ( hCD[0] * HRR_INT__d_d_d_s[ibra * 6 + 4] );
805 
806                     HRR_INT__d_d_d_p[ibra * 18 + 13] = HRR_INT__d_d_f_s[ibra * 10 + 7] + ( hCD[1] * HRR_INT__d_d_d_s[ibra * 6 + 4] );
807 
808                     HRR_INT__d_d_d_p[ibra * 18 + 14] = HRR_INT__d_d_f_s[ibra * 10 + 8] + ( hCD[2] * HRR_INT__d_d_d_s[ibra * 6 + 4] );
809 
810                     HRR_INT__d_d_d_p[ibra * 18 + 15] = HRR_INT__d_d_f_s[ibra * 10 + 5] + ( hCD[0] * HRR_INT__d_d_d_s[ibra * 6 + 5] );
811 
812                     HRR_INT__d_d_d_p[ibra * 18 + 16] = HRR_INT__d_d_f_s[ibra * 10 + 8] + ( hCD[1] * HRR_INT__d_d_d_s[ibra * 6 + 5] );
813 
814                     HRR_INT__d_d_d_p[ibra * 18 + 17] = HRR_INT__d_d_f_s[ibra * 10 + 9] + ( hCD[2] * HRR_INT__d_d_d_s[ibra * 6 + 5] );
815 
816                 }
817 
818 
819             }  // close HRR loop
820 
821 
822         }   // close loop cdbatch
823 
824         istart = iend;
825     }  // close loop over ab
826 
827     return P.nshell12_clip * Q.nshell12_clip;
828 }
829 
830