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_f_k_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_f_k_d)8 int ostei_g_f_k_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_f_k_d)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__g_f_k_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_k_s = work + (SIMINT_NSHELL_SIMD * 0);
30     double * const INT__g_s_l_s = work + (SIMINT_NSHELL_SIMD * 540);
31     double * const INT__g_s_m_s = work + (SIMINT_NSHELL_SIMD * 1215);
32     double * const INT__h_s_k_s = work + (SIMINT_NSHELL_SIMD * 2040);
33     double * const INT__h_s_l_s = work + (SIMINT_NSHELL_SIMD * 2796);
34     double * const INT__h_s_m_s = work + (SIMINT_NSHELL_SIMD * 3741);
35     double * const INT__i_s_k_s = work + (SIMINT_NSHELL_SIMD * 4896);
36     double * const INT__i_s_l_s = work + (SIMINT_NSHELL_SIMD * 5904);
37     double * const INT__i_s_m_s = work + (SIMINT_NSHELL_SIMD * 7164);
38     double * const INT__k_s_k_s = work + (SIMINT_NSHELL_SIMD * 8704);
39     double * const INT__k_s_l_s = work + (SIMINT_NSHELL_SIMD * 10000);
40     double * const INT__k_s_m_s = work + (SIMINT_NSHELL_SIMD * 11620);
41     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*13600);
42     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
43     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_p_s = primwork + 17;
44     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_d_s = primwork + 65;
45     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_f_s = primwork + 155;
46     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_g_s = primwork + 295;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_h_s = primwork + 490;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_i_s = primwork + 742;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_k_s = primwork + 1050;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_l_s = primwork + 1410;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_m_s = primwork + 1815;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_p_s = primwork + 2255;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_d_s = primwork + 2318;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_f_s = primwork + 2444;
55     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_g_s = primwork + 2654;
56     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_h_s = primwork + 2969;
57     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_i_s = primwork + 3410;
58     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_k_s = primwork + 3998;
59     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_l_s = primwork + 4754;
60     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_m_s = primwork + 5699;
61     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_d_s = primwork + 6854;
62     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_f_s = primwork + 7070;
63     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_g_s = primwork + 7430;
64     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_h_s = primwork + 7970;
65     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_i_s = primwork + 8726;
66     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_k_s = primwork + 9734;
67     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_l_s = primwork + 11030;
68     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_m_s = primwork + 12650;
69     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_f_s = primwork + 14630;
70     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_g_s = primwork + 15130;
71     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_h_s = primwork + 15880;
72     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_i_s = primwork + 16930;
73     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_k_s = primwork + 18330;
74     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_l_s = primwork + 20130;
75     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_m_s = primwork + 22380;
76     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_g_s = primwork + 25130;
77     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_h_s = primwork + 26030;
78     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_i_s = primwork + 27290;
79     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_k_s = primwork + 28970;
80     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_l_s = primwork + 31130;
81     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_m_s = primwork + 33830;
82     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_h_s = primwork + 37130;
83     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_i_s = primwork + 38453;
84     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_k_s = primwork + 40217;
85     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_l_s = primwork + 42485;
86     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_m_s = primwork + 45320;
87     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_i_s = primwork + 48785;
88     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_k_s = primwork + 50353;
89     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_l_s = primwork + 52369;
90     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_m_s = primwork + 54889;
91     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_k_s = primwork + 57969;
92     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_l_s = primwork + 59265;
93     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_m_s = primwork + 60885;
94     double * const hrrwork = (double *)(primwork + 62865);
95     double * const HRR_INT__g_p_k_s = hrrwork + 0;
96     double * const HRR_INT__g_p_l_s = hrrwork + 1620;
97     double * const HRR_INT__g_p_m_s = hrrwork + 3645;
98     double * const HRR_INT__g_d_k_s = hrrwork + 6120;
99     double * const HRR_INT__g_d_l_s = hrrwork + 9360;
100     double * const HRR_INT__g_d_m_s = hrrwork + 13410;
101     double * const HRR_INT__g_f_k_s = hrrwork + 18360;
102     double * const HRR_INT__g_f_k_p = hrrwork + 23760;
103     double * const HRR_INT__g_f_l_s = hrrwork + 39960;
104     double * const HRR_INT__g_f_l_p = hrrwork + 46710;
105     double * const HRR_INT__g_f_m_s = hrrwork + 66960;
106     double * const HRR_INT__h_p_k_s = hrrwork + 75210;
107     double * const HRR_INT__h_p_l_s = hrrwork + 77478;
108     double * const HRR_INT__h_p_m_s = hrrwork + 80313;
109     double * const HRR_INT__h_d_k_s = hrrwork + 83778;
110     double * const HRR_INT__h_d_l_s = hrrwork + 88314;
111     double * const HRR_INT__h_d_m_s = hrrwork + 93984;
112     double * const HRR_INT__i_p_k_s = hrrwork + 100914;
113     double * const HRR_INT__i_p_l_s = hrrwork + 103938;
114     double * const HRR_INT__i_p_m_s = hrrwork + 107718;
115 
116 
117     // Create constants
118     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
119     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
120     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
121     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
122     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
123     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
124     const SIMINT_DBLTYPE const_7 = SIMINT_DBLSET1(7);
125     const SIMINT_DBLTYPE const_8 = SIMINT_DBLSET1(8);
126     const SIMINT_DBLTYPE const_9 = SIMINT_DBLSET1(9);
127     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
128 
129 
130     ////////////////////////////////////////
131     // Loop over shells and primitives
132     ////////////////////////////////////////
133 
134     real_abcd = 0;
135     istart = 0;
136     for(ab = 0; ab < P.nshell12_clip; ++ab)
137     {
138         const int iend = istart + P.nprim12[ab];
139 
140         cd = 0;
141         jstart = 0;
142 
143         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
144         {
145             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
146             int jend = jstart;
147             for(i = 0; i < nshellbatch; i++)
148                 jend += Q.nprim12[cd+i];
149 
150             // Clear the beginning of the workspace (where we are accumulating integrals)
151             memset(work, 0, SIMINT_NSHELL_SIMD * 13600 * sizeof(double));
152             abcd = 0;
153 
154 
155             for(i = istart; i < iend; ++i)
156             {
157                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
158 
159                 if(check_screen)
160                 {
161                     // Skip this whole thing if always insignificant
162                     if((P.screen[i] * Q.screen_max) < screen_tol)
163                         continue;
164                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
165                 }
166 
167                 icd = 0;
168                 iprimcd = 0;
169                 nprim_icd = Q.nprim12[cd];
170                 double * restrict PRIM_PTR_INT__g_s_k_s = INT__g_s_k_s + abcd * 540;
171                 double * restrict PRIM_PTR_INT__g_s_l_s = INT__g_s_l_s + abcd * 675;
172                 double * restrict PRIM_PTR_INT__g_s_m_s = INT__g_s_m_s + abcd * 825;
173                 double * restrict PRIM_PTR_INT__h_s_k_s = INT__h_s_k_s + abcd * 756;
174                 double * restrict PRIM_PTR_INT__h_s_l_s = INT__h_s_l_s + abcd * 945;
175                 double * restrict PRIM_PTR_INT__h_s_m_s = INT__h_s_m_s + abcd * 1155;
176                 double * restrict PRIM_PTR_INT__i_s_k_s = INT__i_s_k_s + abcd * 1008;
177                 double * restrict PRIM_PTR_INT__i_s_l_s = INT__i_s_l_s + abcd * 1260;
178                 double * restrict PRIM_PTR_INT__i_s_m_s = INT__i_s_m_s + abcd * 1540;
179                 double * restrict PRIM_PTR_INT__k_s_k_s = INT__k_s_k_s + abcd * 1296;
180                 double * restrict PRIM_PTR_INT__k_s_l_s = INT__k_s_l_s + abcd * 1620;
181                 double * restrict PRIM_PTR_INT__k_s_m_s = INT__k_s_m_s + abcd * 1980;
182 
183 
184 
185                 // Load these one per loop over i
186                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
187                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
188                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
189 
190                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
191 
192                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
193                 {
194                     // calculate the shell offsets
195                     // these are the offset from the shell pointed to by cd
196                     // for each element
197                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
198                     int lastoffset = 0;
199                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
200 
201                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
202                     {
203                         // Handle if the first element of the vector is a new shell
204                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
205                         {
206                             nprim_icd += Q.nprim12[cd + (++icd)];
207                             PRIM_PTR_INT__g_s_k_s += 540;
208                             PRIM_PTR_INT__g_s_l_s += 675;
209                             PRIM_PTR_INT__g_s_m_s += 825;
210                             PRIM_PTR_INT__h_s_k_s += 756;
211                             PRIM_PTR_INT__h_s_l_s += 945;
212                             PRIM_PTR_INT__h_s_m_s += 1155;
213                             PRIM_PTR_INT__i_s_k_s += 1008;
214                             PRIM_PTR_INT__i_s_l_s += 1260;
215                             PRIM_PTR_INT__i_s_m_s += 1540;
216                             PRIM_PTR_INT__k_s_k_s += 1296;
217                             PRIM_PTR_INT__k_s_l_s += 1620;
218                             PRIM_PTR_INT__k_s_m_s += 1980;
219                         }
220                         iprimcd++;
221                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
222                         {
223                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
224                             {
225                                 shelloffsets[n] = shelloffsets[n-1] + 1;
226                                 lastoffset++;
227                                 nprim_icd += Q.nprim12[cd + (++icd)];
228                             }
229                             else
230                                 shelloffsets[n] = shelloffsets[n-1];
231                             iprimcd++;
232                         }
233                     }
234                     else
235                         iprimcd += SIMINT_SIMD_LEN;
236 
237                     // Do we have to compute this vector (or has it been screened out)?
238                     // (not_screened != 0 means we have to do this vector)
239                     if(check_screen)
240                     {
241                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
242                         if(vmax < screen_tol)
243                         {
244                             PRIM_PTR_INT__g_s_k_s += lastoffset*540;
245                             PRIM_PTR_INT__g_s_l_s += lastoffset*675;
246                             PRIM_PTR_INT__g_s_m_s += lastoffset*825;
247                             PRIM_PTR_INT__h_s_k_s += lastoffset*756;
248                             PRIM_PTR_INT__h_s_l_s += lastoffset*945;
249                             PRIM_PTR_INT__h_s_m_s += lastoffset*1155;
250                             PRIM_PTR_INT__i_s_k_s += lastoffset*1008;
251                             PRIM_PTR_INT__i_s_l_s += lastoffset*1260;
252                             PRIM_PTR_INT__i_s_m_s += lastoffset*1540;
253                             PRIM_PTR_INT__k_s_k_s += lastoffset*1296;
254                             PRIM_PTR_INT__k_s_l_s += lastoffset*1620;
255                             PRIM_PTR_INT__k_s_m_s += lastoffset*1980;
256                             continue;
257                         }
258                     }
259 
260                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
261                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
262                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
263                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
264 
265 
266                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
267                     SIMINT_DBLTYPE PQ[3];
268                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
269                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
270                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
271                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
272                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
273                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
274 
275                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
276                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
277                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
278                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
279                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
280                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
281                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
282 
283                     // NOTE: Minus sign!
284                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
285                     SIMINT_DBLTYPE aop_PQ[3];
286                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
287                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
288                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
289 
290                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
291                     SIMINT_DBLTYPE aoq_PQ[3];
292                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
293                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
294                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
295                     // Put a minus sign here so we don't have to in RR routines
296                     a_over_q = SIMINT_NEG(a_over_q);
297 
298 
299                     //////////////////////////////////////////////
300                     // Fjt function section
301                     // Maximum v value: 16
302                     //////////////////////////////////////////////
303                     // The parameter to the Fjt function
304                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
305 
306 
307                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
308 
309 
310                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 16);
311                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
312                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
313                     for(n = 0; n <= 16; n++)
314                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
315 
316                     //////////////////////////////////////////////
317                     // Primitive integrals: Vertical recurrance
318                     //////////////////////////////////////////////
319 
320                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
321                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
322                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
323                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
324                     const SIMINT_DBLTYPE vrr_const_5_over_2p = SIMINT_MUL(const_5, one_over_2p);
325                     const SIMINT_DBLTYPE vrr_const_6_over_2p = SIMINT_MUL(const_6, one_over_2p);
326                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
327                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
328                     const SIMINT_DBLTYPE vrr_const_3_over_2q = SIMINT_MUL(const_3, one_over_2q);
329                     const SIMINT_DBLTYPE vrr_const_4_over_2q = SIMINT_MUL(const_4, one_over_2q);
330                     const SIMINT_DBLTYPE vrr_const_5_over_2q = SIMINT_MUL(const_5, one_over_2q);
331                     const SIMINT_DBLTYPE vrr_const_6_over_2q = SIMINT_MUL(const_6, one_over_2q);
332                     const SIMINT_DBLTYPE vrr_const_7_over_2q = SIMINT_MUL(const_7, one_over_2q);
333                     const SIMINT_DBLTYPE vrr_const_8_over_2q = SIMINT_MUL(const_8, one_over_2q);
334                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
335                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
336                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
337                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
338                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
339                     const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
340                     const SIMINT_DBLTYPE vrr_const_7_over_2pq = SIMINT_MUL(const_7, one_over_2pq);
341                     const SIMINT_DBLTYPE vrr_const_8_over_2pq = SIMINT_MUL(const_8, one_over_2pq);
342                     const SIMINT_DBLTYPE vrr_const_9_over_2pq = SIMINT_MUL(const_9, one_over_2pq);
343 
344 
345 
346                     // Forming PRIM_INT__s_s_p_s[16 * 3];
347                     for(n = 0; n < 16; ++n)  // loop over orders of auxiliary function
348                     {
349 
350                         PRIM_INT__s_s_p_s[n * 3 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
351                         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]);
352 
353                         PRIM_INT__s_s_p_s[n * 3 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
354                         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]);
355 
356                         PRIM_INT__s_s_p_s[n * 3 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
357                         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]);
358 
359                     }
360 
361 
362 
363                     // Forming PRIM_INT__s_s_d_s[15 * 6];
364                     for(n = 0; n < 15; ++n)  // loop over orders of auxiliary function
365                     {
366 
367                         PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_p_s[n * 3 + 0]);
368                         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]);
369                         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]);
370 
371                         PRIM_INT__s_s_d_s[n * 6 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_p_s[n * 3 + 0]);
372                         PRIM_INT__s_s_d_s[n * 6 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_d_s[n * 6 + 1]);
373 
374                         PRIM_INT__s_s_d_s[n * 6 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 0]);
375                         PRIM_INT__s_s_d_s[n * 6 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_d_s[n * 6 + 2]);
376 
377                         PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_p_s[n * 3 + 1]);
378                         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]);
379                         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]);
380 
381                         PRIM_INT__s_s_d_s[n * 6 + 4] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 1]);
382                         PRIM_INT__s_s_d_s[n * 6 + 4] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_d_s[n * 6 + 4]);
383 
384                         PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 2]);
385                         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]);
386                         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]);
387 
388                     }
389 
390 
391 
392                     // Forming PRIM_INT__s_s_f_s[14 * 10];
393                     for(n = 0; n < 14; ++n)  // loop over orders of auxiliary function
394                     {
395 
396                         PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 0]);
397                         PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__s_s_f_s[n * 10 + 0]);
398                         PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_p_s[n * 3 + 0]), PRIM_INT__s_s_f_s[n * 10 + 0]);
399 
400                         PRIM_INT__s_s_f_s[n * 10 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 0]);
401                         PRIM_INT__s_s_f_s[n * 10 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__s_s_f_s[n * 10 + 1]);
402 
403                         PRIM_INT__s_s_f_s[n * 10 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 0]);
404                         PRIM_INT__s_s_f_s[n * 10 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__s_s_f_s[n * 10 + 2]);
405 
406                         PRIM_INT__s_s_f_s[n * 10 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 3]);
407                         PRIM_INT__s_s_f_s[n * 10 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__s_s_f_s[n * 10 + 3]);
408 
409                         PRIM_INT__s_s_f_s[n * 10 + 4] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 1]);
410                         PRIM_INT__s_s_f_s[n * 10 + 4] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__s_s_f_s[n * 10 + 4]);
411 
412                         PRIM_INT__s_s_f_s[n * 10 + 5] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 5]);
413                         PRIM_INT__s_s_f_s[n * 10 + 5] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__s_s_f_s[n * 10 + 5]);
414 
415                         PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 3]);
416                         PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__s_s_f_s[n * 10 + 6]);
417                         PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_p_s[n * 3 + 1]), PRIM_INT__s_s_f_s[n * 10 + 6]);
418 
419                         PRIM_INT__s_s_f_s[n * 10 + 7] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 3]);
420                         PRIM_INT__s_s_f_s[n * 10 + 7] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__s_s_f_s[n * 10 + 7]);
421 
422                         PRIM_INT__s_s_f_s[n * 10 + 8] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 5]);
423                         PRIM_INT__s_s_f_s[n * 10 + 8] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__s_s_f_s[n * 10 + 8]);
424 
425                         PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 5]);
426                         PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__s_s_f_s[n * 10 + 9]);
427                         PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__s_s_p_s[n * 3 + 2]), PRIM_INT__s_s_f_s[n * 10 + 9]);
428 
429                     }
430 
431 
432                     VRR_K_s_s_g_s(
433                             PRIM_INT__s_s_g_s,
434                             PRIM_INT__s_s_f_s,
435                             PRIM_INT__s_s_d_s,
436                             Q_PA,
437                             a_over_q,
438                             aoq_PQ,
439                             one_over_2q,
440                             13);
441 
442 
443                     VRR_K_s_s_h_s(
444                             PRIM_INT__s_s_h_s,
445                             PRIM_INT__s_s_g_s,
446                             PRIM_INT__s_s_f_s,
447                             Q_PA,
448                             a_over_q,
449                             aoq_PQ,
450                             one_over_2q,
451                             12);
452 
453 
454                     ostei_general_vrr1_K(6, 11,
455                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
456                             PRIM_INT__s_s_h_s, PRIM_INT__s_s_g_s, PRIM_INT__s_s_i_s);
457 
458 
459                     ostei_general_vrr1_K(7, 10,
460                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
461                             PRIM_INT__s_s_i_s, PRIM_INT__s_s_h_s, PRIM_INT__s_s_k_s);
462 
463 
464                     ostei_general_vrr_I(1, 0, 7, 0, 7,
465                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
466                             PRIM_INT__s_s_k_s, NULL, NULL, PRIM_INT__s_s_i_s, NULL, PRIM_INT__p_s_k_s);
467 
468 
469                     ostei_general_vrr_I(1, 0, 6, 0, 7,
470                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
471                             PRIM_INT__s_s_i_s, NULL, NULL, PRIM_INT__s_s_h_s, NULL, PRIM_INT__p_s_i_s);
472 
473 
474                     ostei_general_vrr_I(2, 0, 7, 0, 6,
475                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
476                             PRIM_INT__p_s_k_s, PRIM_INT__s_s_k_s, NULL, PRIM_INT__p_s_i_s, NULL, PRIM_INT__d_s_k_s);
477 
478 
479                     ostei_general_vrr_I(1, 0, 5, 0, 7,
480                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
481                             PRIM_INT__s_s_h_s, NULL, NULL, PRIM_INT__s_s_g_s, NULL, PRIM_INT__p_s_h_s);
482 
483 
484                     ostei_general_vrr_I(2, 0, 6, 0, 6,
485                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
486                             PRIM_INT__p_s_i_s, PRIM_INT__s_s_i_s, NULL, PRIM_INT__p_s_h_s, NULL, PRIM_INT__d_s_i_s);
487 
488 
489                     ostei_general_vrr_I(3, 0, 7, 0, 5,
490                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
491                             PRIM_INT__d_s_k_s, PRIM_INT__p_s_k_s, NULL, PRIM_INT__d_s_i_s, NULL, PRIM_INT__f_s_k_s);
492 
493 
494                     VRR_I_p_s_g_s(
495                             PRIM_INT__p_s_g_s,
496                             PRIM_INT__s_s_g_s,
497                             PRIM_INT__s_s_f_s,
498                             P_PA,
499                             aop_PQ,
500                             one_over_2pq,
501                             7);
502 
503 
504                     ostei_general_vrr_I(2, 0, 5, 0, 6,
505                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
506                             PRIM_INT__p_s_h_s, PRIM_INT__s_s_h_s, NULL, PRIM_INT__p_s_g_s, NULL, PRIM_INT__d_s_h_s);
507 
508 
509                     ostei_general_vrr_I(3, 0, 6, 0, 5,
510                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
511                             PRIM_INT__d_s_i_s, PRIM_INT__p_s_i_s, NULL, PRIM_INT__d_s_h_s, NULL, PRIM_INT__f_s_i_s);
512 
513 
514                     ostei_general_vrr_I(4, 0, 7, 0, 4,
515                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
516                             PRIM_INT__f_s_k_s, PRIM_INT__d_s_k_s, NULL, PRIM_INT__f_s_i_s, NULL, PRIM_INT__g_s_k_s);
517 
518 
519                     ostei_general_vrr1_K(8, 9,
520                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
521                             PRIM_INT__s_s_k_s, PRIM_INT__s_s_i_s, PRIM_INT__s_s_l_s);
522 
523 
524                     ostei_general_vrr_I(1, 0, 8, 0, 7,
525                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
526                             PRIM_INT__s_s_l_s, NULL, NULL, PRIM_INT__s_s_k_s, NULL, PRIM_INT__p_s_l_s);
527 
528 
529                     ostei_general_vrr_I(2, 0, 8, 0, 6,
530                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
531                             PRIM_INT__p_s_l_s, PRIM_INT__s_s_l_s, NULL, PRIM_INT__p_s_k_s, NULL, PRIM_INT__d_s_l_s);
532 
533 
534                     ostei_general_vrr_I(3, 0, 8, 0, 5,
535                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
536                             PRIM_INT__d_s_l_s, PRIM_INT__p_s_l_s, NULL, PRIM_INT__d_s_k_s, NULL, PRIM_INT__f_s_l_s);
537 
538 
539                     ostei_general_vrr_I(4, 0, 8, 0, 4,
540                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
541                             PRIM_INT__f_s_l_s, PRIM_INT__d_s_l_s, NULL, PRIM_INT__f_s_k_s, NULL, PRIM_INT__g_s_l_s);
542 
543 
544                     VRR_I_p_s_f_s(
545                             PRIM_INT__p_s_f_s,
546                             PRIM_INT__s_s_f_s,
547                             PRIM_INT__s_s_d_s,
548                             P_PA,
549                             aop_PQ,
550                             one_over_2pq,
551                             7);
552 
553 
554                     ostei_general_vrr_I(2, 0, 4, 0, 6,
555                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
556                             PRIM_INT__p_s_g_s, PRIM_INT__s_s_g_s, NULL, PRIM_INT__p_s_f_s, NULL, PRIM_INT__d_s_g_s);
557 
558 
559                     ostei_general_vrr_I(3, 0, 5, 0, 5,
560                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
561                             PRIM_INT__d_s_h_s, PRIM_INT__p_s_h_s, NULL, PRIM_INT__d_s_g_s, NULL, PRIM_INT__f_s_h_s);
562 
563 
564                     ostei_general_vrr_I(4, 0, 6, 0, 4,
565                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
566                             PRIM_INT__f_s_i_s, PRIM_INT__d_s_i_s, NULL, PRIM_INT__f_s_h_s, NULL, PRIM_INT__g_s_i_s);
567 
568 
569                     ostei_general_vrr_I(5, 0, 7, 0, 3,
570                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
571                             PRIM_INT__g_s_k_s, PRIM_INT__f_s_k_s, NULL, PRIM_INT__g_s_i_s, NULL, PRIM_INT__h_s_k_s);
572 
573 
574                     ostei_general_vrr1_K(9, 8,
575                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
576                             PRIM_INT__s_s_l_s, PRIM_INT__s_s_k_s, PRIM_INT__s_s_m_s);
577 
578 
579                     ostei_general_vrr_I(1, 0, 9, 0, 7,
580                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
581                             PRIM_INT__s_s_m_s, NULL, NULL, PRIM_INT__s_s_l_s, NULL, PRIM_INT__p_s_m_s);
582 
583 
584                     ostei_general_vrr_I(2, 0, 9, 0, 6,
585                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
586                             PRIM_INT__p_s_m_s, PRIM_INT__s_s_m_s, NULL, PRIM_INT__p_s_l_s, NULL, PRIM_INT__d_s_m_s);
587 
588 
589                     ostei_general_vrr_I(3, 0, 9, 0, 5,
590                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
591                             PRIM_INT__d_s_m_s, PRIM_INT__p_s_m_s, NULL, PRIM_INT__d_s_l_s, NULL, PRIM_INT__f_s_m_s);
592 
593 
594                     ostei_general_vrr_I(4, 0, 9, 0, 4,
595                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
596                             PRIM_INT__f_s_m_s, PRIM_INT__d_s_m_s, NULL, PRIM_INT__f_s_l_s, NULL, PRIM_INT__g_s_m_s);
597 
598 
599                     ostei_general_vrr_I(5, 0, 8, 0, 3,
600                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
601                             PRIM_INT__g_s_l_s, PRIM_INT__f_s_l_s, NULL, PRIM_INT__g_s_k_s, NULL, PRIM_INT__h_s_l_s);
602 
603 
604 
605                     // Forming PRIM_INT__p_s_d_s[7 * 18];
606                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
607                     {
608 
609                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 0]);
610                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__p_s_d_s[n * 18 + 0]);
611                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_d_s[n * 18 + 0]);
612 
613                         PRIM_INT__p_s_d_s[n * 18 + 1] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 1]);
614                         PRIM_INT__p_s_d_s[n * 18 + 1] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__p_s_d_s[n * 18 + 1]);
615                         PRIM_INT__p_s_d_s[n * 18 + 1] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_d_s[n * 18 + 1]);
616 
617                         PRIM_INT__p_s_d_s[n * 18 + 2] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 2]);
618                         PRIM_INT__p_s_d_s[n * 18 + 2] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 2], PRIM_INT__p_s_d_s[n * 18 + 2]);
619                         PRIM_INT__p_s_d_s[n * 18 + 2] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_d_s[n * 18 + 2]);
620 
621                         PRIM_INT__p_s_d_s[n * 18 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 3]);
622                         PRIM_INT__p_s_d_s[n * 18 + 3] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__p_s_d_s[n * 18 + 3]);
623 
624                         PRIM_INT__p_s_d_s[n * 18 + 4] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 4]);
625                         PRIM_INT__p_s_d_s[n * 18 + 4] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 4], PRIM_INT__p_s_d_s[n * 18 + 4]);
626 
627                         PRIM_INT__p_s_d_s[n * 18 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 5]);
628                         PRIM_INT__p_s_d_s[n * 18 + 5] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__p_s_d_s[n * 18 + 5]);
629 
630                         PRIM_INT__p_s_d_s[n * 18 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 0]);
631                         PRIM_INT__p_s_d_s[n * 18 + 6] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__p_s_d_s[n * 18 + 6]);
632 
633                         PRIM_INT__p_s_d_s[n * 18 + 7] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 1]);
634                         PRIM_INT__p_s_d_s[n * 18 + 7] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__p_s_d_s[n * 18 + 7]);
635                         PRIM_INT__p_s_d_s[n * 18 + 7] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_d_s[n * 18 + 7]);
636 
637                         PRIM_INT__p_s_d_s[n * 18 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 2]);
638                         PRIM_INT__p_s_d_s[n * 18 + 8] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 2], PRIM_INT__p_s_d_s[n * 18 + 8]);
639 
640                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 3]);
641                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__p_s_d_s[n * 18 + 9]);
642                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_d_s[n * 18 + 9]);
643 
644                         PRIM_INT__p_s_d_s[n * 18 + 10] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 4]);
645                         PRIM_INT__p_s_d_s[n * 18 + 10] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 4], PRIM_INT__p_s_d_s[n * 18 + 10]);
646                         PRIM_INT__p_s_d_s[n * 18 + 10] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_d_s[n * 18 + 10]);
647 
648                         PRIM_INT__p_s_d_s[n * 18 + 11] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 5]);
649                         PRIM_INT__p_s_d_s[n * 18 + 11] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__p_s_d_s[n * 18 + 11]);
650 
651                         PRIM_INT__p_s_d_s[n * 18 + 12] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 0]);
652                         PRIM_INT__p_s_d_s[n * 18 + 12] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__p_s_d_s[n * 18 + 12]);
653 
654                         PRIM_INT__p_s_d_s[n * 18 + 13] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 1]);
655                         PRIM_INT__p_s_d_s[n * 18 + 13] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__p_s_d_s[n * 18 + 13]);
656 
657                         PRIM_INT__p_s_d_s[n * 18 + 14] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 2]);
658                         PRIM_INT__p_s_d_s[n * 18 + 14] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 2], PRIM_INT__p_s_d_s[n * 18 + 14]);
659                         PRIM_INT__p_s_d_s[n * 18 + 14] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_d_s[n * 18 + 14]);
660 
661                         PRIM_INT__p_s_d_s[n * 18 + 15] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 3]);
662                         PRIM_INT__p_s_d_s[n * 18 + 15] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__p_s_d_s[n * 18 + 15]);
663 
664                         PRIM_INT__p_s_d_s[n * 18 + 16] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 4]);
665                         PRIM_INT__p_s_d_s[n * 18 + 16] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 4], PRIM_INT__p_s_d_s[n * 18 + 16]);
666                         PRIM_INT__p_s_d_s[n * 18 + 16] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_d_s[n * 18 + 16]);
667 
668                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 5]);
669                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__p_s_d_s[n * 18 + 17]);
670                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_d_s[n * 18 + 17]);
671 
672                     }
673 
674 
675                     VRR_I_d_s_f_s(
676                             PRIM_INT__d_s_f_s,
677                             PRIM_INT__p_s_f_s,
678                             PRIM_INT__s_s_f_s,
679                             PRIM_INT__p_s_d_s,
680                             P_PA,
681                             a_over_p,
682                             aop_PQ,
683                             one_over_2p,
684                             one_over_2pq,
685                             6);
686 
687 
688                     ostei_general_vrr_I(3, 0, 4, 0, 5,
689                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
690                             PRIM_INT__d_s_g_s, PRIM_INT__p_s_g_s, NULL, PRIM_INT__d_s_f_s, NULL, PRIM_INT__f_s_g_s);
691 
692 
693                     ostei_general_vrr_I(4, 0, 5, 0, 4,
694                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
695                             PRIM_INT__f_s_h_s, PRIM_INT__d_s_h_s, NULL, PRIM_INT__f_s_g_s, NULL, PRIM_INT__g_s_h_s);
696 
697 
698                     ostei_general_vrr_I(5, 0, 6, 0, 3,
699                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
700                             PRIM_INT__g_s_i_s, PRIM_INT__f_s_i_s, NULL, PRIM_INT__g_s_h_s, NULL, PRIM_INT__h_s_i_s);
701 
702 
703                     ostei_general_vrr_I(6, 0, 7, 0, 2,
704                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
705                             PRIM_INT__h_s_k_s, PRIM_INT__g_s_k_s, NULL, PRIM_INT__h_s_i_s, NULL, PRIM_INT__i_s_k_s);
706 
707 
708                     ostei_general_vrr_I(5, 0, 9, 0, 3,
709                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
710                             PRIM_INT__g_s_m_s, PRIM_INT__f_s_m_s, NULL, PRIM_INT__g_s_l_s, NULL, PRIM_INT__h_s_m_s);
711 
712 
713                     ostei_general_vrr_I(6, 0, 8, 0, 2,
714                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
715                             PRIM_INT__h_s_l_s, PRIM_INT__g_s_l_s, NULL, PRIM_INT__h_s_k_s, NULL, PRIM_INT__i_s_l_s);
716 
717 
718 
719                     // Forming PRIM_INT__p_s_p_s[7 * 9];
720                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
721                     {
722 
723                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_p_s[n * 3 + 0]);
724                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 0]);
725                         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]);
726 
727                         PRIM_INT__p_s_p_s[n * 9 + 1] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_p_s[n * 3 + 1]);
728                         PRIM_INT__p_s_p_s[n * 9 + 1] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 1]);
729 
730                         PRIM_INT__p_s_p_s[n * 9 + 2] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_p_s[n * 3 + 2]);
731                         PRIM_INT__p_s_p_s[n * 9 + 2] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 2]);
732 
733                         PRIM_INT__p_s_p_s[n * 9 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_p_s[n * 3 + 0]);
734                         PRIM_INT__p_s_p_s[n * 9 + 3] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 3]);
735 
736                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_p_s[n * 3 + 1]);
737                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 4]);
738                         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]);
739 
740                         PRIM_INT__p_s_p_s[n * 9 + 5] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_p_s[n * 3 + 2]);
741                         PRIM_INT__p_s_p_s[n * 9 + 5] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 5]);
742 
743                         PRIM_INT__p_s_p_s[n * 9 + 6] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_p_s[n * 3 + 0]);
744                         PRIM_INT__p_s_p_s[n * 9 + 6] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 6]);
745 
746                         PRIM_INT__p_s_p_s[n * 9 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_p_s[n * 3 + 1]);
747                         PRIM_INT__p_s_p_s[n * 9 + 7] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 7]);
748 
749                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_p_s[n * 3 + 2]);
750                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 8]);
751                         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]);
752 
753                     }
754 
755 
756                     VRR_I_d_s_d_s(
757                             PRIM_INT__d_s_d_s,
758                             PRIM_INT__p_s_d_s,
759                             PRIM_INT__s_s_d_s,
760                             PRIM_INT__p_s_p_s,
761                             P_PA,
762                             a_over_p,
763                             aop_PQ,
764                             one_over_2p,
765                             one_over_2pq,
766                             6);
767 
768 
769                     ostei_general_vrr_I(3, 0, 3, 0, 5,
770                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
771                             PRIM_INT__d_s_f_s, PRIM_INT__p_s_f_s, NULL, PRIM_INT__d_s_d_s, NULL, PRIM_INT__f_s_f_s);
772 
773 
774                     ostei_general_vrr_I(4, 0, 4, 0, 4,
775                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
776                             PRIM_INT__f_s_g_s, PRIM_INT__d_s_g_s, NULL, PRIM_INT__f_s_f_s, NULL, PRIM_INT__g_s_g_s);
777 
778 
779                     ostei_general_vrr_I(5, 0, 5, 0, 3,
780                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
781                             PRIM_INT__g_s_h_s, PRIM_INT__f_s_h_s, NULL, PRIM_INT__g_s_g_s, NULL, PRIM_INT__h_s_h_s);
782 
783 
784                     ostei_general_vrr_I(6, 0, 6, 0, 2,
785                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
786                             PRIM_INT__h_s_i_s, PRIM_INT__g_s_i_s, NULL, PRIM_INT__h_s_h_s, NULL, PRIM_INT__i_s_i_s);
787 
788 
789                     ostei_general_vrr_I(7, 0, 7, 0, 1,
790                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
791                             PRIM_INT__i_s_k_s, PRIM_INT__h_s_k_s, NULL, PRIM_INT__i_s_i_s, NULL, PRIM_INT__k_s_k_s);
792 
793 
794                     ostei_general_vrr_I(6, 0, 9, 0, 2,
795                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
796                             PRIM_INT__h_s_m_s, PRIM_INT__g_s_m_s, NULL, PRIM_INT__h_s_l_s, NULL, PRIM_INT__i_s_m_s);
797 
798 
799                     ostei_general_vrr_I(7, 0, 8, 0, 1,
800                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
801                             PRIM_INT__i_s_l_s, PRIM_INT__h_s_l_s, NULL, PRIM_INT__i_s_k_s, NULL, PRIM_INT__k_s_l_s);
802 
803 
804                     ostei_general_vrr_I(7, 0, 9, 0, 1,
805                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
806                             PRIM_INT__i_s_m_s, PRIM_INT__h_s_m_s, NULL, PRIM_INT__i_s_l_s, NULL, PRIM_INT__k_s_m_s);
807 
808 
809 
810 
811                     ////////////////////////////////////
812                     // Accumulate contracted integrals
813                     ////////////////////////////////////
814                     if(lastoffset == 0)
815                     {
816                         contract_all(540, PRIM_INT__g_s_k_s, PRIM_PTR_INT__g_s_k_s);
817                         contract_all(675, PRIM_INT__g_s_l_s, PRIM_PTR_INT__g_s_l_s);
818                         contract_all(825, PRIM_INT__g_s_m_s, PRIM_PTR_INT__g_s_m_s);
819                         contract_all(756, PRIM_INT__h_s_k_s, PRIM_PTR_INT__h_s_k_s);
820                         contract_all(945, PRIM_INT__h_s_l_s, PRIM_PTR_INT__h_s_l_s);
821                         contract_all(1155, PRIM_INT__h_s_m_s, PRIM_PTR_INT__h_s_m_s);
822                         contract_all(1008, PRIM_INT__i_s_k_s, PRIM_PTR_INT__i_s_k_s);
823                         contract_all(1260, PRIM_INT__i_s_l_s, PRIM_PTR_INT__i_s_l_s);
824                         contract_all(1540, PRIM_INT__i_s_m_s, PRIM_PTR_INT__i_s_m_s);
825                         contract_all(1296, PRIM_INT__k_s_k_s, PRIM_PTR_INT__k_s_k_s);
826                         contract_all(1620, PRIM_INT__k_s_l_s, PRIM_PTR_INT__k_s_l_s);
827                         contract_all(1980, PRIM_INT__k_s_m_s, PRIM_PTR_INT__k_s_m_s);
828                     }
829                     else
830                     {
831                         contract(540, shelloffsets, PRIM_INT__g_s_k_s, PRIM_PTR_INT__g_s_k_s);
832                         contract(675, shelloffsets, PRIM_INT__g_s_l_s, PRIM_PTR_INT__g_s_l_s);
833                         contract(825, shelloffsets, PRIM_INT__g_s_m_s, PRIM_PTR_INT__g_s_m_s);
834                         contract(756, shelloffsets, PRIM_INT__h_s_k_s, PRIM_PTR_INT__h_s_k_s);
835                         contract(945, shelloffsets, PRIM_INT__h_s_l_s, PRIM_PTR_INT__h_s_l_s);
836                         contract(1155, shelloffsets, PRIM_INT__h_s_m_s, PRIM_PTR_INT__h_s_m_s);
837                         contract(1008, shelloffsets, PRIM_INT__i_s_k_s, PRIM_PTR_INT__i_s_k_s);
838                         contract(1260, shelloffsets, PRIM_INT__i_s_l_s, PRIM_PTR_INT__i_s_l_s);
839                         contract(1540, shelloffsets, PRIM_INT__i_s_m_s, PRIM_PTR_INT__i_s_m_s);
840                         contract(1296, shelloffsets, PRIM_INT__k_s_k_s, PRIM_PTR_INT__k_s_k_s);
841                         contract(1620, shelloffsets, PRIM_INT__k_s_l_s, PRIM_PTR_INT__k_s_l_s);
842                         contract(1980, shelloffsets, PRIM_INT__k_s_m_s, PRIM_PTR_INT__k_s_m_s);
843                         PRIM_PTR_INT__g_s_k_s += lastoffset*540;
844                         PRIM_PTR_INT__g_s_l_s += lastoffset*675;
845                         PRIM_PTR_INT__g_s_m_s += lastoffset*825;
846                         PRIM_PTR_INT__h_s_k_s += lastoffset*756;
847                         PRIM_PTR_INT__h_s_l_s += lastoffset*945;
848                         PRIM_PTR_INT__h_s_m_s += lastoffset*1155;
849                         PRIM_PTR_INT__i_s_k_s += lastoffset*1008;
850                         PRIM_PTR_INT__i_s_l_s += lastoffset*1260;
851                         PRIM_PTR_INT__i_s_m_s += lastoffset*1540;
852                         PRIM_PTR_INT__k_s_k_s += lastoffset*1296;
853                         PRIM_PTR_INT__k_s_l_s += lastoffset*1620;
854                         PRIM_PTR_INT__k_s_m_s += lastoffset*1980;
855                     }
856 
857                 }  // close loop over j
858             }  // close loop over i
859 
860             //Advance to the next batch
861             jstart = SIMINT_SIMD_ROUND(jend);
862 
863             //////////////////////////////////////////////
864             // Contracted integrals: Horizontal recurrance
865             //////////////////////////////////////////////
866 
867 
868             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
869 
870 
871             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
872             {
873                 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
874 
875                 // set up HRR pointers
876                 double const * restrict HRR_INT__g_s_k_s = INT__g_s_k_s + abcd * 540;
877                 double const * restrict HRR_INT__g_s_l_s = INT__g_s_l_s + abcd * 675;
878                 double const * restrict HRR_INT__g_s_m_s = INT__g_s_m_s + abcd * 825;
879                 double const * restrict HRR_INT__h_s_k_s = INT__h_s_k_s + abcd * 756;
880                 double const * restrict HRR_INT__h_s_l_s = INT__h_s_l_s + abcd * 945;
881                 double const * restrict HRR_INT__h_s_m_s = INT__h_s_m_s + abcd * 1155;
882                 double const * restrict HRR_INT__i_s_k_s = INT__i_s_k_s + abcd * 1008;
883                 double const * restrict HRR_INT__i_s_l_s = INT__i_s_l_s + abcd * 1260;
884                 double const * restrict HRR_INT__i_s_m_s = INT__i_s_m_s + abcd * 1540;
885                 double const * restrict HRR_INT__k_s_k_s = INT__k_s_k_s + abcd * 1296;
886                 double const * restrict HRR_INT__k_s_l_s = INT__k_s_l_s + abcd * 1620;
887                 double const * restrict HRR_INT__k_s_m_s = INT__k_s_m_s + abcd * 1980;
888                 double * restrict HRR_INT__g_f_k_d = INT__g_f_k_d + real_abcd * 32400;
889 
890                 // form INT__g_p_k_s
891                 HRR_J_g_p(
892                     HRR_INT__g_p_k_s,
893                     HRR_INT__g_s_k_s,
894                     HRR_INT__h_s_k_s,
895                     hAB, 36);
896 
897                 // form INT__g_p_l_s
898                 HRR_J_g_p(
899                     HRR_INT__g_p_l_s,
900                     HRR_INT__g_s_l_s,
901                     HRR_INT__h_s_l_s,
902                     hAB, 45);
903 
904                 // form INT__g_p_m_s
905                 HRR_J_g_p(
906                     HRR_INT__g_p_m_s,
907                     HRR_INT__g_s_m_s,
908                     HRR_INT__h_s_m_s,
909                     hAB, 55);
910 
911                 // form INT__h_p_k_s
912                 ostei_general_hrr_J(5, 1, 7, 0, hAB, HRR_INT__i_s_k_s, HRR_INT__h_s_k_s, HRR_INT__h_p_k_s);
913 
914                 // form INT__h_p_l_s
915                 ostei_general_hrr_J(5, 1, 8, 0, hAB, HRR_INT__i_s_l_s, HRR_INT__h_s_l_s, HRR_INT__h_p_l_s);
916 
917                 // form INT__h_p_m_s
918                 ostei_general_hrr_J(5, 1, 9, 0, hAB, HRR_INT__i_s_m_s, HRR_INT__h_s_m_s, HRR_INT__h_p_m_s);
919 
920                 // form INT__i_p_k_s
921                 ostei_general_hrr_J(6, 1, 7, 0, hAB, HRR_INT__k_s_k_s, HRR_INT__i_s_k_s, HRR_INT__i_p_k_s);
922 
923                 // form INT__i_p_l_s
924                 ostei_general_hrr_J(6, 1, 8, 0, hAB, HRR_INT__k_s_l_s, HRR_INT__i_s_l_s, HRR_INT__i_p_l_s);
925 
926                 // form INT__i_p_m_s
927                 ostei_general_hrr_J(6, 1, 9, 0, hAB, HRR_INT__k_s_m_s, HRR_INT__i_s_m_s, HRR_INT__i_p_m_s);
928 
929                 // form INT__g_d_k_s
930                 ostei_general_hrr_J(4, 2, 7, 0, hAB, HRR_INT__h_p_k_s, HRR_INT__g_p_k_s, HRR_INT__g_d_k_s);
931 
932                 // form INT__g_d_l_s
933                 ostei_general_hrr_J(4, 2, 8, 0, hAB, HRR_INT__h_p_l_s, HRR_INT__g_p_l_s, HRR_INT__g_d_l_s);
934 
935                 // form INT__g_d_m_s
936                 ostei_general_hrr_J(4, 2, 9, 0, hAB, HRR_INT__h_p_m_s, HRR_INT__g_p_m_s, HRR_INT__g_d_m_s);
937 
938                 // form INT__h_d_k_s
939                 ostei_general_hrr_J(5, 2, 7, 0, hAB, HRR_INT__i_p_k_s, HRR_INT__h_p_k_s, HRR_INT__h_d_k_s);
940 
941                 // form INT__h_d_l_s
942                 ostei_general_hrr_J(5, 2, 8, 0, hAB, HRR_INT__i_p_l_s, HRR_INT__h_p_l_s, HRR_INT__h_d_l_s);
943 
944                 // form INT__h_d_m_s
945                 ostei_general_hrr_J(5, 2, 9, 0, hAB, HRR_INT__i_p_m_s, HRR_INT__h_p_m_s, HRR_INT__h_d_m_s);
946 
947                 // form INT__g_f_k_s
948                 ostei_general_hrr_J(4, 3, 7, 0, hAB, HRR_INT__h_d_k_s, HRR_INT__g_d_k_s, HRR_INT__g_f_k_s);
949 
950                 // form INT__g_f_l_s
951                 ostei_general_hrr_J(4, 3, 8, 0, hAB, HRR_INT__h_d_l_s, HRR_INT__g_d_l_s, HRR_INT__g_f_l_s);
952 
953                 // form INT__g_f_m_s
954                 ostei_general_hrr_J(4, 3, 9, 0, hAB, HRR_INT__h_d_m_s, HRR_INT__g_d_m_s, HRR_INT__g_f_m_s);
955 
956                 // form INT__g_f_k_p
957                 ostei_general_hrr_L(4, 3, 7, 1, hCD, HRR_INT__g_f_l_s, HRR_INT__g_f_k_s, HRR_INT__g_f_k_p);
958 
959                 // form INT__g_f_l_p
960                 ostei_general_hrr_L(4, 3, 8, 1, hCD, HRR_INT__g_f_m_s, HRR_INT__g_f_l_s, HRR_INT__g_f_l_p);
961 
962                 // form INT__g_f_k_d
963                 ostei_general_hrr_L(4, 3, 7, 2, hCD, HRR_INT__g_f_l_p, HRR_INT__g_f_k_p, HRR_INT__g_f_k_d);
964 
965 
966             }  // close HRR loop
967 
968 
969         }   // close loop cdbatch
970 
971         istart = iend;
972     }  // close loop over ab
973 
974     return P.nshell12_clip * Q.nshell12_clip;
975 }
976 
ostei_f_g_k_d(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__f_g_k_d)977 int ostei_f_g_k_d(struct simint_multi_shellpair const P,
978                   struct simint_multi_shellpair const Q,
979                   double screen_tol,
980                   double * const restrict work,
981                   double * const restrict INT__f_g_k_d)
982 {
983     double P_AB[3*P.nshell12];
984     struct simint_multi_shellpair P_tmp = P;
985     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
986     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
987     P_tmp.AB_x = P_AB;
988     P_tmp.AB_y = P_AB + P.nshell12;
989     P_tmp.AB_z = P_AB + 2*P.nshell12;
990 
991     for(int i = 0; i < P.nshell12; i++)
992     {
993         P_tmp.AB_x[i] = -P.AB_x[i];
994         P_tmp.AB_y[i] = -P.AB_y[i];
995         P_tmp.AB_z[i] = -P.AB_z[i];
996     }
997 
998     int ret = ostei_g_f_k_d(P_tmp, Q, screen_tol, work, INT__f_g_k_d);
999     double buffer[32400] SIMINT_ALIGN_ARRAY_DBL;
1000 
1001     for(int q = 0; q < ret; q++)
1002     {
1003         int idx = 0;
1004         for(int a = 0; a < 10; ++a)
1005         for(int b = 0; b < 15; ++b)
1006         for(int c = 0; c < 36; ++c)
1007         for(int d = 0; d < 6; ++d)
1008             buffer[idx++] = INT__f_g_k_d[q*32400+b*2160+a*216+c*6+d];
1009 
1010         memcpy(INT__f_g_k_d+q*32400, buffer, 32400*sizeof(double));
1011     }
1012 
1013     return ret;
1014 }
1015 
ostei_g_f_d_k(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__g_f_d_k)1016 int ostei_g_f_d_k(struct simint_multi_shellpair const P,
1017                   struct simint_multi_shellpair const Q,
1018                   double screen_tol,
1019                   double * const restrict work,
1020                   double * const restrict INT__g_f_d_k)
1021 {
1022     double Q_AB[3*Q.nshell12];
1023     struct simint_multi_shellpair Q_tmp = Q;
1024     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
1025     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
1026     Q_tmp.AB_x = Q_AB;
1027     Q_tmp.AB_y = Q_AB + Q.nshell12;
1028     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
1029 
1030     for(int i = 0; i < Q.nshell12; i++)
1031     {
1032         Q_tmp.AB_x[i] = -Q.AB_x[i];
1033         Q_tmp.AB_y[i] = -Q.AB_y[i];
1034         Q_tmp.AB_z[i] = -Q.AB_z[i];
1035     }
1036 
1037     int ret = ostei_g_f_k_d(P, Q_tmp, screen_tol, work, INT__g_f_d_k);
1038     double buffer[32400] SIMINT_ALIGN_ARRAY_DBL;
1039 
1040     for(int q = 0; q < ret; q++)
1041     {
1042         int idx = 0;
1043         for(int a = 0; a < 15; ++a)
1044         for(int b = 0; b < 10; ++b)
1045         for(int c = 0; c < 6; ++c)
1046         for(int d = 0; d < 36; ++d)
1047             buffer[idx++] = INT__g_f_d_k[q*32400+a*2160+b*216+d*6+c];
1048 
1049         memcpy(INT__g_f_d_k+q*32400, buffer, 32400*sizeof(double));
1050     }
1051 
1052     return ret;
1053 }
1054 
ostei_f_g_d_k(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__f_g_d_k)1055 int ostei_f_g_d_k(struct simint_multi_shellpair const P,
1056                   struct simint_multi_shellpair const Q,
1057                   double screen_tol,
1058                   double * const restrict work,
1059                   double * const restrict INT__f_g_d_k)
1060 {
1061     double P_AB[3*P.nshell12];
1062     struct simint_multi_shellpair P_tmp = P;
1063     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
1064     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
1065     P_tmp.AB_x = P_AB;
1066     P_tmp.AB_y = P_AB + P.nshell12;
1067     P_tmp.AB_z = P_AB + 2*P.nshell12;
1068 
1069     for(int i = 0; i < P.nshell12; i++)
1070     {
1071         P_tmp.AB_x[i] = -P.AB_x[i];
1072         P_tmp.AB_y[i] = -P.AB_y[i];
1073         P_tmp.AB_z[i] = -P.AB_z[i];
1074     }
1075 
1076     double Q_AB[3*Q.nshell12];
1077     struct simint_multi_shellpair Q_tmp = Q;
1078     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
1079     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
1080     Q_tmp.AB_x = Q_AB;
1081     Q_tmp.AB_y = Q_AB + Q.nshell12;
1082     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
1083 
1084     for(int i = 0; i < Q.nshell12; i++)
1085     {
1086         Q_tmp.AB_x[i] = -Q.AB_x[i];
1087         Q_tmp.AB_y[i] = -Q.AB_y[i];
1088         Q_tmp.AB_z[i] = -Q.AB_z[i];
1089     }
1090 
1091     int ret = ostei_g_f_k_d(P_tmp, Q_tmp, screen_tol, work, INT__f_g_d_k);
1092     double buffer[32400] SIMINT_ALIGN_ARRAY_DBL;
1093 
1094     for(int q = 0; q < ret; q++)
1095     {
1096         int idx = 0;
1097         for(int a = 0; a < 10; ++a)
1098         for(int b = 0; b < 15; ++b)
1099         for(int c = 0; c < 6; ++c)
1100         for(int d = 0; d < 36; ++d)
1101             buffer[idx++] = INT__f_g_d_k[q*32400+b*2160+a*216+d*6+c];
1102 
1103         memcpy(INT__f_g_d_k+q*32400, buffer, 32400*sizeof(double));
1104     }
1105 
1106     return ret;
1107 }
1108 
1109