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