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