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_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__d_f_d_d)8 int ostei_d_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__d_f_d_d)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__d_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__s_f_d_s = work + (SIMINT_NSHELL_SIMD * 0);
30     double * const INT__s_f_f_s = work + (SIMINT_NSHELL_SIMD * 60);
31     double * const INT__s_f_g_s = work + (SIMINT_NSHELL_SIMD * 160);
32     double * const INT__s_g_d_s = work + (SIMINT_NSHELL_SIMD * 310);
33     double * const INT__s_g_f_s = work + (SIMINT_NSHELL_SIMD * 400);
34     double * const INT__s_g_g_s = work + (SIMINT_NSHELL_SIMD * 550);
35     double * const INT__s_h_d_s = work + (SIMINT_NSHELL_SIMD * 775);
36     double * const INT__s_h_f_s = work + (SIMINT_NSHELL_SIMD * 901);
37     double * const INT__s_h_g_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_p_s_s = primwork + 22;
42     SIMINT_DBLTYPE * const restrict PRIM_INT__s_p_p_s = primwork + 49;
43     SIMINT_DBLTYPE * const restrict PRIM_INT__s_p_d_s = primwork + 85;
44     SIMINT_DBLTYPE * const restrict PRIM_INT__s_d_s_s = primwork + 139;
45     SIMINT_DBLTYPE * const restrict PRIM_INT__s_d_p_s = primwork + 187;
46     SIMINT_DBLTYPE * const restrict PRIM_INT__s_d_d_s = primwork + 259;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__s_d_f_s = primwork + 367;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__s_f_s_s = primwork + 487;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__s_f_p_s = primwork + 557;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__s_f_d_s = primwork + 677;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__s_f_f_s = primwork + 857;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__s_f_g_s = primwork + 1057;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__s_g_s_s = primwork + 1207;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__s_g_p_s = primwork + 1297;
55     SIMINT_DBLTYPE * const restrict PRIM_INT__s_g_d_s = primwork + 1477;
56     SIMINT_DBLTYPE * const restrict PRIM_INT__s_g_f_s = primwork + 1747;
57     SIMINT_DBLTYPE * const restrict PRIM_INT__s_g_g_s = primwork + 2047;
58     SIMINT_DBLTYPE * const restrict PRIM_INT__s_h_s_s = primwork + 2272;
59     SIMINT_DBLTYPE * const restrict PRIM_INT__s_h_p_s = primwork + 2377;
60     SIMINT_DBLTYPE * const restrict PRIM_INT__s_h_d_s = primwork + 2629;
61     SIMINT_DBLTYPE * const restrict PRIM_INT__s_h_f_s = primwork + 3007;
62     SIMINT_DBLTYPE * const restrict PRIM_INT__s_h_g_s = primwork + 3427;
63     double * const hrrwork = (double *)(primwork + 3742);
64     double * const HRR_INT__p_f_d_s = hrrwork + 0;
65     double * const HRR_INT__p_f_f_s = hrrwork + 180;
66     double * const HRR_INT__p_f_g_s = hrrwork + 480;
67     double * const HRR_INT__p_g_d_s = hrrwork + 930;
68     double * const HRR_INT__p_g_f_s = hrrwork + 1200;
69     double * const HRR_INT__p_g_g_s = hrrwork + 1650;
70     double * const HRR_INT__d_f_d_s = hrrwork + 2325;
71     double * const HRR_INT__d_f_d_p = hrrwork + 2685;
72     double * const HRR_INT__d_f_f_s = hrrwork + 3765;
73     double * const HRR_INT__d_f_f_p = hrrwork + 4365;
74     double * const HRR_INT__d_f_g_s = hrrwork + 6165;
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__s_f_d_s = INT__s_f_d_s + abcd * 60;
127                 double * restrict PRIM_PTR_INT__s_f_f_s = INT__s_f_f_s + abcd * 100;
128                 double * restrict PRIM_PTR_INT__s_f_g_s = INT__s_f_g_s + abcd * 150;
129                 double * restrict PRIM_PTR_INT__s_g_d_s = INT__s_g_d_s + abcd * 90;
130                 double * restrict PRIM_PTR_INT__s_g_f_s = INT__s_g_f_s + abcd * 150;
131                 double * restrict PRIM_PTR_INT__s_g_g_s = INT__s_g_g_s + abcd * 225;
132                 double * restrict PRIM_PTR_INT__s_h_d_s = INT__s_h_d_s + abcd * 126;
133                 double * restrict PRIM_PTR_INT__s_h_f_s = INT__s_h_f_s + abcd * 210;
134                 double * restrict PRIM_PTR_INT__s_h_g_s = INT__s_h_g_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_PB[3] = { SIMINT_DBLSET1(P.PB_x[i]), SIMINT_DBLSET1(P.PB_y[i]), SIMINT_DBLSET1(P.PB_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__s_f_d_s += 60;
161                             PRIM_PTR_INT__s_f_f_s += 100;
162                             PRIM_PTR_INT__s_f_g_s += 150;
163                             PRIM_PTR_INT__s_g_d_s += 90;
164                             PRIM_PTR_INT__s_g_f_s += 150;
165                             PRIM_PTR_INT__s_g_g_s += 225;
166                             PRIM_PTR_INT__s_h_d_s += 126;
167                             PRIM_PTR_INT__s_h_f_s += 210;
168                             PRIM_PTR_INT__s_h_g_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__s_f_d_s += lastoffset*60;
195                             PRIM_PTR_INT__s_f_f_s += lastoffset*100;
196                             PRIM_PTR_INT__s_f_g_s += lastoffset*150;
197                             PRIM_PTR_INT__s_g_d_s += lastoffset*90;
198                             PRIM_PTR_INT__s_g_f_s += lastoffset*150;
199                             PRIM_PTR_INT__s_g_g_s += lastoffset*225;
200                             PRIM_PTR_INT__s_h_d_s += lastoffset*126;
201                             PRIM_PTR_INT__s_h_f_s += lastoffset*210;
202                             PRIM_PTR_INT__s_h_g_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_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
271                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
272                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
273                     const SIMINT_DBLTYPE vrr_const_3_over_2q = SIMINT_MUL(const_3, 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_p_s_s[9 * 3];
283                     for(n = 0; n < 9; ++n)  // loop over orders of auxiliary function
284                     {
285 
286                         PRIM_INT__s_p_s_s[n * 3 + 0] = SIMINT_MUL(P_PB[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
287                         PRIM_INT__s_p_s_s[n * 3 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_p_s_s[n * 3 + 0]);
288 
289                         PRIM_INT__s_p_s_s[n * 3 + 1] = SIMINT_MUL(P_PB[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
290                         PRIM_INT__s_p_s_s[n * 3 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_p_s_s[n * 3 + 1]);
291 
292                         PRIM_INT__s_p_s_s[n * 3 + 2] = SIMINT_MUL(P_PB[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
293                         PRIM_INT__s_p_s_s[n * 3 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_p_s_s[n * 3 + 2]);
294 
295                     }
296 
297 
298 
299                     // Forming PRIM_INT__s_d_s_s[8 * 6];
300                     for(n = 0; n < 8; ++n)  // loop over orders of auxiliary function
301                     {
302 
303                         PRIM_INT__s_d_s_s[n * 6 + 0] = SIMINT_MUL(P_PB[0], PRIM_INT__s_p_s_s[n * 3 + 0]);
304                         PRIM_INT__s_d_s_s[n * 6 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_d_s_s[n * 6 + 0]);
305                         PRIM_INT__s_d_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__s_d_s_s[n * 6 + 0]);
306 
307                         PRIM_INT__s_d_s_s[n * 6 + 1] = SIMINT_MUL(P_PB[1], PRIM_INT__s_p_s_s[n * 3 + 0]);
308                         PRIM_INT__s_d_s_s[n * 6 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_d_s_s[n * 6 + 1]);
309 
310                         PRIM_INT__s_d_s_s[n * 6 + 2] = SIMINT_MUL(P_PB[2], PRIM_INT__s_p_s_s[n * 3 + 0]);
311                         PRIM_INT__s_d_s_s[n * 6 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_d_s_s[n * 6 + 2]);
312 
313                         PRIM_INT__s_d_s_s[n * 6 + 3] = SIMINT_MUL(P_PB[1], PRIM_INT__s_p_s_s[n * 3 + 1]);
314                         PRIM_INT__s_d_s_s[n * 6 + 3] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_p_s_s[(n+1) * 3 + 1], PRIM_INT__s_d_s_s[n * 6 + 3]);
315                         PRIM_INT__s_d_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__s_d_s_s[n * 6 + 3]);
316 
317                         PRIM_INT__s_d_s_s[n * 6 + 4] = SIMINT_MUL(P_PB[2], PRIM_INT__s_p_s_s[n * 3 + 1]);
318                         PRIM_INT__s_d_s_s[n * 6 + 4] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_p_s_s[(n+1) * 3 + 1], PRIM_INT__s_d_s_s[n * 6 + 4]);
319 
320                         PRIM_INT__s_d_s_s[n * 6 + 5] = SIMINT_MUL(P_PB[2], PRIM_INT__s_p_s_s[n * 3 + 2]);
321                         PRIM_INT__s_d_s_s[n * 6 + 5] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_p_s_s[(n+1) * 3 + 2], PRIM_INT__s_d_s_s[n * 6 + 5]);
322                         PRIM_INT__s_d_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__s_d_s_s[n * 6 + 5]);
323 
324                     }
325 
326 
327 
328                     // Forming PRIM_INT__s_f_s_s[7 * 10];
329                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
330                     {
331 
332                         PRIM_INT__s_f_s_s[n * 10 + 0] = SIMINT_MUL(P_PB[0], PRIM_INT__s_d_s_s[n * 6 + 0]);
333                         PRIM_INT__s_f_s_s[n * 10 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_d_s_s[(n+1) * 6 + 0], PRIM_INT__s_f_s_s[n * 10 + 0]);
334                         PRIM_INT__s_f_s_s[n * 10 + 0] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_p_s_s[n * 3 + 0]), PRIM_INT__s_f_s_s[n * 10 + 0]);
335 
336                         PRIM_INT__s_f_s_s[n * 10 + 1] = SIMINT_MUL(P_PB[1], PRIM_INT__s_d_s_s[n * 6 + 0]);
337                         PRIM_INT__s_f_s_s[n * 10 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_d_s_s[(n+1) * 6 + 0], PRIM_INT__s_f_s_s[n * 10 + 1]);
338 
339                         PRIM_INT__s_f_s_s[n * 10 + 2] = SIMINT_MUL(P_PB[2], PRIM_INT__s_d_s_s[n * 6 + 0]);
340                         PRIM_INT__s_f_s_s[n * 10 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_d_s_s[(n+1) * 6 + 0], PRIM_INT__s_f_s_s[n * 10 + 2]);
341 
342                         PRIM_INT__s_f_s_s[n * 10 + 3] = SIMINT_MUL(P_PB[0], PRIM_INT__s_d_s_s[n * 6 + 3]);
343                         PRIM_INT__s_f_s_s[n * 10 + 3] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_d_s_s[(n+1) * 6 + 3], PRIM_INT__s_f_s_s[n * 10 + 3]);
344 
345                         PRIM_INT__s_f_s_s[n * 10 + 4] = SIMINT_MUL(P_PB[2], PRIM_INT__s_d_s_s[n * 6 + 1]);
346                         PRIM_INT__s_f_s_s[n * 10 + 4] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_d_s_s[(n+1) * 6 + 1], PRIM_INT__s_f_s_s[n * 10 + 4]);
347 
348                         PRIM_INT__s_f_s_s[n * 10 + 5] = SIMINT_MUL(P_PB[0], PRIM_INT__s_d_s_s[n * 6 + 5]);
349                         PRIM_INT__s_f_s_s[n * 10 + 5] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_d_s_s[(n+1) * 6 + 5], PRIM_INT__s_f_s_s[n * 10 + 5]);
350 
351                         PRIM_INT__s_f_s_s[n * 10 + 6] = SIMINT_MUL(P_PB[1], PRIM_INT__s_d_s_s[n * 6 + 3]);
352                         PRIM_INT__s_f_s_s[n * 10 + 6] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_d_s_s[(n+1) * 6 + 3], PRIM_INT__s_f_s_s[n * 10 + 6]);
353                         PRIM_INT__s_f_s_s[n * 10 + 6] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_p_s_s[(n+1) * 3 + 1], PRIM_INT__s_p_s_s[n * 3 + 1]), PRIM_INT__s_f_s_s[n * 10 + 6]);
354 
355                         PRIM_INT__s_f_s_s[n * 10 + 7] = SIMINT_MUL(P_PB[2], PRIM_INT__s_d_s_s[n * 6 + 3]);
356                         PRIM_INT__s_f_s_s[n * 10 + 7] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_d_s_s[(n+1) * 6 + 3], PRIM_INT__s_f_s_s[n * 10 + 7]);
357 
358                         PRIM_INT__s_f_s_s[n * 10 + 8] = SIMINT_MUL(P_PB[1], PRIM_INT__s_d_s_s[n * 6 + 5]);
359                         PRIM_INT__s_f_s_s[n * 10 + 8] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_d_s_s[(n+1) * 6 + 5], PRIM_INT__s_f_s_s[n * 10 + 8]);
360 
361                         PRIM_INT__s_f_s_s[n * 10 + 9] = SIMINT_MUL(P_PB[2], PRIM_INT__s_d_s_s[n * 6 + 5]);
362                         PRIM_INT__s_f_s_s[n * 10 + 9] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_d_s_s[(n+1) * 6 + 5], PRIM_INT__s_f_s_s[n * 10 + 9]);
363                         PRIM_INT__s_f_s_s[n * 10 + 9] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_p_s_s[(n+1) * 3 + 2], PRIM_INT__s_p_s_s[n * 3 + 2]), PRIM_INT__s_f_s_s[n * 10 + 9]);
364 
365                     }
366 
367 
368                     VRR_K_s_f_p_s(
369                             PRIM_INT__s_f_p_s,
370                             PRIM_INT__s_f_s_s,
371                             PRIM_INT__s_d_s_s,
372                             Q_PA,
373                             aoq_PQ,
374                             one_over_2pq,
375                             4);
376 
377 
378 
379                     // Forming PRIM_INT__s_d_p_s[4 * 18];
380                     for(n = 0; n < 4; ++n)  // loop over orders of auxiliary function
381                     {
382 
383                         PRIM_INT__s_d_p_s[n * 18 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_d_s_s[n * 6 + 0]);
384                         PRIM_INT__s_d_p_s[n * 18 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_d_s_s[(n+1) * 6 + 0], PRIM_INT__s_d_p_s[n * 18 + 0]);
385                         PRIM_INT__s_d_p_s[n * 18 + 0] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_d_p_s[n * 18 + 0]);
386 
387                         PRIM_INT__s_d_p_s[n * 18 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_d_s_s[n * 6 + 0]);
388                         PRIM_INT__s_d_p_s[n * 18 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_d_s_s[(n+1) * 6 + 0], PRIM_INT__s_d_p_s[n * 18 + 1]);
389 
390                         PRIM_INT__s_d_p_s[n * 18 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_d_s_s[n * 6 + 0]);
391                         PRIM_INT__s_d_p_s[n * 18 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_d_s_s[(n+1) * 6 + 0], PRIM_INT__s_d_p_s[n * 18 + 2]);
392 
393                         PRIM_INT__s_d_p_s[n * 18 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_d_s_s[n * 6 + 1]);
394                         PRIM_INT__s_d_p_s[n * 18 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_d_s_s[(n+1) * 6 + 1], PRIM_INT__s_d_p_s[n * 18 + 3]);
395                         PRIM_INT__s_d_p_s[n * 18 + 3] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_p_s_s[(n+1) * 3 + 1], PRIM_INT__s_d_p_s[n * 18 + 3]);
396 
397                         PRIM_INT__s_d_p_s[n * 18 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_d_s_s[n * 6 + 1]);
398                         PRIM_INT__s_d_p_s[n * 18 + 4] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_d_s_s[(n+1) * 6 + 1], PRIM_INT__s_d_p_s[n * 18 + 4]);
399                         PRIM_INT__s_d_p_s[n * 18 + 4] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_d_p_s[n * 18 + 4]);
400 
401                         PRIM_INT__s_d_p_s[n * 18 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_d_s_s[n * 6 + 1]);
402                         PRIM_INT__s_d_p_s[n * 18 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_d_s_s[(n+1) * 6 + 1], PRIM_INT__s_d_p_s[n * 18 + 5]);
403 
404                         PRIM_INT__s_d_p_s[n * 18 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_d_s_s[n * 6 + 2]);
405                         PRIM_INT__s_d_p_s[n * 18 + 6] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_d_s_s[(n+1) * 6 + 2], PRIM_INT__s_d_p_s[n * 18 + 6]);
406                         PRIM_INT__s_d_p_s[n * 18 + 6] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_p_s_s[(n+1) * 3 + 2], PRIM_INT__s_d_p_s[n * 18 + 6]);
407 
408                         PRIM_INT__s_d_p_s[n * 18 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_d_s_s[n * 6 + 2]);
409                         PRIM_INT__s_d_p_s[n * 18 + 7] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_d_s_s[(n+1) * 6 + 2], PRIM_INT__s_d_p_s[n * 18 + 7]);
410 
411                         PRIM_INT__s_d_p_s[n * 18 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_d_s_s[n * 6 + 2]);
412                         PRIM_INT__s_d_p_s[n * 18 + 8] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_d_s_s[(n+1) * 6 + 2], PRIM_INT__s_d_p_s[n * 18 + 8]);
413                         PRIM_INT__s_d_p_s[n * 18 + 8] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_d_p_s[n * 18 + 8]);
414 
415                         PRIM_INT__s_d_p_s[n * 18 + 9] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_d_s_s[n * 6 + 3]);
416                         PRIM_INT__s_d_p_s[n * 18 + 9] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_d_s_s[(n+1) * 6 + 3], PRIM_INT__s_d_p_s[n * 18 + 9]);
417 
418                         PRIM_INT__s_d_p_s[n * 18 + 10] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_d_s_s[n * 6 + 3]);
419                         PRIM_INT__s_d_p_s[n * 18 + 10] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_d_s_s[(n+1) * 6 + 3], PRIM_INT__s_d_p_s[n * 18 + 10]);
420                         PRIM_INT__s_d_p_s[n * 18 + 10] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_p_s_s[(n+1) * 3 + 1], PRIM_INT__s_d_p_s[n * 18 + 10]);
421 
422                         PRIM_INT__s_d_p_s[n * 18 + 11] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_d_s_s[n * 6 + 3]);
423                         PRIM_INT__s_d_p_s[n * 18 + 11] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_d_s_s[(n+1) * 6 + 3], PRIM_INT__s_d_p_s[n * 18 + 11]);
424 
425                         PRIM_INT__s_d_p_s[n * 18 + 12] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_d_s_s[n * 6 + 4]);
426                         PRIM_INT__s_d_p_s[n * 18 + 12] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_d_s_s[(n+1) * 6 + 4], PRIM_INT__s_d_p_s[n * 18 + 12]);
427 
428                         PRIM_INT__s_d_p_s[n * 18 + 13] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_d_s_s[n * 6 + 4]);
429                         PRIM_INT__s_d_p_s[n * 18 + 13] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_d_s_s[(n+1) * 6 + 4], PRIM_INT__s_d_p_s[n * 18 + 13]);
430                         PRIM_INT__s_d_p_s[n * 18 + 13] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_p_s_s[(n+1) * 3 + 2], PRIM_INT__s_d_p_s[n * 18 + 13]);
431 
432                         PRIM_INT__s_d_p_s[n * 18 + 14] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_d_s_s[n * 6 + 4]);
433                         PRIM_INT__s_d_p_s[n * 18 + 14] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_d_s_s[(n+1) * 6 + 4], PRIM_INT__s_d_p_s[n * 18 + 14]);
434                         PRIM_INT__s_d_p_s[n * 18 + 14] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_p_s_s[(n+1) * 3 + 1], PRIM_INT__s_d_p_s[n * 18 + 14]);
435 
436                         PRIM_INT__s_d_p_s[n * 18 + 15] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_d_s_s[n * 6 + 5]);
437                         PRIM_INT__s_d_p_s[n * 18 + 15] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_d_s_s[(n+1) * 6 + 5], PRIM_INT__s_d_p_s[n * 18 + 15]);
438 
439                         PRIM_INT__s_d_p_s[n * 18 + 16] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_d_s_s[n * 6 + 5]);
440                         PRIM_INT__s_d_p_s[n * 18 + 16] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_d_s_s[(n+1) * 6 + 5], PRIM_INT__s_d_p_s[n * 18 + 16]);
441 
442                         PRIM_INT__s_d_p_s[n * 18 + 17] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_d_s_s[n * 6 + 5]);
443                         PRIM_INT__s_d_p_s[n * 18 + 17] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_d_s_s[(n+1) * 6 + 5], PRIM_INT__s_d_p_s[n * 18 + 17]);
444                         PRIM_INT__s_d_p_s[n * 18 + 17] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_p_s_s[(n+1) * 3 + 2], PRIM_INT__s_d_p_s[n * 18 + 17]);
445 
446                     }
447 
448 
449                     VRR_K_s_f_d_s(
450                             PRIM_INT__s_f_d_s,
451                             PRIM_INT__s_f_p_s,
452                             PRIM_INT__s_f_s_s,
453                             PRIM_INT__s_d_p_s,
454                             Q_PA,
455                             a_over_q,
456                             aoq_PQ,
457                             one_over_2pq,
458                             one_over_2q,
459                             3);
460 
461 
462 
463                     // Forming PRIM_INT__s_p_p_s[4 * 9];
464                     for(n = 0; n < 4; ++n)  // loop over orders of auxiliary function
465                     {
466 
467                         PRIM_INT__s_p_p_s[n * 9 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_p_s_s[n * 3 + 0]);
468                         PRIM_INT__s_p_p_s[n * 9 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_p_p_s[n * 9 + 0]);
469                         PRIM_INT__s_p_p_s[n * 9 + 0] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_p_p_s[n * 9 + 0]);
470 
471                         PRIM_INT__s_p_p_s[n * 9 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_p_s_s[n * 3 + 0]);
472                         PRIM_INT__s_p_p_s[n * 9 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_p_p_s[n * 9 + 1]);
473 
474                         PRIM_INT__s_p_p_s[n * 9 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_p_s_s[n * 3 + 0]);
475                         PRIM_INT__s_p_p_s[n * 9 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_p_p_s[n * 9 + 2]);
476 
477                         PRIM_INT__s_p_p_s[n * 9 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_p_s_s[n * 3 + 1]);
478                         PRIM_INT__s_p_p_s[n * 9 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_p_s_s[(n+1) * 3 + 1], PRIM_INT__s_p_p_s[n * 9 + 3]);
479 
480                         PRIM_INT__s_p_p_s[n * 9 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_p_s_s[n * 3 + 1]);
481                         PRIM_INT__s_p_p_s[n * 9 + 4] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_p_s_s[(n+1) * 3 + 1], PRIM_INT__s_p_p_s[n * 9 + 4]);
482                         PRIM_INT__s_p_p_s[n * 9 + 4] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_p_p_s[n * 9 + 4]);
483 
484                         PRIM_INT__s_p_p_s[n * 9 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_p_s_s[n * 3 + 1]);
485                         PRIM_INT__s_p_p_s[n * 9 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_p_s_s[(n+1) * 3 + 1], PRIM_INT__s_p_p_s[n * 9 + 5]);
486 
487                         PRIM_INT__s_p_p_s[n * 9 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_p_s_s[n * 3 + 2]);
488                         PRIM_INT__s_p_p_s[n * 9 + 6] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_p_s_s[(n+1) * 3 + 2], PRIM_INT__s_p_p_s[n * 9 + 6]);
489 
490                         PRIM_INT__s_p_p_s[n * 9 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_p_s_s[n * 3 + 2]);
491                         PRIM_INT__s_p_p_s[n * 9 + 7] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_p_s_s[(n+1) * 3 + 2], PRIM_INT__s_p_p_s[n * 9 + 7]);
492 
493                         PRIM_INT__s_p_p_s[n * 9 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_p_s_s[n * 3 + 2]);
494                         PRIM_INT__s_p_p_s[n * 9 + 8] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_p_s_s[(n+1) * 3 + 2], PRIM_INT__s_p_p_s[n * 9 + 8]);
495                         PRIM_INT__s_p_p_s[n * 9 + 8] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_p_p_s[n * 9 + 8]);
496 
497                     }
498 
499 
500                     VRR_K_s_d_d_s(
501                             PRIM_INT__s_d_d_s,
502                             PRIM_INT__s_d_p_s,
503                             PRIM_INT__s_d_s_s,
504                             PRIM_INT__s_p_p_s,
505                             Q_PA,
506                             a_over_q,
507                             aoq_PQ,
508                             one_over_2pq,
509                             one_over_2q,
510                             3);
511 
512 
513                     ostei_general_vrr_K(0, 3, 3, 0, 2,
514                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
515                             PRIM_INT__s_f_d_s, PRIM_INT__s_f_p_s, NULL, NULL, PRIM_INT__s_d_d_s, PRIM_INT__s_f_f_s);
516 
517 
518                     VRR_J_s_g_s_s(
519                             PRIM_INT__s_g_s_s,
520                             PRIM_INT__s_f_s_s,
521                             PRIM_INT__s_d_s_s,
522                             P_PB,
523                             a_over_p,
524                             aop_PQ,
525                             one_over_2p,
526                             6);
527 
528 
529                     VRR_K_s_g_p_s(
530                             PRIM_INT__s_g_p_s,
531                             PRIM_INT__s_g_s_s,
532                             PRIM_INT__s_f_s_s,
533                             Q_PA,
534                             aoq_PQ,
535                             one_over_2pq,
536                             4);
537 
538 
539                     ostei_general_vrr_K(0, 4, 2, 0, 3,
540                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
541                             PRIM_INT__s_g_p_s, PRIM_INT__s_g_s_s, NULL, NULL, PRIM_INT__s_f_p_s, PRIM_INT__s_g_d_s);
542 
543 
544 
545                     // Forming PRIM_INT__s_s_p_s[4 * 3];
546                     for(n = 0; n < 4; ++n)  // loop over orders of auxiliary function
547                     {
548 
549                         PRIM_INT__s_s_p_s[n * 3 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
550                         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]);
551 
552                         PRIM_INT__s_s_p_s[n * 3 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
553                         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]);
554 
555                         PRIM_INT__s_s_p_s[n * 3 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
556                         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]);
557 
558                     }
559 
560 
561 
562                     // Forming PRIM_INT__s_p_d_s[3 * 18];
563                     for(n = 0; n < 3; ++n)  // loop over orders of auxiliary function
564                     {
565 
566                         PRIM_INT__s_p_d_s[n * 18 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_p_p_s[n * 9 + 0]);
567                         PRIM_INT__s_p_d_s[n * 18 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_p_p_s[(n+1) * 9 + 0], PRIM_INT__s_p_d_s[n * 18 + 0]);
568                         PRIM_INT__s_p_d_s[n * 18 + 0] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_p_s_s[n * 3 + 0]), PRIM_INT__s_p_d_s[n * 18 + 0]);
569                         PRIM_INT__s_p_d_s[n * 18 + 0] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_p_d_s[n * 18 + 0]);
570 
571                         PRIM_INT__s_p_d_s[n * 18 + 3] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_p_p_s[n * 9 + 1]);
572                         PRIM_INT__s_p_d_s[n * 18 + 3] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_p_p_s[(n+1) * 9 + 1], PRIM_INT__s_p_d_s[n * 18 + 3]);
573                         PRIM_INT__s_p_d_s[n * 18 + 3] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_p_s_s[n * 3 + 0]), PRIM_INT__s_p_d_s[n * 18 + 3]);
574 
575                         PRIM_INT__s_p_d_s[n * 18 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_p_p_s[n * 9 + 2]);
576                         PRIM_INT__s_p_d_s[n * 18 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_p_p_s[(n+1) * 9 + 2], PRIM_INT__s_p_d_s[n * 18 + 5]);
577                         PRIM_INT__s_p_d_s[n * 18 + 5] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_p_s_s[n * 3 + 0]), PRIM_INT__s_p_d_s[n * 18 + 5]);
578 
579                         PRIM_INT__s_p_d_s[n * 18 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_p_p_s[n * 9 + 3]);
580                         PRIM_INT__s_p_d_s[n * 18 + 6] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_p_p_s[(n+1) * 9 + 3], PRIM_INT__s_p_d_s[n * 18 + 6]);
581                         PRIM_INT__s_p_d_s[n * 18 + 6] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_p_s_s[(n+1) * 3 + 1], PRIM_INT__s_p_s_s[n * 3 + 1]), PRIM_INT__s_p_d_s[n * 18 + 6]);
582 
583                         PRIM_INT__s_p_d_s[n * 18 + 9] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_p_p_s[n * 9 + 4]);
584                         PRIM_INT__s_p_d_s[n * 18 + 9] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_p_p_s[(n+1) * 9 + 4], PRIM_INT__s_p_d_s[n * 18 + 9]);
585                         PRIM_INT__s_p_d_s[n * 18 + 9] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_p_s_s[(n+1) * 3 + 1], PRIM_INT__s_p_s_s[n * 3 + 1]), PRIM_INT__s_p_d_s[n * 18 + 9]);
586                         PRIM_INT__s_p_d_s[n * 18 + 9] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_p_d_s[n * 18 + 9]);
587 
588                         PRIM_INT__s_p_d_s[n * 18 + 11] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_p_p_s[n * 9 + 5]);
589                         PRIM_INT__s_p_d_s[n * 18 + 11] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_p_p_s[(n+1) * 9 + 5], PRIM_INT__s_p_d_s[n * 18 + 11]);
590                         PRIM_INT__s_p_d_s[n * 18 + 11] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_p_s_s[(n+1) * 3 + 1], PRIM_INT__s_p_s_s[n * 3 + 1]), PRIM_INT__s_p_d_s[n * 18 + 11]);
591 
592                         PRIM_INT__s_p_d_s[n * 18 + 12] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_p_p_s[n * 9 + 6]);
593                         PRIM_INT__s_p_d_s[n * 18 + 12] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_p_p_s[(n+1) * 9 + 6], PRIM_INT__s_p_d_s[n * 18 + 12]);
594                         PRIM_INT__s_p_d_s[n * 18 + 12] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_p_s_s[(n+1) * 3 + 2], PRIM_INT__s_p_s_s[n * 3 + 2]), PRIM_INT__s_p_d_s[n * 18 + 12]);
595 
596                         PRIM_INT__s_p_d_s[n * 18 + 15] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_p_p_s[n * 9 + 7]);
597                         PRIM_INT__s_p_d_s[n * 18 + 15] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_p_p_s[(n+1) * 9 + 7], PRIM_INT__s_p_d_s[n * 18 + 15]);
598                         PRIM_INT__s_p_d_s[n * 18 + 15] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_p_s_s[(n+1) * 3 + 2], PRIM_INT__s_p_s_s[n * 3 + 2]), PRIM_INT__s_p_d_s[n * 18 + 15]);
599 
600                         PRIM_INT__s_p_d_s[n * 18 + 17] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_p_p_s[n * 9 + 8]);
601                         PRIM_INT__s_p_d_s[n * 18 + 17] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_p_p_s[(n+1) * 9 + 8], PRIM_INT__s_p_d_s[n * 18 + 17]);
602                         PRIM_INT__s_p_d_s[n * 18 + 17] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_p_s_s[(n+1) * 3 + 2], PRIM_INT__s_p_s_s[n * 3 + 2]), PRIM_INT__s_p_d_s[n * 18 + 17]);
603                         PRIM_INT__s_p_d_s[n * 18 + 17] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__s_p_d_s[n * 18 + 17]);
604 
605                     }
606 
607 
608                     VRR_K_s_d_f_s(
609                             PRIM_INT__s_d_f_s,
610                             PRIM_INT__s_d_d_s,
611                             PRIM_INT__s_d_p_s,
612                             PRIM_INT__s_p_d_s,
613                             Q_PA,
614                             a_over_q,
615                             aoq_PQ,
616                             one_over_2pq,
617                             one_over_2q,
618                             2);
619 
620 
621                     ostei_general_vrr_K(0, 3, 4, 0, 1,
622                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
623                             PRIM_INT__s_f_f_s, PRIM_INT__s_f_d_s, NULL, NULL, PRIM_INT__s_d_f_s, PRIM_INT__s_f_g_s);
624 
625 
626                     ostei_general_vrr_K(0, 4, 3, 0, 2,
627                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
628                             PRIM_INT__s_g_d_s, PRIM_INT__s_g_p_s, NULL, NULL, PRIM_INT__s_f_d_s, PRIM_INT__s_g_f_s);
629 
630 
631                     VRR_J_s_h_s_s(
632                             PRIM_INT__s_h_s_s,
633                             PRIM_INT__s_g_s_s,
634                             PRIM_INT__s_f_s_s,
635                             P_PB,
636                             a_over_p,
637                             aop_PQ,
638                             one_over_2p,
639                             5);
640 
641 
642                     ostei_general_vrr_K(0, 5, 1, 0, 4,
643                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
644                             PRIM_INT__s_h_s_s, NULL, NULL, NULL, PRIM_INT__s_g_s_s, PRIM_INT__s_h_p_s);
645 
646 
647                     ostei_general_vrr_K(0, 5, 2, 0, 3,
648                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
649                             PRIM_INT__s_h_p_s, PRIM_INT__s_h_s_s, NULL, NULL, PRIM_INT__s_g_p_s, PRIM_INT__s_h_d_s);
650 
651 
652                     ostei_general_vrr_K(0, 4, 4, 0, 1,
653                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
654                             PRIM_INT__s_g_f_s, PRIM_INT__s_g_d_s, NULL, NULL, PRIM_INT__s_f_f_s, PRIM_INT__s_g_g_s);
655 
656 
657                     ostei_general_vrr_K(0, 5, 3, 0, 2,
658                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
659                             PRIM_INT__s_h_d_s, PRIM_INT__s_h_p_s, NULL, NULL, PRIM_INT__s_g_d_s, PRIM_INT__s_h_f_s);
660 
661 
662                     ostei_general_vrr_K(0, 5, 4, 0, 1,
663                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
664                             PRIM_INT__s_h_f_s, PRIM_INT__s_h_d_s, NULL, NULL, PRIM_INT__s_g_f_s, PRIM_INT__s_h_g_s);
665 
666 
667 
668 
669                     ////////////////////////////////////
670                     // Accumulate contracted integrals
671                     ////////////////////////////////////
672                     if(lastoffset == 0)
673                     {
674                         contract_all(60, PRIM_INT__s_f_d_s, PRIM_PTR_INT__s_f_d_s);
675                         contract_all(100, PRIM_INT__s_f_f_s, PRIM_PTR_INT__s_f_f_s);
676                         contract_all(150, PRIM_INT__s_f_g_s, PRIM_PTR_INT__s_f_g_s);
677                         contract_all(90, PRIM_INT__s_g_d_s, PRIM_PTR_INT__s_g_d_s);
678                         contract_all(150, PRIM_INT__s_g_f_s, PRIM_PTR_INT__s_g_f_s);
679                         contract_all(225, PRIM_INT__s_g_g_s, PRIM_PTR_INT__s_g_g_s);
680                         contract_all(126, PRIM_INT__s_h_d_s, PRIM_PTR_INT__s_h_d_s);
681                         contract_all(210, PRIM_INT__s_h_f_s, PRIM_PTR_INT__s_h_f_s);
682                         contract_all(315, PRIM_INT__s_h_g_s, PRIM_PTR_INT__s_h_g_s);
683                     }
684                     else
685                     {
686                         contract(60, shelloffsets, PRIM_INT__s_f_d_s, PRIM_PTR_INT__s_f_d_s);
687                         contract(100, shelloffsets, PRIM_INT__s_f_f_s, PRIM_PTR_INT__s_f_f_s);
688                         contract(150, shelloffsets, PRIM_INT__s_f_g_s, PRIM_PTR_INT__s_f_g_s);
689                         contract(90, shelloffsets, PRIM_INT__s_g_d_s, PRIM_PTR_INT__s_g_d_s);
690                         contract(150, shelloffsets, PRIM_INT__s_g_f_s, PRIM_PTR_INT__s_g_f_s);
691                         contract(225, shelloffsets, PRIM_INT__s_g_g_s, PRIM_PTR_INT__s_g_g_s);
692                         contract(126, shelloffsets, PRIM_INT__s_h_d_s, PRIM_PTR_INT__s_h_d_s);
693                         contract(210, shelloffsets, PRIM_INT__s_h_f_s, PRIM_PTR_INT__s_h_f_s);
694                         contract(315, shelloffsets, PRIM_INT__s_h_g_s, PRIM_PTR_INT__s_h_g_s);
695                         PRIM_PTR_INT__s_f_d_s += lastoffset*60;
696                         PRIM_PTR_INT__s_f_f_s += lastoffset*100;
697                         PRIM_PTR_INT__s_f_g_s += lastoffset*150;
698                         PRIM_PTR_INT__s_g_d_s += lastoffset*90;
699                         PRIM_PTR_INT__s_g_f_s += lastoffset*150;
700                         PRIM_PTR_INT__s_g_g_s += lastoffset*225;
701                         PRIM_PTR_INT__s_h_d_s += lastoffset*126;
702                         PRIM_PTR_INT__s_h_f_s += lastoffset*210;
703                         PRIM_PTR_INT__s_h_g_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__s_f_d_s = INT__s_f_d_s + abcd * 60;
726                 double const * restrict HRR_INT__s_f_f_s = INT__s_f_f_s + abcd * 100;
727                 double const * restrict HRR_INT__s_f_g_s = INT__s_f_g_s + abcd * 150;
728                 double const * restrict HRR_INT__s_g_d_s = INT__s_g_d_s + abcd * 90;
729                 double const * restrict HRR_INT__s_g_f_s = INT__s_g_f_s + abcd * 150;
730                 double const * restrict HRR_INT__s_g_g_s = INT__s_g_g_s + abcd * 225;
731                 double const * restrict HRR_INT__s_h_d_s = INT__s_h_d_s + abcd * 126;
732                 double const * restrict HRR_INT__s_h_f_s = INT__s_h_f_s + abcd * 210;
733                 double const * restrict HRR_INT__s_h_g_s = INT__s_h_g_s + abcd * 315;
734                 double * restrict HRR_INT__d_f_d_d = INT__d_f_d_d + real_abcd * 2160;
735 
736                 // form INT__p_f_d_s
737                 HRR_I_p_f(
738                     HRR_INT__p_f_d_s,
739                     HRR_INT__s_f_d_s,
740                     HRR_INT__s_g_d_s,
741                     hAB, 6);
742 
743                 // form INT__p_f_f_s
744                 HRR_I_p_f(
745                     HRR_INT__p_f_f_s,
746                     HRR_INT__s_f_f_s,
747                     HRR_INT__s_g_f_s,
748                     hAB, 10);
749 
750                 // form INT__p_f_g_s
751                 HRR_I_p_f(
752                     HRR_INT__p_f_g_s,
753                     HRR_INT__s_f_g_s,
754                     HRR_INT__s_g_g_s,
755                     hAB, 15);
756 
757                 // form INT__p_g_d_s
758                 HRR_I_p_g(
759                     HRR_INT__p_g_d_s,
760                     HRR_INT__s_g_d_s,
761                     HRR_INT__s_h_d_s,
762                     hAB, 6);
763 
764                 // form INT__p_g_f_s
765                 HRR_I_p_g(
766                     HRR_INT__p_g_f_s,
767                     HRR_INT__s_g_f_s,
768                     HRR_INT__s_h_f_s,
769                     hAB, 10);
770 
771                 // form INT__p_g_g_s
772                 HRR_I_p_g(
773                     HRR_INT__p_g_g_s,
774                     HRR_INT__s_g_g_s,
775                     HRR_INT__s_h_g_s,
776                     hAB, 15);
777 
778                 // form INT__d_f_d_s
779                 HRR_I_d_f(
780                     HRR_INT__d_f_d_s,
781                     HRR_INT__p_f_d_s,
782                     HRR_INT__p_g_d_s,
783                     hAB, 6);
784 
785                 // form INT__d_f_f_s
786                 HRR_I_d_f(
787                     HRR_INT__d_f_f_s,
788                     HRR_INT__p_f_f_s,
789                     HRR_INT__p_g_f_s,
790                     hAB, 10);
791 
792                 // form INT__d_f_g_s
793                 HRR_I_d_f(
794                     HRR_INT__d_f_g_s,
795                     HRR_INT__p_f_g_s,
796                     HRR_INT__p_g_g_s,
797                     hAB, 15);
798 
799                 // form INT__d_f_d_p
800                 for(ibra = 0; ibra < 60; ++ibra)
801                 {
802                     HRR_INT__d_f_d_p[ibra * 18 + 0] = HRR_INT__d_f_f_s[ibra * 10 + 0] + ( hCD[0] * HRR_INT__d_f_d_s[ibra * 6 + 0] );
803 
804                     HRR_INT__d_f_d_p[ibra * 18 + 1] = HRR_INT__d_f_f_s[ibra * 10 + 1] + ( hCD[1] * HRR_INT__d_f_d_s[ibra * 6 + 0] );
805 
806                     HRR_INT__d_f_d_p[ibra * 18 + 2] = HRR_INT__d_f_f_s[ibra * 10 + 2] + ( hCD[2] * HRR_INT__d_f_d_s[ibra * 6 + 0] );
807 
808                     HRR_INT__d_f_d_p[ibra * 18 + 3] = HRR_INT__d_f_f_s[ibra * 10 + 1] + ( hCD[0] * HRR_INT__d_f_d_s[ibra * 6 + 1] );
809 
810                     HRR_INT__d_f_d_p[ibra * 18 + 4] = HRR_INT__d_f_f_s[ibra * 10 + 3] + ( hCD[1] * HRR_INT__d_f_d_s[ibra * 6 + 1] );
811 
812                     HRR_INT__d_f_d_p[ibra * 18 + 5] = HRR_INT__d_f_f_s[ibra * 10 + 4] + ( hCD[2] * HRR_INT__d_f_d_s[ibra * 6 + 1] );
813 
814                     HRR_INT__d_f_d_p[ibra * 18 + 6] = HRR_INT__d_f_f_s[ibra * 10 + 2] + ( hCD[0] * HRR_INT__d_f_d_s[ibra * 6 + 2] );
815 
816                     HRR_INT__d_f_d_p[ibra * 18 + 7] = HRR_INT__d_f_f_s[ibra * 10 + 4] + ( hCD[1] * HRR_INT__d_f_d_s[ibra * 6 + 2] );
817 
818                     HRR_INT__d_f_d_p[ibra * 18 + 8] = HRR_INT__d_f_f_s[ibra * 10 + 5] + ( hCD[2] * HRR_INT__d_f_d_s[ibra * 6 + 2] );
819 
820                     HRR_INT__d_f_d_p[ibra * 18 + 9] = HRR_INT__d_f_f_s[ibra * 10 + 3] + ( hCD[0] * HRR_INT__d_f_d_s[ibra * 6 + 3] );
821 
822                     HRR_INT__d_f_d_p[ibra * 18 + 10] = HRR_INT__d_f_f_s[ibra * 10 + 6] + ( hCD[1] * HRR_INT__d_f_d_s[ibra * 6 + 3] );
823 
824                     HRR_INT__d_f_d_p[ibra * 18 + 11] = HRR_INT__d_f_f_s[ibra * 10 + 7] + ( hCD[2] * HRR_INT__d_f_d_s[ibra * 6 + 3] );
825 
826                     HRR_INT__d_f_d_p[ibra * 18 + 12] = HRR_INT__d_f_f_s[ibra * 10 + 4] + ( hCD[0] * HRR_INT__d_f_d_s[ibra * 6 + 4] );
827 
828                     HRR_INT__d_f_d_p[ibra * 18 + 13] = HRR_INT__d_f_f_s[ibra * 10 + 7] + ( hCD[1] * HRR_INT__d_f_d_s[ibra * 6 + 4] );
829 
830                     HRR_INT__d_f_d_p[ibra * 18 + 14] = HRR_INT__d_f_f_s[ibra * 10 + 8] + ( hCD[2] * HRR_INT__d_f_d_s[ibra * 6 + 4] );
831 
832                     HRR_INT__d_f_d_p[ibra * 18 + 15] = HRR_INT__d_f_f_s[ibra * 10 + 5] + ( hCD[0] * HRR_INT__d_f_d_s[ibra * 6 + 5] );
833 
834                     HRR_INT__d_f_d_p[ibra * 18 + 16] = HRR_INT__d_f_f_s[ibra * 10 + 8] + ( hCD[1] * HRR_INT__d_f_d_s[ibra * 6 + 5] );
835 
836                     HRR_INT__d_f_d_p[ibra * 18 + 17] = HRR_INT__d_f_f_s[ibra * 10 + 9] + ( hCD[2] * HRR_INT__d_f_d_s[ibra * 6 + 5] );
837 
838                 }
839 
840                 // form INT__d_f_f_p
841                 HRR_L_f_p(
842                     HRR_INT__d_f_f_p,
843                     HRR_INT__d_f_f_s,
844                     HRR_INT__d_f_g_s,
845                     hCD, 60);
846 
847                 // form INT__d_f_d_d
848                 HRR_L_d_d(
849                     HRR_INT__d_f_d_d,
850                     HRR_INT__d_f_d_p,
851                     HRR_INT__d_f_f_p,
852                     hCD, 60);
853 
854 
855             }  // close HRR loop
856 
857 
858         }   // close loop cdbatch
859 
860         istart = iend;
861     }  // close loop over ab
862 
863     return P.nshell12_clip * Q.nshell12_clip;
864 }
865 
866