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_g_d_g_d(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__g_d_g_d)8 int ostei_g_d_g_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__g_d_g_d)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__g_d_g_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__g_s_g_s = work + (SIMINT_NSHELL_SIMD * 0);
30     double * const INT__g_s_h_s = work + (SIMINT_NSHELL_SIMD * 225);
31     double * const INT__g_s_i_s = work + (SIMINT_NSHELL_SIMD * 540);
32     double * const INT__h_s_g_s = work + (SIMINT_NSHELL_SIMD * 960);
33     double * const INT__h_s_h_s = work + (SIMINT_NSHELL_SIMD * 1275);
34     double * const INT__h_s_i_s = work + (SIMINT_NSHELL_SIMD * 1716);
35     double * const INT__i_s_g_s = work + (SIMINT_NSHELL_SIMD * 2304);
36     double * const INT__i_s_h_s = work + (SIMINT_NSHELL_SIMD * 2724);
37     double * const INT__i_s_i_s = work + (SIMINT_NSHELL_SIMD * 3312);
38     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*4096);
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 + 13;
41     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_d_s = primwork + 31;
42     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_s = primwork + 61;
43     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_p_s = primwork + 97;
44     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_d_s = primwork + 151;
45     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_f_s = primwork + 241;
46     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 361;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_p_s = primwork + 427;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_d_s = primwork + 535;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_f_s = primwork + 715;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_g_s = primwork + 955;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_s = primwork + 1225;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_p_s = primwork + 1325;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_d_s = primwork + 1505;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_f_s = primwork + 1805;
55     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_g_s = primwork + 2205;
56     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_h_s = primwork + 2655;
57     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_s_s = primwork + 3075;
58     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_p_s = primwork + 3210;
59     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_d_s = primwork + 3480;
60     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_f_s = primwork + 3930;
61     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_g_s = primwork + 4530;
62     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_h_s = primwork + 5205;
63     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_i_s = primwork + 5835;
64     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_s_s = primwork + 6255;
65     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_p_s = primwork + 6423;
66     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_d_s = primwork + 6801;
67     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_f_s = primwork + 7431;
68     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_g_s = primwork + 8271;
69     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_h_s = primwork + 9216;
70     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_i_s = primwork + 10098;
71     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_s_s = primwork + 10686;
72     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_p_s = primwork + 10882;
73     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_d_s = primwork + 11386;
74     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_f_s = primwork + 12226;
75     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_g_s = primwork + 13346;
76     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_h_s = primwork + 14606;
77     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_i_s = primwork + 15782;
78     double * const hrrwork = (double *)(primwork + 16566);
79     double * const HRR_INT__g_p_g_s = hrrwork + 0;
80     double * const HRR_INT__g_p_h_s = hrrwork + 675;
81     double * const HRR_INT__g_p_i_s = hrrwork + 1620;
82     double * const HRR_INT__g_d_g_s = hrrwork + 2880;
83     double * const HRR_INT__g_d_g_p = hrrwork + 4230;
84     double * const HRR_INT__g_d_h_s = hrrwork + 8280;
85     double * const HRR_INT__g_d_h_p = hrrwork + 10170;
86     double * const HRR_INT__g_d_i_s = hrrwork + 15840;
87     double * const HRR_INT__h_p_g_s = hrrwork + 18360;
88     double * const HRR_INT__h_p_h_s = hrrwork + 19305;
89     double * const HRR_INT__h_p_i_s = hrrwork + 20628;
90 
91 
92     // Create constants
93     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
94     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
95     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
96     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
97     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
98     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
99     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
100 
101 
102     ////////////////////////////////////////
103     // Loop over shells and primitives
104     ////////////////////////////////////////
105 
106     real_abcd = 0;
107     istart = 0;
108     for(ab = 0; ab < P.nshell12_clip; ++ab)
109     {
110         const int iend = istart + P.nprim12[ab];
111 
112         cd = 0;
113         jstart = 0;
114 
115         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
116         {
117             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
118             int jend = jstart;
119             for(i = 0; i < nshellbatch; i++)
120                 jend += Q.nprim12[cd+i];
121 
122             // Clear the beginning of the workspace (where we are accumulating integrals)
123             memset(work, 0, SIMINT_NSHELL_SIMD * 4096 * sizeof(double));
124             abcd = 0;
125 
126 
127             for(i = istart; i < iend; ++i)
128             {
129                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
130 
131                 if(check_screen)
132                 {
133                     // Skip this whole thing if always insignificant
134                     if((P.screen[i] * Q.screen_max) < screen_tol)
135                         continue;
136                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
137                 }
138 
139                 icd = 0;
140                 iprimcd = 0;
141                 nprim_icd = Q.nprim12[cd];
142                 double * restrict PRIM_PTR_INT__g_s_g_s = INT__g_s_g_s + abcd * 225;
143                 double * restrict PRIM_PTR_INT__g_s_h_s = INT__g_s_h_s + abcd * 315;
144                 double * restrict PRIM_PTR_INT__g_s_i_s = INT__g_s_i_s + abcd * 420;
145                 double * restrict PRIM_PTR_INT__h_s_g_s = INT__h_s_g_s + abcd * 315;
146                 double * restrict PRIM_PTR_INT__h_s_h_s = INT__h_s_h_s + abcd * 441;
147                 double * restrict PRIM_PTR_INT__h_s_i_s = INT__h_s_i_s + abcd * 588;
148                 double * restrict PRIM_PTR_INT__i_s_g_s = INT__i_s_g_s + abcd * 420;
149                 double * restrict PRIM_PTR_INT__i_s_h_s = INT__i_s_h_s + abcd * 588;
150                 double * restrict PRIM_PTR_INT__i_s_i_s = INT__i_s_i_s + abcd * 784;
151 
152 
153 
154                 // Load these one per loop over i
155                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
156                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
157                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
158 
159                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
160 
161                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
162                 {
163                     // calculate the shell offsets
164                     // these are the offset from the shell pointed to by cd
165                     // for each element
166                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
167                     int lastoffset = 0;
168                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
169 
170                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
171                     {
172                         // Handle if the first element of the vector is a new shell
173                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
174                         {
175                             nprim_icd += Q.nprim12[cd + (++icd)];
176                             PRIM_PTR_INT__g_s_g_s += 225;
177                             PRIM_PTR_INT__g_s_h_s += 315;
178                             PRIM_PTR_INT__g_s_i_s += 420;
179                             PRIM_PTR_INT__h_s_g_s += 315;
180                             PRIM_PTR_INT__h_s_h_s += 441;
181                             PRIM_PTR_INT__h_s_i_s += 588;
182                             PRIM_PTR_INT__i_s_g_s += 420;
183                             PRIM_PTR_INT__i_s_h_s += 588;
184                             PRIM_PTR_INT__i_s_i_s += 784;
185                         }
186                         iprimcd++;
187                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
188                         {
189                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
190                             {
191                                 shelloffsets[n] = shelloffsets[n-1] + 1;
192                                 lastoffset++;
193                                 nprim_icd += Q.nprim12[cd + (++icd)];
194                             }
195                             else
196                                 shelloffsets[n] = shelloffsets[n-1];
197                             iprimcd++;
198                         }
199                     }
200                     else
201                         iprimcd += SIMINT_SIMD_LEN;
202 
203                     // Do we have to compute this vector (or has it been screened out)?
204                     // (not_screened != 0 means we have to do this vector)
205                     if(check_screen)
206                     {
207                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
208                         if(vmax < screen_tol)
209                         {
210                             PRIM_PTR_INT__g_s_g_s += lastoffset*225;
211                             PRIM_PTR_INT__g_s_h_s += lastoffset*315;
212                             PRIM_PTR_INT__g_s_i_s += lastoffset*420;
213                             PRIM_PTR_INT__h_s_g_s += lastoffset*315;
214                             PRIM_PTR_INT__h_s_h_s += lastoffset*441;
215                             PRIM_PTR_INT__h_s_i_s += lastoffset*588;
216                             PRIM_PTR_INT__i_s_g_s += lastoffset*420;
217                             PRIM_PTR_INT__i_s_h_s += lastoffset*588;
218                             PRIM_PTR_INT__i_s_i_s += lastoffset*784;
219                             continue;
220                         }
221                     }
222 
223                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
224                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
225                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
226                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
227 
228 
229                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
230                     SIMINT_DBLTYPE PQ[3];
231                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
232                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
233                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
234                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
235                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
236                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
237 
238                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
239                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
240                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
241                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
242                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
243                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
244                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
245 
246                     // NOTE: Minus sign!
247                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
248                     SIMINT_DBLTYPE aop_PQ[3];
249                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
250                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
251                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
252 
253                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
254                     SIMINT_DBLTYPE aoq_PQ[3];
255                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
256                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
257                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
258                     // Put a minus sign here so we don't have to in RR routines
259                     a_over_q = SIMINT_NEG(a_over_q);
260 
261 
262                     //////////////////////////////////////////////
263                     // Fjt function section
264                     // Maximum v value: 12
265                     //////////////////////////////////////////////
266                     // The parameter to the Fjt function
267                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
268 
269 
270                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
271 
272 
273                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 12);
274                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
275                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
276                     for(n = 0; n <= 12; n++)
277                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
278 
279                     //////////////////////////////////////////////
280                     // Primitive integrals: Vertical recurrance
281                     //////////////////////////////////////////////
282 
283                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
284                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
285                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
286                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
287                     const SIMINT_DBLTYPE vrr_const_5_over_2p = SIMINT_MUL(const_5, one_over_2p);
288                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
289                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
290                     const SIMINT_DBLTYPE vrr_const_3_over_2q = SIMINT_MUL(const_3, one_over_2q);
291                     const SIMINT_DBLTYPE vrr_const_4_over_2q = SIMINT_MUL(const_4, one_over_2q);
292                     const SIMINT_DBLTYPE vrr_const_5_over_2q = SIMINT_MUL(const_5, one_over_2q);
293                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
294                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
295                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
296                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
297                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
298                     const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
299 
300 
301 
302                     // Forming PRIM_INT__p_s_s_s[12 * 3];
303                     for(n = 0; n < 12; ++n)  // loop over orders of auxiliary function
304                     {
305 
306                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
307                         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]);
308 
309                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
310                         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]);
311 
312                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
313                         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]);
314 
315                     }
316 
317 
318 
319                     // Forming PRIM_INT__d_s_s_s[11 * 6];
320                     for(n = 0; n < 11; ++n)  // loop over orders of auxiliary function
321                     {
322 
323                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
324                         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]);
325                         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]);
326 
327                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
328                         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]);
329 
330                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
331                         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]);
332 
333                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
334                         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]);
335                         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]);
336 
337                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
338                         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]);
339 
340                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
341                         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]);
342                         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]);
343 
344                     }
345 
346 
347 
348                     // Forming PRIM_INT__f_s_s_s[10 * 10];
349                     for(n = 0; n < 10; ++n)  // loop over orders of auxiliary function
350                     {
351 
352                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
353                         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]);
354                         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]);
355 
356                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
357                         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]);
358 
359                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
360                         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]);
361 
362                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
363                         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]);
364 
365                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
366                         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]);
367 
368                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
369                         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]);
370 
371                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
372                         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]);
373                         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]);
374 
375                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
376                         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]);
377 
378                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
379                         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]);
380 
381                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
382                         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]);
383                         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]);
384 
385                     }
386 
387 
388                     VRR_I_g_s_s_s(
389                             PRIM_INT__g_s_s_s,
390                             PRIM_INT__f_s_s_s,
391                             PRIM_INT__d_s_s_s,
392                             P_PA,
393                             a_over_p,
394                             aop_PQ,
395                             one_over_2p,
396                             9);
397 
398 
399                     VRR_K_g_s_p_s(
400                             PRIM_INT__g_s_p_s,
401                             PRIM_INT__g_s_s_s,
402                             PRIM_INT__f_s_s_s,
403                             Q_PA,
404                             aoq_PQ,
405                             one_over_2pq,
406                             6);
407 
408 
409                     VRR_K_f_s_p_s(
410                             PRIM_INT__f_s_p_s,
411                             PRIM_INT__f_s_s_s,
412                             PRIM_INT__d_s_s_s,
413                             Q_PA,
414                             aoq_PQ,
415                             one_over_2pq,
416                             6);
417 
418 
419                     ostei_general_vrr_K(4, 0, 2, 0, 5,
420                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
421                             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);
422 
423 
424 
425                     // Forming PRIM_INT__d_s_p_s[6 * 18];
426                     for(n = 0; n < 6; ++n)  // loop over orders of auxiliary function
427                     {
428 
429                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
430                         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]);
431                         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]);
432 
433                         PRIM_INT__d_s_p_s[n * 18 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
434                         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]);
435 
436                         PRIM_INT__d_s_p_s[n * 18 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
437                         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]);
438 
439                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 1]);
440                         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]);
441                         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]);
442 
443                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 1]);
444                         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]);
445                         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]);
446 
447                         PRIM_INT__d_s_p_s[n * 18 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
448                         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]);
449 
450                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 2]);
451                         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]);
452                         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]);
453 
454                         PRIM_INT__d_s_p_s[n * 18 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 2]);
455                         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]);
456 
457                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 2]);
458                         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]);
459                         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]);
460 
461                         PRIM_INT__d_s_p_s[n * 18 + 9] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
462                         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]);
463 
464                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
465                         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]);
466                         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]);
467 
468                         PRIM_INT__d_s_p_s[n * 18 + 11] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
469                         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]);
470 
471                         PRIM_INT__d_s_p_s[n * 18 + 12] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 4]);
472                         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]);
473 
474                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 4]);
475                         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]);
476                         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]);
477 
478                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 4]);
479                         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]);
480                         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]);
481 
482                         PRIM_INT__d_s_p_s[n * 18 + 15] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
483                         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]);
484 
485                         PRIM_INT__d_s_p_s[n * 18 + 16] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
486                         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]);
487 
488                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
489                         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]);
490                         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]);
491 
492                     }
493 
494 
495                     VRR_K_f_s_d_s(
496                             PRIM_INT__f_s_d_s,
497                             PRIM_INT__f_s_p_s,
498                             PRIM_INT__f_s_s_s,
499                             PRIM_INT__d_s_p_s,
500                             Q_PA,
501                             a_over_q,
502                             aoq_PQ,
503                             one_over_2pq,
504                             one_over_2q,
505                             5);
506 
507 
508                     ostei_general_vrr_K(4, 0, 3, 0, 4,
509                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
510                             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);
511 
512 
513 
514                     // Forming PRIM_INT__p_s_p_s[6 * 9];
515                     for(n = 0; n < 6; ++n)  // loop over orders of auxiliary function
516                     {
517 
518                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
519                         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]);
520                         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]);
521 
522                         PRIM_INT__p_s_p_s[n * 9 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
523                         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]);
524 
525                         PRIM_INT__p_s_p_s[n * 9 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
526                         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]);
527 
528                         PRIM_INT__p_s_p_s[n * 9 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 1]);
529                         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]);
530 
531                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
532                         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]);
533                         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]);
534 
535                         PRIM_INT__p_s_p_s[n * 9 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
536                         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]);
537 
538                         PRIM_INT__p_s_p_s[n * 9 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 2]);
539                         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]);
540 
541                         PRIM_INT__p_s_p_s[n * 9 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 2]);
542                         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]);
543 
544                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
545                         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]);
546                         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]);
547 
548                     }
549 
550 
551                     VRR_K_d_s_d_s(
552                             PRIM_INT__d_s_d_s,
553                             PRIM_INT__d_s_p_s,
554                             PRIM_INT__d_s_s_s,
555                             PRIM_INT__p_s_p_s,
556                             Q_PA,
557                             a_over_q,
558                             aoq_PQ,
559                             one_over_2pq,
560                             one_over_2q,
561                             5);
562 
563 
564                     ostei_general_vrr_K(3, 0, 3, 0, 4,
565                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
566                             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);
567 
568 
569                     ostei_general_vrr_K(4, 0, 4, 0, 3,
570                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
571                             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);
572 
573 
574 
575                     // Forming PRIM_INT__s_s_p_s[6 * 3];
576                     for(n = 0; n < 6; ++n)  // loop over orders of auxiliary function
577                     {
578 
579                         PRIM_INT__s_s_p_s[n * 3 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
580                         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]);
581 
582                         PRIM_INT__s_s_p_s[n * 3 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
583                         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]);
584 
585                         PRIM_INT__s_s_p_s[n * 3 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
586                         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]);
587 
588                     }
589 
590 
591 
592                     // Forming PRIM_INT__p_s_d_s[5 * 18];
593                     for(n = 0; n < 5; ++n)  // loop over orders of auxiliary function
594                     {
595 
596                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_p_s[n * 9 + 0]);
597                         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]);
598                         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]);
599                         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]);
600 
601                         PRIM_INT__p_s_d_s[n * 18 + 3] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_p_s[n * 9 + 1]);
602                         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]);
603                         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]);
604 
605                         PRIM_INT__p_s_d_s[n * 18 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_p_s[n * 9 + 2]);
606                         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]);
607                         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]);
608 
609                         PRIM_INT__p_s_d_s[n * 18 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_p_s[n * 9 + 3]);
610                         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]);
611                         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]);
612 
613                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_p_s[n * 9 + 4]);
614                         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]);
615                         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]);
616                         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]);
617 
618                         PRIM_INT__p_s_d_s[n * 18 + 11] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_p_s[n * 9 + 5]);
619                         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]);
620                         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]);
621 
622                         PRIM_INT__p_s_d_s[n * 18 + 12] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_p_s[n * 9 + 6]);
623                         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]);
624                         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]);
625 
626                         PRIM_INT__p_s_d_s[n * 18 + 15] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_p_s[n * 9 + 7]);
627                         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]);
628                         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]);
629 
630                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_p_s[n * 9 + 8]);
631                         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]);
632                         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]);
633                         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]);
634 
635                     }
636 
637 
638                     VRR_K_d_s_f_s(
639                             PRIM_INT__d_s_f_s,
640                             PRIM_INT__d_s_d_s,
641                             PRIM_INT__d_s_p_s,
642                             PRIM_INT__p_s_d_s,
643                             Q_PA,
644                             a_over_q,
645                             aoq_PQ,
646                             one_over_2pq,
647                             one_over_2q,
648                             4);
649 
650 
651                     ostei_general_vrr_K(3, 0, 4, 0, 3,
652                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
653                             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);
654 
655 
656                     ostei_general_vrr_K(4, 0, 5, 0, 2,
657                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
658                             PRIM_INT__g_s_g_s, PRIM_INT__g_s_f_s, NULL, PRIM_INT__f_s_g_s, NULL, PRIM_INT__g_s_h_s);
659 
660 
661                     VRR_I_h_s_s_s(
662                             PRIM_INT__h_s_s_s,
663                             PRIM_INT__g_s_s_s,
664                             PRIM_INT__f_s_s_s,
665                             P_PA,
666                             a_over_p,
667                             aop_PQ,
668                             one_over_2p,
669                             8);
670 
671 
672                     ostei_general_vrr_K(5, 0, 1, 0, 6,
673                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
674                             PRIM_INT__h_s_s_s, NULL, NULL, PRIM_INT__g_s_s_s, NULL, PRIM_INT__h_s_p_s);
675 
676 
677                     ostei_general_vrr_K(5, 0, 2, 0, 5,
678                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
679                             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);
680 
681 
682                     ostei_general_vrr_K(5, 0, 3, 0, 4,
683                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
684                             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);
685 
686 
687                     ostei_general_vrr_K(5, 0, 4, 0, 3,
688                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
689                             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);
690 
691 
692 
693                     // Forming PRIM_INT__s_s_d_s[5 * 6];
694                     for(n = 0; n < 5; ++n)  // loop over orders of auxiliary function
695                     {
696 
697                         PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_p_s[n * 3 + 0]);
698                         PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_d_s[n * 6 + 0]);
699                         PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__s_s_d_s[n * 6 + 0]);
700 
701                         PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_p_s[n * 3 + 1]);
702                         PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_d_s[n * 6 + 3]);
703                         PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__s_s_d_s[n * 6 + 3]);
704 
705                         PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 2]);
706                         PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__s_s_d_s[n * 6 + 5]);
707                         PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__s_s_d_s[n * 6 + 5]);
708 
709                     }
710 
711 
712                     VRR_K_p_s_f_s(
713                             PRIM_INT__p_s_f_s,
714                             PRIM_INT__p_s_d_s,
715                             PRIM_INT__p_s_p_s,
716                             PRIM_INT__s_s_d_s,
717                             Q_PA,
718                             a_over_q,
719                             aoq_PQ,
720                             one_over_2pq,
721                             one_over_2q,
722                             4);
723 
724 
725                     ostei_general_vrr_K(2, 0, 4, 0, 3,
726                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
727                             PRIM_INT__d_s_f_s, PRIM_INT__d_s_d_s, NULL, PRIM_INT__p_s_f_s, NULL, PRIM_INT__d_s_g_s);
728 
729 
730                     ostei_general_vrr_K(3, 0, 5, 0, 2,
731                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
732                             PRIM_INT__f_s_g_s, PRIM_INT__f_s_f_s, NULL, PRIM_INT__d_s_g_s, NULL, PRIM_INT__f_s_h_s);
733 
734 
735                     ostei_general_vrr_K(4, 0, 6, 0, 1,
736                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
737                             PRIM_INT__g_s_h_s, PRIM_INT__g_s_g_s, NULL, PRIM_INT__f_s_h_s, NULL, PRIM_INT__g_s_i_s);
738 
739 
740                     ostei_general_vrr_K(5, 0, 5, 0, 2,
741                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
742                             PRIM_INT__h_s_g_s, PRIM_INT__h_s_f_s, NULL, PRIM_INT__g_s_g_s, NULL, PRIM_INT__h_s_h_s);
743 
744 
745                     ostei_general_vrr1_I(6, 7,
746                             one_over_2p, a_over_p, aop_PQ, P_PA,
747                             PRIM_INT__h_s_s_s, PRIM_INT__g_s_s_s, PRIM_INT__i_s_s_s);
748 
749 
750                     ostei_general_vrr_K(6, 0, 1, 0, 6,
751                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
752                             PRIM_INT__i_s_s_s, NULL, NULL, PRIM_INT__h_s_s_s, NULL, PRIM_INT__i_s_p_s);
753 
754 
755                     ostei_general_vrr_K(6, 0, 2, 0, 5,
756                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
757                             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);
758 
759 
760                     ostei_general_vrr_K(6, 0, 3, 0, 4,
761                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
762                             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);
763 
764 
765                     ostei_general_vrr_K(6, 0, 4, 0, 3,
766                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
767                             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);
768 
769 
770                     ostei_general_vrr_K(5, 0, 6, 0, 1,
771                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
772                             PRIM_INT__h_s_h_s, PRIM_INT__h_s_g_s, NULL, PRIM_INT__g_s_h_s, NULL, PRIM_INT__h_s_i_s);
773 
774 
775                     ostei_general_vrr_K(6, 0, 5, 0, 2,
776                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
777                             PRIM_INT__i_s_g_s, PRIM_INT__i_s_f_s, NULL, PRIM_INT__h_s_g_s, NULL, PRIM_INT__i_s_h_s);
778 
779 
780                     ostei_general_vrr_K(6, 0, 6, 0, 1,
781                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
782                             PRIM_INT__i_s_h_s, PRIM_INT__i_s_g_s, NULL, PRIM_INT__h_s_h_s, NULL, PRIM_INT__i_s_i_s);
783 
784 
785 
786 
787                     ////////////////////////////////////
788                     // Accumulate contracted integrals
789                     ////////////////////////////////////
790                     if(lastoffset == 0)
791                     {
792                         contract_all(225, PRIM_INT__g_s_g_s, PRIM_PTR_INT__g_s_g_s);
793                         contract_all(315, PRIM_INT__g_s_h_s, PRIM_PTR_INT__g_s_h_s);
794                         contract_all(420, PRIM_INT__g_s_i_s, PRIM_PTR_INT__g_s_i_s);
795                         contract_all(315, PRIM_INT__h_s_g_s, PRIM_PTR_INT__h_s_g_s);
796                         contract_all(441, PRIM_INT__h_s_h_s, PRIM_PTR_INT__h_s_h_s);
797                         contract_all(588, PRIM_INT__h_s_i_s, PRIM_PTR_INT__h_s_i_s);
798                         contract_all(420, PRIM_INT__i_s_g_s, PRIM_PTR_INT__i_s_g_s);
799                         contract_all(588, PRIM_INT__i_s_h_s, PRIM_PTR_INT__i_s_h_s);
800                         contract_all(784, PRIM_INT__i_s_i_s, PRIM_PTR_INT__i_s_i_s);
801                     }
802                     else
803                     {
804                         contract(225, shelloffsets, PRIM_INT__g_s_g_s, PRIM_PTR_INT__g_s_g_s);
805                         contract(315, shelloffsets, PRIM_INT__g_s_h_s, PRIM_PTR_INT__g_s_h_s);
806                         contract(420, shelloffsets, PRIM_INT__g_s_i_s, PRIM_PTR_INT__g_s_i_s);
807                         contract(315, shelloffsets, PRIM_INT__h_s_g_s, PRIM_PTR_INT__h_s_g_s);
808                         contract(441, shelloffsets, PRIM_INT__h_s_h_s, PRIM_PTR_INT__h_s_h_s);
809                         contract(588, shelloffsets, PRIM_INT__h_s_i_s, PRIM_PTR_INT__h_s_i_s);
810                         contract(420, shelloffsets, PRIM_INT__i_s_g_s, PRIM_PTR_INT__i_s_g_s);
811                         contract(588, shelloffsets, PRIM_INT__i_s_h_s, PRIM_PTR_INT__i_s_h_s);
812                         contract(784, shelloffsets, PRIM_INT__i_s_i_s, PRIM_PTR_INT__i_s_i_s);
813                         PRIM_PTR_INT__g_s_g_s += lastoffset*225;
814                         PRIM_PTR_INT__g_s_h_s += lastoffset*315;
815                         PRIM_PTR_INT__g_s_i_s += lastoffset*420;
816                         PRIM_PTR_INT__h_s_g_s += lastoffset*315;
817                         PRIM_PTR_INT__h_s_h_s += lastoffset*441;
818                         PRIM_PTR_INT__h_s_i_s += lastoffset*588;
819                         PRIM_PTR_INT__i_s_g_s += lastoffset*420;
820                         PRIM_PTR_INT__i_s_h_s += lastoffset*588;
821                         PRIM_PTR_INT__i_s_i_s += lastoffset*784;
822                     }
823 
824                 }  // close loop over j
825             }  // close loop over i
826 
827             //Advance to the next batch
828             jstart = SIMINT_SIMD_ROUND(jend);
829 
830             //////////////////////////////////////////////
831             // Contracted integrals: Horizontal recurrance
832             //////////////////////////////////////////////
833 
834 
835             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
836 
837 
838             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
839             {
840                 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
841 
842                 // set up HRR pointers
843                 double const * restrict HRR_INT__g_s_g_s = INT__g_s_g_s + abcd * 225;
844                 double const * restrict HRR_INT__g_s_h_s = INT__g_s_h_s + abcd * 315;
845                 double const * restrict HRR_INT__g_s_i_s = INT__g_s_i_s + abcd * 420;
846                 double const * restrict HRR_INT__h_s_g_s = INT__h_s_g_s + abcd * 315;
847                 double const * restrict HRR_INT__h_s_h_s = INT__h_s_h_s + abcd * 441;
848                 double const * restrict HRR_INT__h_s_i_s = INT__h_s_i_s + abcd * 588;
849                 double const * restrict HRR_INT__i_s_g_s = INT__i_s_g_s + abcd * 420;
850                 double const * restrict HRR_INT__i_s_h_s = INT__i_s_h_s + abcd * 588;
851                 double const * restrict HRR_INT__i_s_i_s = INT__i_s_i_s + abcd * 784;
852                 double * restrict HRR_INT__g_d_g_d = INT__g_d_g_d + real_abcd * 8100;
853 
854                 // form INT__g_p_g_s
855                 HRR_J_g_p(
856                     HRR_INT__g_p_g_s,
857                     HRR_INT__g_s_g_s,
858                     HRR_INT__h_s_g_s,
859                     hAB, 15);
860 
861                 // form INT__g_p_h_s
862                 HRR_J_g_p(
863                     HRR_INT__g_p_h_s,
864                     HRR_INT__g_s_h_s,
865                     HRR_INT__h_s_h_s,
866                     hAB, 21);
867 
868                 // form INT__g_p_i_s
869                 HRR_J_g_p(
870                     HRR_INT__g_p_i_s,
871                     HRR_INT__g_s_i_s,
872                     HRR_INT__h_s_i_s,
873                     hAB, 28);
874 
875                 // form INT__h_p_g_s
876                 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);
877 
878                 // form INT__h_p_h_s
879                 ostei_general_hrr_J(5, 1, 5, 0, hAB, HRR_INT__i_s_h_s, HRR_INT__h_s_h_s, HRR_INT__h_p_h_s);
880 
881                 // form INT__h_p_i_s
882                 ostei_general_hrr_J(5, 1, 6, 0, hAB, HRR_INT__i_s_i_s, HRR_INT__h_s_i_s, HRR_INT__h_p_i_s);
883 
884                 // form INT__g_d_g_s
885                 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);
886 
887                 // form INT__g_d_h_s
888                 ostei_general_hrr_J(4, 2, 5, 0, hAB, HRR_INT__h_p_h_s, HRR_INT__g_p_h_s, HRR_INT__g_d_h_s);
889 
890                 // form INT__g_d_i_s
891                 ostei_general_hrr_J(4, 2, 6, 0, hAB, HRR_INT__h_p_i_s, HRR_INT__g_p_i_s, HRR_INT__g_d_i_s);
892 
893                 // form INT__g_d_g_p
894                 HRR_L_g_p(
895                     HRR_INT__g_d_g_p,
896                     HRR_INT__g_d_g_s,
897                     HRR_INT__g_d_h_s,
898                     hCD, 90);
899 
900                 // form INT__g_d_h_p
901                 ostei_general_hrr_L(4, 2, 5, 1, hCD, HRR_INT__g_d_i_s, HRR_INT__g_d_h_s, HRR_INT__g_d_h_p);
902 
903                 // form INT__g_d_g_d
904                 ostei_general_hrr_L(4, 2, 4, 2, hCD, HRR_INT__g_d_h_p, HRR_INT__g_d_g_p, HRR_INT__g_d_g_d);
905 
906 
907             }  // close HRR loop
908 
909 
910         }   // close loop cdbatch
911 
912         istart = iend;
913     }  // close loop over ab
914 
915     return P.nshell12_clip * Q.nshell12_clip;
916 }
917 
ostei_d_g_g_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_g_g_d)918 int ostei_d_g_g_d(struct simint_multi_shellpair const P,
919                   struct simint_multi_shellpair const Q,
920                   double screen_tol,
921                   double * const restrict work,
922                   double * const restrict INT__d_g_g_d)
923 {
924     double P_AB[3*P.nshell12];
925     struct simint_multi_shellpair P_tmp = P;
926     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
927     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
928     P_tmp.AB_x = P_AB;
929     P_tmp.AB_y = P_AB + P.nshell12;
930     P_tmp.AB_z = P_AB + 2*P.nshell12;
931 
932     for(int i = 0; i < P.nshell12; i++)
933     {
934         P_tmp.AB_x[i] = -P.AB_x[i];
935         P_tmp.AB_y[i] = -P.AB_y[i];
936         P_tmp.AB_z[i] = -P.AB_z[i];
937     }
938 
939     int ret = ostei_g_d_g_d(P_tmp, Q, screen_tol, work, INT__d_g_g_d);
940     double buffer[8100] SIMINT_ALIGN_ARRAY_DBL;
941 
942     for(int q = 0; q < ret; q++)
943     {
944         int idx = 0;
945         for(int a = 0; a < 6; ++a)
946         for(int b = 0; b < 15; ++b)
947         for(int c = 0; c < 15; ++c)
948         for(int d = 0; d < 6; ++d)
949             buffer[idx++] = INT__d_g_g_d[q*8100+b*540+a*90+c*6+d];
950 
951         memcpy(INT__d_g_g_d+q*8100, buffer, 8100*sizeof(double));
952     }
953 
954     return ret;
955 }
956 
ostei_g_d_d_g(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__g_d_d_g)957 int ostei_g_d_d_g(struct simint_multi_shellpair const P,
958                   struct simint_multi_shellpair const Q,
959                   double screen_tol,
960                   double * const restrict work,
961                   double * const restrict INT__g_d_d_g)
962 {
963     double Q_AB[3*Q.nshell12];
964     struct simint_multi_shellpair Q_tmp = Q;
965     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
966     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
967     Q_tmp.AB_x = Q_AB;
968     Q_tmp.AB_y = Q_AB + Q.nshell12;
969     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
970 
971     for(int i = 0; i < Q.nshell12; i++)
972     {
973         Q_tmp.AB_x[i] = -Q.AB_x[i];
974         Q_tmp.AB_y[i] = -Q.AB_y[i];
975         Q_tmp.AB_z[i] = -Q.AB_z[i];
976     }
977 
978     int ret = ostei_g_d_g_d(P, Q_tmp, screen_tol, work, INT__g_d_d_g);
979     double buffer[8100] SIMINT_ALIGN_ARRAY_DBL;
980 
981     for(int q = 0; q < ret; q++)
982     {
983         int idx = 0;
984         for(int a = 0; a < 15; ++a)
985         for(int b = 0; b < 6; ++b)
986         for(int c = 0; c < 6; ++c)
987         for(int d = 0; d < 15; ++d)
988             buffer[idx++] = INT__g_d_d_g[q*8100+a*540+b*90+d*6+c];
989 
990         memcpy(INT__g_d_d_g+q*8100, buffer, 8100*sizeof(double));
991     }
992 
993     return ret;
994 }
995 
ostei_d_g_d_g(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__d_g_d_g)996 int ostei_d_g_d_g(struct simint_multi_shellpair const P,
997                   struct simint_multi_shellpair const Q,
998                   double screen_tol,
999                   double * const restrict work,
1000                   double * const restrict INT__d_g_d_g)
1001 {
1002     double P_AB[3*P.nshell12];
1003     struct simint_multi_shellpair P_tmp = P;
1004     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
1005     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
1006     P_tmp.AB_x = P_AB;
1007     P_tmp.AB_y = P_AB + P.nshell12;
1008     P_tmp.AB_z = P_AB + 2*P.nshell12;
1009 
1010     for(int i = 0; i < P.nshell12; i++)
1011     {
1012         P_tmp.AB_x[i] = -P.AB_x[i];
1013         P_tmp.AB_y[i] = -P.AB_y[i];
1014         P_tmp.AB_z[i] = -P.AB_z[i];
1015     }
1016 
1017     double Q_AB[3*Q.nshell12];
1018     struct simint_multi_shellpair Q_tmp = Q;
1019     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
1020     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
1021     Q_tmp.AB_x = Q_AB;
1022     Q_tmp.AB_y = Q_AB + Q.nshell12;
1023     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
1024 
1025     for(int i = 0; i < Q.nshell12; i++)
1026     {
1027         Q_tmp.AB_x[i] = -Q.AB_x[i];
1028         Q_tmp.AB_y[i] = -Q.AB_y[i];
1029         Q_tmp.AB_z[i] = -Q.AB_z[i];
1030     }
1031 
1032     int ret = ostei_g_d_g_d(P_tmp, Q_tmp, screen_tol, work, INT__d_g_d_g);
1033     double buffer[8100] SIMINT_ALIGN_ARRAY_DBL;
1034 
1035     for(int q = 0; q < ret; q++)
1036     {
1037         int idx = 0;
1038         for(int a = 0; a < 6; ++a)
1039         for(int b = 0; b < 15; ++b)
1040         for(int c = 0; c < 6; ++c)
1041         for(int d = 0; d < 15; ++d)
1042             buffer[idx++] = INT__d_g_d_g[q*8100+b*540+a*90+d*6+c];
1043 
1044         memcpy(INT__d_g_d_g+q*8100, buffer, 8100*sizeof(double));
1045     }
1046 
1047     return ret;
1048 }
1049 
1050