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_h_g_f_d(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__h_g_f_d)8 int ostei_h_g_f_d(struct simint_multi_shellpair const P,
9                   struct simint_multi_shellpair const Q,
10                   double screen_tol,
11                   double * const restrict work,
12                   double * const restrict INT__h_g_f_d)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__h_g_f_d);
17     int ab, cd, abcd;
18     int istart, jstart;
19     int iprimcd, nprim_icd, icd;
20     const int check_screen = (screen_tol > 0.0);
21     int i, j;
22     int n;
23     int not_screened;
24     int real_abcd;
25     int iket;
26     int ibra;
27 
28     // partition workspace
29     double * const INT__h_s_f_s = work + (SIMINT_NSHELL_SIMD * 0);
30     double * const INT__h_s_g_s = work + (SIMINT_NSHELL_SIMD * 210);
31     double * const INT__h_s_h_s = work + (SIMINT_NSHELL_SIMD * 525);
32     double * const INT__i_s_f_s = work + (SIMINT_NSHELL_SIMD * 966);
33     double * const INT__i_s_g_s = work + (SIMINT_NSHELL_SIMD * 1246);
34     double * const INT__i_s_h_s = work + (SIMINT_NSHELL_SIMD * 1666);
35     double * const INT__k_s_f_s = work + (SIMINT_NSHELL_SIMD * 2254);
36     double * const INT__k_s_g_s = work + (SIMINT_NSHELL_SIMD * 2614);
37     double * const INT__k_s_h_s = work + (SIMINT_NSHELL_SIMD * 3154);
38     double * const INT__l_s_f_s = work + (SIMINT_NSHELL_SIMD * 3910);
39     double * const INT__l_s_g_s = work + (SIMINT_NSHELL_SIMD * 4360);
40     double * const INT__l_s_h_s = work + (SIMINT_NSHELL_SIMD * 5035);
41     double * const INT__m_s_f_s = work + (SIMINT_NSHELL_SIMD * 5980);
42     double * const INT__m_s_g_s = work + (SIMINT_NSHELL_SIMD * 6530);
43     double * const INT__m_s_h_s = work + (SIMINT_NSHELL_SIMD * 7355);
44     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*8510);
45     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
46     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_s = primwork + 15;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_p_s = primwork + 57;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 102;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_p_s = primwork + 180;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_d_s = primwork + 270;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_s = primwork + 414;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_p_s = primwork + 534;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_d_s = primwork + 684;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_f_s = primwork + 924;
55     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_s_s = primwork + 1224;
56     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_p_s = primwork + 1389;
57     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_d_s = primwork + 1614;
58     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_f_s = primwork + 1974;
59     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_g_s = primwork + 2424;
60     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_s_s = primwork + 2874;
61     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_p_s = primwork + 3084;
62     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_d_s = primwork + 3399;
63     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_f_s = primwork + 3903;
64     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_g_s = primwork + 4533;
65     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_h_s = primwork + 5163;
66     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_s_s = primwork + 5604;
67     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_p_s = primwork + 5856;
68     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_d_s = primwork + 6276;
69     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_f_s = primwork + 6948;
70     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_g_s = primwork + 7788;
71     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_h_s = primwork + 8628;
72     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_s_s = primwork + 9216;
73     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_p_s = primwork + 9504;
74     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_d_s = primwork + 10044;
75     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_f_s = primwork + 10908;
76     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_g_s = primwork + 11988;
77     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_h_s = primwork + 13068;
78     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_s_s = primwork + 13824;
79     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_p_s = primwork + 14139;
80     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_d_s = primwork + 14814;
81     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_f_s = primwork + 15894;
82     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_g_s = primwork + 17244;
83     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_h_s = primwork + 18594;
84     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_s_s = primwork + 19539;
85     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_p_s = primwork + 19869;
86     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_d_s = primwork + 20694;
87     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_f_s = primwork + 22014;
88     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_g_s = primwork + 23664;
89     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_h_s = primwork + 25314;
90     double * const hrrwork = (double *)(primwork + 26469);
91     double * const HRR_INT__h_p_f_s = hrrwork + 0;
92     double * const HRR_INT__h_p_g_s = hrrwork + 630;
93     double * const HRR_INT__h_p_h_s = hrrwork + 1575;
94     double * const HRR_INT__h_d_f_s = hrrwork + 2898;
95     double * const HRR_INT__h_d_g_s = hrrwork + 4158;
96     double * const HRR_INT__h_d_h_s = hrrwork + 6048;
97     double * const HRR_INT__h_f_f_s = hrrwork + 8694;
98     double * const HRR_INT__h_f_g_s = hrrwork + 10794;
99     double * const HRR_INT__h_f_h_s = hrrwork + 13944;
100     double * const HRR_INT__h_g_f_s = hrrwork + 18354;
101     double * const HRR_INT__h_g_f_p = hrrwork + 21504;
102     double * const HRR_INT__h_g_g_s = hrrwork + 30954;
103     double * const HRR_INT__h_g_g_p = hrrwork + 35679;
104     double * const HRR_INT__h_g_h_s = hrrwork + 49854;
105     double * const HRR_INT__i_p_f_s = hrrwork + 56469;
106     double * const HRR_INT__i_p_g_s = hrrwork + 57309;
107     double * const HRR_INT__i_p_h_s = hrrwork + 58569;
108     double * const HRR_INT__i_d_f_s = hrrwork + 60333;
109     double * const HRR_INT__i_d_g_s = hrrwork + 62013;
110     double * const HRR_INT__i_d_h_s = hrrwork + 64533;
111     double * const HRR_INT__i_f_f_s = hrrwork + 68061;
112     double * const HRR_INT__i_f_g_s = hrrwork + 70861;
113     double * const HRR_INT__i_f_h_s = hrrwork + 75061;
114     double * const HRR_INT__k_p_f_s = hrrwork + 80941;
115     double * const HRR_INT__k_p_g_s = hrrwork + 82021;
116     double * const HRR_INT__k_p_h_s = hrrwork + 83641;
117     double * const HRR_INT__k_d_f_s = hrrwork + 85909;
118     double * const HRR_INT__k_d_g_s = hrrwork + 88069;
119     double * const HRR_INT__k_d_h_s = hrrwork + 91309;
120     double * const HRR_INT__l_p_f_s = hrrwork + 95845;
121     double * const HRR_INT__l_p_g_s = hrrwork + 97195;
122     double * const HRR_INT__l_p_h_s = hrrwork + 99220;
123 
124 
125     // Create constants
126     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
127     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
128     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
129     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
130     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
131     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
132     const SIMINT_DBLTYPE const_7 = SIMINT_DBLSET1(7);
133     const SIMINT_DBLTYPE const_8 = SIMINT_DBLSET1(8);
134     const SIMINT_DBLTYPE const_9 = SIMINT_DBLSET1(9);
135     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
136 
137 
138     ////////////////////////////////////////
139     // Loop over shells and primitives
140     ////////////////////////////////////////
141 
142     real_abcd = 0;
143     istart = 0;
144     for(ab = 0; ab < P.nshell12_clip; ++ab)
145     {
146         const int iend = istart + P.nprim12[ab];
147 
148         cd = 0;
149         jstart = 0;
150 
151         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
152         {
153             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
154             int jend = jstart;
155             for(i = 0; i < nshellbatch; i++)
156                 jend += Q.nprim12[cd+i];
157 
158             // Clear the beginning of the workspace (where we are accumulating integrals)
159             memset(work, 0, SIMINT_NSHELL_SIMD * 8510 * sizeof(double));
160             abcd = 0;
161 
162 
163             for(i = istart; i < iend; ++i)
164             {
165                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
166 
167                 if(check_screen)
168                 {
169                     // Skip this whole thing if always insignificant
170                     if((P.screen[i] * Q.screen_max) < screen_tol)
171                         continue;
172                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
173                 }
174 
175                 icd = 0;
176                 iprimcd = 0;
177                 nprim_icd = Q.nprim12[cd];
178                 double * restrict PRIM_PTR_INT__h_s_f_s = INT__h_s_f_s + abcd * 210;
179                 double * restrict PRIM_PTR_INT__h_s_g_s = INT__h_s_g_s + abcd * 315;
180                 double * restrict PRIM_PTR_INT__h_s_h_s = INT__h_s_h_s + abcd * 441;
181                 double * restrict PRIM_PTR_INT__i_s_f_s = INT__i_s_f_s + abcd * 280;
182                 double * restrict PRIM_PTR_INT__i_s_g_s = INT__i_s_g_s + abcd * 420;
183                 double * restrict PRIM_PTR_INT__i_s_h_s = INT__i_s_h_s + abcd * 588;
184                 double * restrict PRIM_PTR_INT__k_s_f_s = INT__k_s_f_s + abcd * 360;
185                 double * restrict PRIM_PTR_INT__k_s_g_s = INT__k_s_g_s + abcd * 540;
186                 double * restrict PRIM_PTR_INT__k_s_h_s = INT__k_s_h_s + abcd * 756;
187                 double * restrict PRIM_PTR_INT__l_s_f_s = INT__l_s_f_s + abcd * 450;
188                 double * restrict PRIM_PTR_INT__l_s_g_s = INT__l_s_g_s + abcd * 675;
189                 double * restrict PRIM_PTR_INT__l_s_h_s = INT__l_s_h_s + abcd * 945;
190                 double * restrict PRIM_PTR_INT__m_s_f_s = INT__m_s_f_s + abcd * 550;
191                 double * restrict PRIM_PTR_INT__m_s_g_s = INT__m_s_g_s + abcd * 825;
192                 double * restrict PRIM_PTR_INT__m_s_h_s = INT__m_s_h_s + abcd * 1155;
193 
194 
195 
196                 // Load these one per loop over i
197                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
198                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
199                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
200 
201                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
202 
203                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
204                 {
205                     // calculate the shell offsets
206                     // these are the offset from the shell pointed to by cd
207                     // for each element
208                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
209                     int lastoffset = 0;
210                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
211 
212                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
213                     {
214                         // Handle if the first element of the vector is a new shell
215                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
216                         {
217                             nprim_icd += Q.nprim12[cd + (++icd)];
218                             PRIM_PTR_INT__h_s_f_s += 210;
219                             PRIM_PTR_INT__h_s_g_s += 315;
220                             PRIM_PTR_INT__h_s_h_s += 441;
221                             PRIM_PTR_INT__i_s_f_s += 280;
222                             PRIM_PTR_INT__i_s_g_s += 420;
223                             PRIM_PTR_INT__i_s_h_s += 588;
224                             PRIM_PTR_INT__k_s_f_s += 360;
225                             PRIM_PTR_INT__k_s_g_s += 540;
226                             PRIM_PTR_INT__k_s_h_s += 756;
227                             PRIM_PTR_INT__l_s_f_s += 450;
228                             PRIM_PTR_INT__l_s_g_s += 675;
229                             PRIM_PTR_INT__l_s_h_s += 945;
230                             PRIM_PTR_INT__m_s_f_s += 550;
231                             PRIM_PTR_INT__m_s_g_s += 825;
232                             PRIM_PTR_INT__m_s_h_s += 1155;
233                         }
234                         iprimcd++;
235                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
236                         {
237                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
238                             {
239                                 shelloffsets[n] = shelloffsets[n-1] + 1;
240                                 lastoffset++;
241                                 nprim_icd += Q.nprim12[cd + (++icd)];
242                             }
243                             else
244                                 shelloffsets[n] = shelloffsets[n-1];
245                             iprimcd++;
246                         }
247                     }
248                     else
249                         iprimcd += SIMINT_SIMD_LEN;
250 
251                     // Do we have to compute this vector (or has it been screened out)?
252                     // (not_screened != 0 means we have to do this vector)
253                     if(check_screen)
254                     {
255                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
256                         if(vmax < screen_tol)
257                         {
258                             PRIM_PTR_INT__h_s_f_s += lastoffset*210;
259                             PRIM_PTR_INT__h_s_g_s += lastoffset*315;
260                             PRIM_PTR_INT__h_s_h_s += lastoffset*441;
261                             PRIM_PTR_INT__i_s_f_s += lastoffset*280;
262                             PRIM_PTR_INT__i_s_g_s += lastoffset*420;
263                             PRIM_PTR_INT__i_s_h_s += lastoffset*588;
264                             PRIM_PTR_INT__k_s_f_s += lastoffset*360;
265                             PRIM_PTR_INT__k_s_g_s += lastoffset*540;
266                             PRIM_PTR_INT__k_s_h_s += lastoffset*756;
267                             PRIM_PTR_INT__l_s_f_s += lastoffset*450;
268                             PRIM_PTR_INT__l_s_g_s += lastoffset*675;
269                             PRIM_PTR_INT__l_s_h_s += lastoffset*945;
270                             PRIM_PTR_INT__m_s_f_s += lastoffset*550;
271                             PRIM_PTR_INT__m_s_g_s += lastoffset*825;
272                             PRIM_PTR_INT__m_s_h_s += lastoffset*1155;
273                             continue;
274                         }
275                     }
276 
277                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
278                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
279                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
280                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
281 
282 
283                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
284                     SIMINT_DBLTYPE PQ[3];
285                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
286                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
287                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
288                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
289                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
290                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
291 
292                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
293                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
294                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
295                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
296                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
297                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
298                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
299 
300                     // NOTE: Minus sign!
301                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
302                     SIMINT_DBLTYPE aop_PQ[3];
303                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
304                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
305                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
306 
307                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
308                     SIMINT_DBLTYPE aoq_PQ[3];
309                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
310                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
311                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
312                     // Put a minus sign here so we don't have to in RR routines
313                     a_over_q = SIMINT_NEG(a_over_q);
314 
315 
316                     //////////////////////////////////////////////
317                     // Fjt function section
318                     // Maximum v value: 14
319                     //////////////////////////////////////////////
320                     // The parameter to the Fjt function
321                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
322 
323 
324                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
325 
326 
327                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 14);
328                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
329                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
330                     for(n = 0; n <= 14; n++)
331                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
332 
333                     //////////////////////////////////////////////
334                     // Primitive integrals: Vertical recurrance
335                     //////////////////////////////////////////////
336 
337                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
338                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
339                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
340                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
341                     const SIMINT_DBLTYPE vrr_const_5_over_2p = SIMINT_MUL(const_5, one_over_2p);
342                     const SIMINT_DBLTYPE vrr_const_6_over_2p = SIMINT_MUL(const_6, one_over_2p);
343                     const SIMINT_DBLTYPE vrr_const_7_over_2p = SIMINT_MUL(const_7, one_over_2p);
344                     const SIMINT_DBLTYPE vrr_const_8_over_2p = SIMINT_MUL(const_8, one_over_2p);
345                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
346                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
347                     const SIMINT_DBLTYPE vrr_const_3_over_2q = SIMINT_MUL(const_3, one_over_2q);
348                     const SIMINT_DBLTYPE vrr_const_4_over_2q = SIMINT_MUL(const_4, one_over_2q);
349                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
350                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
351                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
352                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
353                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
354                     const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
355                     const SIMINT_DBLTYPE vrr_const_7_over_2pq = SIMINT_MUL(const_7, one_over_2pq);
356                     const SIMINT_DBLTYPE vrr_const_8_over_2pq = SIMINT_MUL(const_8, one_over_2pq);
357                     const SIMINT_DBLTYPE vrr_const_9_over_2pq = SIMINT_MUL(const_9, one_over_2pq);
358 
359 
360 
361                     // Forming PRIM_INT__p_s_s_s[14 * 3];
362                     for(n = 0; n < 14; ++n)  // loop over orders of auxiliary function
363                     {
364 
365                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
366                         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]);
367 
368                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
369                         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]);
370 
371                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
372                         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]);
373 
374                     }
375 
376 
377 
378                     // Forming PRIM_INT__d_s_s_s[13 * 6];
379                     for(n = 0; n < 13; ++n)  // loop over orders of auxiliary function
380                     {
381 
382                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
383                         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]);
384                         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]);
385 
386                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
387                         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]);
388 
389                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
390                         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]);
391 
392                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
393                         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]);
394                         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]);
395 
396                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
397                         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]);
398 
399                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
400                         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]);
401                         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]);
402 
403                     }
404 
405 
406 
407                     // Forming PRIM_INT__f_s_s_s[12 * 10];
408                     for(n = 0; n < 12; ++n)  // loop over orders of auxiliary function
409                     {
410 
411                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
412                         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]);
413                         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]);
414 
415                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
416                         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]);
417 
418                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
419                         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]);
420 
421                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
422                         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]);
423 
424                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
425                         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]);
426 
427                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
428                         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]);
429 
430                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
431                         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]);
432                         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]);
433 
434                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
435                         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]);
436 
437                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
438                         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]);
439 
440                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
441                         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]);
442                         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]);
443 
444                     }
445 
446 
447                     VRR_I_g_s_s_s(
448                             PRIM_INT__g_s_s_s,
449                             PRIM_INT__f_s_s_s,
450                             PRIM_INT__d_s_s_s,
451                             P_PA,
452                             a_over_p,
453                             aop_PQ,
454                             one_over_2p,
455                             11);
456 
457 
458                     VRR_I_h_s_s_s(
459                             PRIM_INT__h_s_s_s,
460                             PRIM_INT__g_s_s_s,
461                             PRIM_INT__f_s_s_s,
462                             P_PA,
463                             a_over_p,
464                             aop_PQ,
465                             one_over_2p,
466                             10);
467 
468 
469                     ostei_general_vrr_K(5, 0, 1, 0, 5,
470                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
471                             PRIM_INT__h_s_s_s, NULL, NULL, PRIM_INT__g_s_s_s, NULL, PRIM_INT__h_s_p_s);
472 
473 
474                     VRR_K_g_s_p_s(
475                             PRIM_INT__g_s_p_s,
476                             PRIM_INT__g_s_s_s,
477                             PRIM_INT__f_s_s_s,
478                             Q_PA,
479                             aoq_PQ,
480                             one_over_2pq,
481                             5);
482 
483 
484                     ostei_general_vrr_K(5, 0, 2, 0, 4,
485                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
486                             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);
487 
488 
489                     VRR_K_f_s_p_s(
490                             PRIM_INT__f_s_p_s,
491                             PRIM_INT__f_s_s_s,
492                             PRIM_INT__d_s_s_s,
493                             Q_PA,
494                             aoq_PQ,
495                             one_over_2pq,
496                             5);
497 
498 
499                     ostei_general_vrr_K(4, 0, 2, 0, 4,
500                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
501                             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);
502 
503 
504                     ostei_general_vrr_K(5, 0, 3, 0, 3,
505                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
506                             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);
507 
508 
509 
510                     // Forming PRIM_INT__d_s_p_s[5 * 18];
511                     for(n = 0; n < 5; ++n)  // loop over orders of auxiliary function
512                     {
513 
514                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
515                         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]);
516                         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]);
517 
518                         PRIM_INT__d_s_p_s[n * 18 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
519                         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]);
520 
521                         PRIM_INT__d_s_p_s[n * 18 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
522                         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]);
523 
524                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 1]);
525                         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]);
526                         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]);
527 
528                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 1]);
529                         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]);
530                         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]);
531 
532                         PRIM_INT__d_s_p_s[n * 18 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
533                         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]);
534 
535                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 2]);
536                         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]);
537                         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]);
538 
539                         PRIM_INT__d_s_p_s[n * 18 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 2]);
540                         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]);
541 
542                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 2]);
543                         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]);
544                         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]);
545 
546                         PRIM_INT__d_s_p_s[n * 18 + 9] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
547                         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]);
548 
549                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
550                         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]);
551                         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]);
552 
553                         PRIM_INT__d_s_p_s[n * 18 + 11] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
554                         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]);
555 
556                         PRIM_INT__d_s_p_s[n * 18 + 12] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 4]);
557                         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]);
558 
559                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 4]);
560                         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]);
561                         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]);
562 
563                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 4]);
564                         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]);
565                         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]);
566 
567                         PRIM_INT__d_s_p_s[n * 18 + 15] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
568                         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]);
569 
570                         PRIM_INT__d_s_p_s[n * 18 + 16] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
571                         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]);
572 
573                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
574                         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]);
575                         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]);
576 
577                     }
578 
579 
580                     VRR_K_f_s_d_s(
581                             PRIM_INT__f_s_d_s,
582                             PRIM_INT__f_s_p_s,
583                             PRIM_INT__f_s_s_s,
584                             PRIM_INT__d_s_p_s,
585                             Q_PA,
586                             a_over_q,
587                             aoq_PQ,
588                             one_over_2pq,
589                             one_over_2q,
590                             4);
591 
592 
593                     ostei_general_vrr_K(4, 0, 3, 0, 3,
594                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
595                             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);
596 
597 
598                     ostei_general_vrr_K(5, 0, 4, 0, 2,
599                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
600                             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);
601 
602 
603                     ostei_general_vrr1_I(6, 9,
604                             one_over_2p, a_over_p, aop_PQ, P_PA,
605                             PRIM_INT__h_s_s_s, PRIM_INT__g_s_s_s, PRIM_INT__i_s_s_s);
606 
607 
608                     ostei_general_vrr_K(6, 0, 1, 0, 5,
609                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
610                             PRIM_INT__i_s_s_s, NULL, NULL, PRIM_INT__h_s_s_s, NULL, PRIM_INT__i_s_p_s);
611 
612 
613                     ostei_general_vrr_K(6, 0, 2, 0, 4,
614                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
615                             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);
616 
617 
618                     ostei_general_vrr_K(6, 0, 3, 0, 3,
619                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
620                             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);
621 
622 
623 
624                     // Forming PRIM_INT__p_s_p_s[5 * 9];
625                     for(n = 0; n < 5; ++n)  // loop over orders of auxiliary function
626                     {
627 
628                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
629                         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]);
630                         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]);
631 
632                         PRIM_INT__p_s_p_s[n * 9 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
633                         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]);
634 
635                         PRIM_INT__p_s_p_s[n * 9 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
636                         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]);
637 
638                         PRIM_INT__p_s_p_s[n * 9 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 1]);
639                         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]);
640 
641                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
642                         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]);
643                         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]);
644 
645                         PRIM_INT__p_s_p_s[n * 9 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
646                         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]);
647 
648                         PRIM_INT__p_s_p_s[n * 9 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 2]);
649                         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]);
650 
651                         PRIM_INT__p_s_p_s[n * 9 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 2]);
652                         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]);
653 
654                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
655                         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]);
656                         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]);
657 
658                     }
659 
660 
661                     VRR_K_d_s_d_s(
662                             PRIM_INT__d_s_d_s,
663                             PRIM_INT__d_s_p_s,
664                             PRIM_INT__d_s_s_s,
665                             PRIM_INT__p_s_p_s,
666                             Q_PA,
667                             a_over_q,
668                             aoq_PQ,
669                             one_over_2pq,
670                             one_over_2q,
671                             4);
672 
673 
674                     ostei_general_vrr_K(3, 0, 3, 0, 3,
675                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
676                             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);
677 
678 
679                     ostei_general_vrr_K(4, 0, 4, 0, 2,
680                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
681                             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);
682 
683 
684                     ostei_general_vrr_K(5, 0, 5, 0, 1,
685                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
686                             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);
687 
688 
689                     ostei_general_vrr_K(6, 0, 4, 0, 2,
690                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
691                             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);
692 
693 
694                     ostei_general_vrr1_I(7, 8,
695                             one_over_2p, a_over_p, aop_PQ, P_PA,
696                             PRIM_INT__i_s_s_s, PRIM_INT__h_s_s_s, PRIM_INT__k_s_s_s);
697 
698 
699                     ostei_general_vrr_K(7, 0, 1, 0, 5,
700                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
701                             PRIM_INT__k_s_s_s, NULL, NULL, PRIM_INT__i_s_s_s, NULL, PRIM_INT__k_s_p_s);
702 
703 
704                     ostei_general_vrr_K(7, 0, 2, 0, 4,
705                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
706                             PRIM_INT__k_s_p_s, PRIM_INT__k_s_s_s, NULL, PRIM_INT__i_s_p_s, NULL, PRIM_INT__k_s_d_s);
707 
708 
709                     ostei_general_vrr_K(7, 0, 3, 0, 3,
710                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
711                             PRIM_INT__k_s_d_s, PRIM_INT__k_s_p_s, NULL, PRIM_INT__i_s_d_s, NULL, PRIM_INT__k_s_f_s);
712 
713 
714                     ostei_general_vrr_K(6, 0, 5, 0, 1,
715                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
716                             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);
717 
718 
719                     ostei_general_vrr_K(7, 0, 4, 0, 2,
720                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
721                             PRIM_INT__k_s_f_s, PRIM_INT__k_s_d_s, NULL, PRIM_INT__i_s_f_s, NULL, PRIM_INT__k_s_g_s);
722 
723 
724                     ostei_general_vrr1_I(8, 7,
725                             one_over_2p, a_over_p, aop_PQ, P_PA,
726                             PRIM_INT__k_s_s_s, PRIM_INT__i_s_s_s, PRIM_INT__l_s_s_s);
727 
728 
729                     ostei_general_vrr_K(8, 0, 1, 0, 5,
730                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
731                             PRIM_INT__l_s_s_s, NULL, NULL, PRIM_INT__k_s_s_s, NULL, PRIM_INT__l_s_p_s);
732 
733 
734                     ostei_general_vrr_K(8, 0, 2, 0, 4,
735                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
736                             PRIM_INT__l_s_p_s, PRIM_INT__l_s_s_s, NULL, PRIM_INT__k_s_p_s, NULL, PRIM_INT__l_s_d_s);
737 
738 
739                     ostei_general_vrr_K(8, 0, 3, 0, 3,
740                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
741                             PRIM_INT__l_s_d_s, PRIM_INT__l_s_p_s, NULL, PRIM_INT__k_s_d_s, NULL, PRIM_INT__l_s_f_s);
742 
743 
744                     ostei_general_vrr_K(7, 0, 5, 0, 1,
745                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
746                             PRIM_INT__k_s_g_s, PRIM_INT__k_s_f_s, NULL, PRIM_INT__i_s_g_s, NULL, PRIM_INT__k_s_h_s);
747 
748 
749                     ostei_general_vrr_K(8, 0, 4, 0, 2,
750                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
751                             PRIM_INT__l_s_f_s, PRIM_INT__l_s_d_s, NULL, PRIM_INT__k_s_f_s, NULL, PRIM_INT__l_s_g_s);
752 
753 
754                     ostei_general_vrr1_I(9, 6,
755                             one_over_2p, a_over_p, aop_PQ, P_PA,
756                             PRIM_INT__l_s_s_s, PRIM_INT__k_s_s_s, PRIM_INT__m_s_s_s);
757 
758 
759                     ostei_general_vrr_K(9, 0, 1, 0, 5,
760                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
761                             PRIM_INT__m_s_s_s, NULL, NULL, PRIM_INT__l_s_s_s, NULL, PRIM_INT__m_s_p_s);
762 
763 
764                     ostei_general_vrr_K(9, 0, 2, 0, 4,
765                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
766                             PRIM_INT__m_s_p_s, PRIM_INT__m_s_s_s, NULL, PRIM_INT__l_s_p_s, NULL, PRIM_INT__m_s_d_s);
767 
768 
769                     ostei_general_vrr_K(9, 0, 3, 0, 3,
770                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
771                             PRIM_INT__m_s_d_s, PRIM_INT__m_s_p_s, NULL, PRIM_INT__l_s_d_s, NULL, PRIM_INT__m_s_f_s);
772 
773 
774                     ostei_general_vrr_K(8, 0, 5, 0, 1,
775                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
776                             PRIM_INT__l_s_g_s, PRIM_INT__l_s_f_s, NULL, PRIM_INT__k_s_g_s, NULL, PRIM_INT__l_s_h_s);
777 
778 
779                     ostei_general_vrr_K(9, 0, 4, 0, 2,
780                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
781                             PRIM_INT__m_s_f_s, PRIM_INT__m_s_d_s, NULL, PRIM_INT__l_s_f_s, NULL, PRIM_INT__m_s_g_s);
782 
783 
784                     ostei_general_vrr_K(9, 0, 5, 0, 1,
785                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
786                             PRIM_INT__m_s_g_s, PRIM_INT__m_s_f_s, NULL, PRIM_INT__l_s_g_s, NULL, PRIM_INT__m_s_h_s);
787 
788 
789 
790 
791                     ////////////////////////////////////
792                     // Accumulate contracted integrals
793                     ////////////////////////////////////
794                     if(lastoffset == 0)
795                     {
796                         contract_all(210, PRIM_INT__h_s_f_s, PRIM_PTR_INT__h_s_f_s);
797                         contract_all(315, PRIM_INT__h_s_g_s, PRIM_PTR_INT__h_s_g_s);
798                         contract_all(441, PRIM_INT__h_s_h_s, PRIM_PTR_INT__h_s_h_s);
799                         contract_all(280, PRIM_INT__i_s_f_s, PRIM_PTR_INT__i_s_f_s);
800                         contract_all(420, PRIM_INT__i_s_g_s, PRIM_PTR_INT__i_s_g_s);
801                         contract_all(588, PRIM_INT__i_s_h_s, PRIM_PTR_INT__i_s_h_s);
802                         contract_all(360, PRIM_INT__k_s_f_s, PRIM_PTR_INT__k_s_f_s);
803                         contract_all(540, PRIM_INT__k_s_g_s, PRIM_PTR_INT__k_s_g_s);
804                         contract_all(756, PRIM_INT__k_s_h_s, PRIM_PTR_INT__k_s_h_s);
805                         contract_all(450, PRIM_INT__l_s_f_s, PRIM_PTR_INT__l_s_f_s);
806                         contract_all(675, PRIM_INT__l_s_g_s, PRIM_PTR_INT__l_s_g_s);
807                         contract_all(945, PRIM_INT__l_s_h_s, PRIM_PTR_INT__l_s_h_s);
808                         contract_all(550, PRIM_INT__m_s_f_s, PRIM_PTR_INT__m_s_f_s);
809                         contract_all(825, PRIM_INT__m_s_g_s, PRIM_PTR_INT__m_s_g_s);
810                         contract_all(1155, PRIM_INT__m_s_h_s, PRIM_PTR_INT__m_s_h_s);
811                     }
812                     else
813                     {
814                         contract(210, shelloffsets, PRIM_INT__h_s_f_s, PRIM_PTR_INT__h_s_f_s);
815                         contract(315, shelloffsets, PRIM_INT__h_s_g_s, PRIM_PTR_INT__h_s_g_s);
816                         contract(441, shelloffsets, PRIM_INT__h_s_h_s, PRIM_PTR_INT__h_s_h_s);
817                         contract(280, shelloffsets, PRIM_INT__i_s_f_s, PRIM_PTR_INT__i_s_f_s);
818                         contract(420, shelloffsets, PRIM_INT__i_s_g_s, PRIM_PTR_INT__i_s_g_s);
819                         contract(588, shelloffsets, PRIM_INT__i_s_h_s, PRIM_PTR_INT__i_s_h_s);
820                         contract(360, shelloffsets, PRIM_INT__k_s_f_s, PRIM_PTR_INT__k_s_f_s);
821                         contract(540, shelloffsets, PRIM_INT__k_s_g_s, PRIM_PTR_INT__k_s_g_s);
822                         contract(756, shelloffsets, PRIM_INT__k_s_h_s, PRIM_PTR_INT__k_s_h_s);
823                         contract(450, shelloffsets, PRIM_INT__l_s_f_s, PRIM_PTR_INT__l_s_f_s);
824                         contract(675, shelloffsets, PRIM_INT__l_s_g_s, PRIM_PTR_INT__l_s_g_s);
825                         contract(945, shelloffsets, PRIM_INT__l_s_h_s, PRIM_PTR_INT__l_s_h_s);
826                         contract(550, shelloffsets, PRIM_INT__m_s_f_s, PRIM_PTR_INT__m_s_f_s);
827                         contract(825, shelloffsets, PRIM_INT__m_s_g_s, PRIM_PTR_INT__m_s_g_s);
828                         contract(1155, shelloffsets, PRIM_INT__m_s_h_s, PRIM_PTR_INT__m_s_h_s);
829                         PRIM_PTR_INT__h_s_f_s += lastoffset*210;
830                         PRIM_PTR_INT__h_s_g_s += lastoffset*315;
831                         PRIM_PTR_INT__h_s_h_s += lastoffset*441;
832                         PRIM_PTR_INT__i_s_f_s += lastoffset*280;
833                         PRIM_PTR_INT__i_s_g_s += lastoffset*420;
834                         PRIM_PTR_INT__i_s_h_s += lastoffset*588;
835                         PRIM_PTR_INT__k_s_f_s += lastoffset*360;
836                         PRIM_PTR_INT__k_s_g_s += lastoffset*540;
837                         PRIM_PTR_INT__k_s_h_s += lastoffset*756;
838                         PRIM_PTR_INT__l_s_f_s += lastoffset*450;
839                         PRIM_PTR_INT__l_s_g_s += lastoffset*675;
840                         PRIM_PTR_INT__l_s_h_s += lastoffset*945;
841                         PRIM_PTR_INT__m_s_f_s += lastoffset*550;
842                         PRIM_PTR_INT__m_s_g_s += lastoffset*825;
843                         PRIM_PTR_INT__m_s_h_s += lastoffset*1155;
844                     }
845 
846                 }  // close loop over j
847             }  // close loop over i
848 
849             //Advance to the next batch
850             jstart = SIMINT_SIMD_ROUND(jend);
851 
852             //////////////////////////////////////////////
853             // Contracted integrals: Horizontal recurrance
854             //////////////////////////////////////////////
855 
856 
857             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
858 
859 
860             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
861             {
862                 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
863 
864                 // set up HRR pointers
865                 double const * restrict HRR_INT__h_s_f_s = INT__h_s_f_s + abcd * 210;
866                 double const * restrict HRR_INT__h_s_g_s = INT__h_s_g_s + abcd * 315;
867                 double const * restrict HRR_INT__h_s_h_s = INT__h_s_h_s + abcd * 441;
868                 double const * restrict HRR_INT__i_s_f_s = INT__i_s_f_s + abcd * 280;
869                 double const * restrict HRR_INT__i_s_g_s = INT__i_s_g_s + abcd * 420;
870                 double const * restrict HRR_INT__i_s_h_s = INT__i_s_h_s + abcd * 588;
871                 double const * restrict HRR_INT__k_s_f_s = INT__k_s_f_s + abcd * 360;
872                 double const * restrict HRR_INT__k_s_g_s = INT__k_s_g_s + abcd * 540;
873                 double const * restrict HRR_INT__k_s_h_s = INT__k_s_h_s + abcd * 756;
874                 double const * restrict HRR_INT__l_s_f_s = INT__l_s_f_s + abcd * 450;
875                 double const * restrict HRR_INT__l_s_g_s = INT__l_s_g_s + abcd * 675;
876                 double const * restrict HRR_INT__l_s_h_s = INT__l_s_h_s + abcd * 945;
877                 double const * restrict HRR_INT__m_s_f_s = INT__m_s_f_s + abcd * 550;
878                 double const * restrict HRR_INT__m_s_g_s = INT__m_s_g_s + abcd * 825;
879                 double const * restrict HRR_INT__m_s_h_s = INT__m_s_h_s + abcd * 1155;
880                 double * restrict HRR_INT__h_g_f_d = INT__h_g_f_d + real_abcd * 18900;
881 
882                 // form INT__h_p_f_s
883                 ostei_general_hrr_J(5, 1, 3, 0, hAB, HRR_INT__i_s_f_s, HRR_INT__h_s_f_s, HRR_INT__h_p_f_s);
884 
885                 // form INT__h_p_g_s
886                 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);
887 
888                 // form INT__h_p_h_s
889                 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);
890 
891                 // form INT__i_p_f_s
892                 ostei_general_hrr_J(6, 1, 3, 0, hAB, HRR_INT__k_s_f_s, HRR_INT__i_s_f_s, HRR_INT__i_p_f_s);
893 
894                 // form INT__i_p_g_s
895                 ostei_general_hrr_J(6, 1, 4, 0, hAB, HRR_INT__k_s_g_s, HRR_INT__i_s_g_s, HRR_INT__i_p_g_s);
896 
897                 // form INT__i_p_h_s
898                 ostei_general_hrr_J(6, 1, 5, 0, hAB, HRR_INT__k_s_h_s, HRR_INT__i_s_h_s, HRR_INT__i_p_h_s);
899 
900                 // form INT__k_p_f_s
901                 ostei_general_hrr_J(7, 1, 3, 0, hAB, HRR_INT__l_s_f_s, HRR_INT__k_s_f_s, HRR_INT__k_p_f_s);
902 
903                 // form INT__k_p_g_s
904                 ostei_general_hrr_J(7, 1, 4, 0, hAB, HRR_INT__l_s_g_s, HRR_INT__k_s_g_s, HRR_INT__k_p_g_s);
905 
906                 // form INT__k_p_h_s
907                 ostei_general_hrr_J(7, 1, 5, 0, hAB, HRR_INT__l_s_h_s, HRR_INT__k_s_h_s, HRR_INT__k_p_h_s);
908 
909                 // form INT__l_p_f_s
910                 ostei_general_hrr_J(8, 1, 3, 0, hAB, HRR_INT__m_s_f_s, HRR_INT__l_s_f_s, HRR_INT__l_p_f_s);
911 
912                 // form INT__l_p_g_s
913                 ostei_general_hrr_J(8, 1, 4, 0, hAB, HRR_INT__m_s_g_s, HRR_INT__l_s_g_s, HRR_INT__l_p_g_s);
914 
915                 // form INT__l_p_h_s
916                 ostei_general_hrr_J(8, 1, 5, 0, hAB, HRR_INT__m_s_h_s, HRR_INT__l_s_h_s, HRR_INT__l_p_h_s);
917 
918                 // form INT__h_d_f_s
919                 ostei_general_hrr_J(5, 2, 3, 0, hAB, HRR_INT__i_p_f_s, HRR_INT__h_p_f_s, HRR_INT__h_d_f_s);
920 
921                 // form INT__h_d_g_s
922                 ostei_general_hrr_J(5, 2, 4, 0, hAB, HRR_INT__i_p_g_s, HRR_INT__h_p_g_s, HRR_INT__h_d_g_s);
923 
924                 // form INT__h_d_h_s
925                 ostei_general_hrr_J(5, 2, 5, 0, hAB, HRR_INT__i_p_h_s, HRR_INT__h_p_h_s, HRR_INT__h_d_h_s);
926 
927                 // form INT__i_d_f_s
928                 ostei_general_hrr_J(6, 2, 3, 0, hAB, HRR_INT__k_p_f_s, HRR_INT__i_p_f_s, HRR_INT__i_d_f_s);
929 
930                 // form INT__i_d_g_s
931                 ostei_general_hrr_J(6, 2, 4, 0, hAB, HRR_INT__k_p_g_s, HRR_INT__i_p_g_s, HRR_INT__i_d_g_s);
932 
933                 // form INT__i_d_h_s
934                 ostei_general_hrr_J(6, 2, 5, 0, hAB, HRR_INT__k_p_h_s, HRR_INT__i_p_h_s, HRR_INT__i_d_h_s);
935 
936                 // form INT__k_d_f_s
937                 ostei_general_hrr_J(7, 2, 3, 0, hAB, HRR_INT__l_p_f_s, HRR_INT__k_p_f_s, HRR_INT__k_d_f_s);
938 
939                 // form INT__k_d_g_s
940                 ostei_general_hrr_J(7, 2, 4, 0, hAB, HRR_INT__l_p_g_s, HRR_INT__k_p_g_s, HRR_INT__k_d_g_s);
941 
942                 // form INT__k_d_h_s
943                 ostei_general_hrr_J(7, 2, 5, 0, hAB, HRR_INT__l_p_h_s, HRR_INT__k_p_h_s, HRR_INT__k_d_h_s);
944 
945                 // form INT__h_f_f_s
946                 ostei_general_hrr_J(5, 3, 3, 0, hAB, HRR_INT__i_d_f_s, HRR_INT__h_d_f_s, HRR_INT__h_f_f_s);
947 
948                 // form INT__h_f_g_s
949                 ostei_general_hrr_J(5, 3, 4, 0, hAB, HRR_INT__i_d_g_s, HRR_INT__h_d_g_s, HRR_INT__h_f_g_s);
950 
951                 // form INT__h_f_h_s
952                 ostei_general_hrr_J(5, 3, 5, 0, hAB, HRR_INT__i_d_h_s, HRR_INT__h_d_h_s, HRR_INT__h_f_h_s);
953 
954                 // form INT__i_f_f_s
955                 ostei_general_hrr_J(6, 3, 3, 0, hAB, HRR_INT__k_d_f_s, HRR_INT__i_d_f_s, HRR_INT__i_f_f_s);
956 
957                 // form INT__i_f_g_s
958                 ostei_general_hrr_J(6, 3, 4, 0, hAB, HRR_INT__k_d_g_s, HRR_INT__i_d_g_s, HRR_INT__i_f_g_s);
959 
960                 // form INT__i_f_h_s
961                 ostei_general_hrr_J(6, 3, 5, 0, hAB, HRR_INT__k_d_h_s, HRR_INT__i_d_h_s, HRR_INT__i_f_h_s);
962 
963                 // form INT__h_g_f_s
964                 ostei_general_hrr_J(5, 4, 3, 0, hAB, HRR_INT__i_f_f_s, HRR_INT__h_f_f_s, HRR_INT__h_g_f_s);
965 
966                 // form INT__h_g_g_s
967                 ostei_general_hrr_J(5, 4, 4, 0, hAB, HRR_INT__i_f_g_s, HRR_INT__h_f_g_s, HRR_INT__h_g_g_s);
968 
969                 // form INT__h_g_h_s
970                 ostei_general_hrr_J(5, 4, 5, 0, hAB, HRR_INT__i_f_h_s, HRR_INT__h_f_h_s, HRR_INT__h_g_h_s);
971 
972                 // form INT__h_g_f_p
973                 HRR_L_f_p(
974                     HRR_INT__h_g_f_p,
975                     HRR_INT__h_g_f_s,
976                     HRR_INT__h_g_g_s,
977                     hCD, 315);
978 
979                 // form INT__h_g_g_p
980                 HRR_L_g_p(
981                     HRR_INT__h_g_g_p,
982                     HRR_INT__h_g_g_s,
983                     HRR_INT__h_g_h_s,
984                     hCD, 315);
985 
986                 // form INT__h_g_f_d
987                 HRR_L_f_d(
988                     HRR_INT__h_g_f_d,
989                     HRR_INT__h_g_f_p,
990                     HRR_INT__h_g_g_p,
991                     hCD, 315);
992 
993 
994             }  // close HRR loop
995 
996 
997         }   // close loop cdbatch
998 
999         istart = iend;
1000     }  // close loop over ab
1001 
1002     return P.nshell12_clip * Q.nshell12_clip;
1003 }
1004 
ostei_g_h_f_d(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__g_h_f_d)1005 int ostei_g_h_f_d(struct simint_multi_shellpair const P,
1006                   struct simint_multi_shellpair const Q,
1007                   double screen_tol,
1008                   double * const restrict work,
1009                   double * const restrict INT__g_h_f_d)
1010 {
1011     double P_AB[3*P.nshell12];
1012     struct simint_multi_shellpair P_tmp = P;
1013     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
1014     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
1015     P_tmp.AB_x = P_AB;
1016     P_tmp.AB_y = P_AB + P.nshell12;
1017     P_tmp.AB_z = P_AB + 2*P.nshell12;
1018 
1019     for(int i = 0; i < P.nshell12; i++)
1020     {
1021         P_tmp.AB_x[i] = -P.AB_x[i];
1022         P_tmp.AB_y[i] = -P.AB_y[i];
1023         P_tmp.AB_z[i] = -P.AB_z[i];
1024     }
1025 
1026     int ret = ostei_h_g_f_d(P_tmp, Q, screen_tol, work, INT__g_h_f_d);
1027     double buffer[18900] SIMINT_ALIGN_ARRAY_DBL;
1028 
1029     for(int q = 0; q < ret; q++)
1030     {
1031         int idx = 0;
1032         for(int a = 0; a < 15; ++a)
1033         for(int b = 0; b < 21; ++b)
1034         for(int c = 0; c < 10; ++c)
1035         for(int d = 0; d < 6; ++d)
1036             buffer[idx++] = INT__g_h_f_d[q*18900+b*900+a*60+c*6+d];
1037 
1038         memcpy(INT__g_h_f_d+q*18900, buffer, 18900*sizeof(double));
1039     }
1040 
1041     return ret;
1042 }
1043 
ostei_h_g_d_f(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__h_g_d_f)1044 int ostei_h_g_d_f(struct simint_multi_shellpair const P,
1045                   struct simint_multi_shellpair const Q,
1046                   double screen_tol,
1047                   double * const restrict work,
1048                   double * const restrict INT__h_g_d_f)
1049 {
1050     double Q_AB[3*Q.nshell12];
1051     struct simint_multi_shellpair Q_tmp = Q;
1052     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
1053     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
1054     Q_tmp.AB_x = Q_AB;
1055     Q_tmp.AB_y = Q_AB + Q.nshell12;
1056     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
1057 
1058     for(int i = 0; i < Q.nshell12; i++)
1059     {
1060         Q_tmp.AB_x[i] = -Q.AB_x[i];
1061         Q_tmp.AB_y[i] = -Q.AB_y[i];
1062         Q_tmp.AB_z[i] = -Q.AB_z[i];
1063     }
1064 
1065     int ret = ostei_h_g_f_d(P, Q_tmp, screen_tol, work, INT__h_g_d_f);
1066     double buffer[18900] SIMINT_ALIGN_ARRAY_DBL;
1067 
1068     for(int q = 0; q < ret; q++)
1069     {
1070         int idx = 0;
1071         for(int a = 0; a < 21; ++a)
1072         for(int b = 0; b < 15; ++b)
1073         for(int c = 0; c < 6; ++c)
1074         for(int d = 0; d < 10; ++d)
1075             buffer[idx++] = INT__h_g_d_f[q*18900+a*900+b*60+d*6+c];
1076 
1077         memcpy(INT__h_g_d_f+q*18900, buffer, 18900*sizeof(double));
1078     }
1079 
1080     return ret;
1081 }
1082 
ostei_g_h_d_f(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__g_h_d_f)1083 int ostei_g_h_d_f(struct simint_multi_shellpair const P,
1084                   struct simint_multi_shellpair const Q,
1085                   double screen_tol,
1086                   double * const restrict work,
1087                   double * const restrict INT__g_h_d_f)
1088 {
1089     double P_AB[3*P.nshell12];
1090     struct simint_multi_shellpair P_tmp = P;
1091     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
1092     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
1093     P_tmp.AB_x = P_AB;
1094     P_tmp.AB_y = P_AB + P.nshell12;
1095     P_tmp.AB_z = P_AB + 2*P.nshell12;
1096 
1097     for(int i = 0; i < P.nshell12; i++)
1098     {
1099         P_tmp.AB_x[i] = -P.AB_x[i];
1100         P_tmp.AB_y[i] = -P.AB_y[i];
1101         P_tmp.AB_z[i] = -P.AB_z[i];
1102     }
1103 
1104     double Q_AB[3*Q.nshell12];
1105     struct simint_multi_shellpair Q_tmp = Q;
1106     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
1107     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
1108     Q_tmp.AB_x = Q_AB;
1109     Q_tmp.AB_y = Q_AB + Q.nshell12;
1110     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
1111 
1112     for(int i = 0; i < Q.nshell12; i++)
1113     {
1114         Q_tmp.AB_x[i] = -Q.AB_x[i];
1115         Q_tmp.AB_y[i] = -Q.AB_y[i];
1116         Q_tmp.AB_z[i] = -Q.AB_z[i];
1117     }
1118 
1119     int ret = ostei_h_g_f_d(P_tmp, Q_tmp, screen_tol, work, INT__g_h_d_f);
1120     double buffer[18900] SIMINT_ALIGN_ARRAY_DBL;
1121 
1122     for(int q = 0; q < ret; q++)
1123     {
1124         int idx = 0;
1125         for(int a = 0; a < 15; ++a)
1126         for(int b = 0; b < 21; ++b)
1127         for(int c = 0; c < 6; ++c)
1128         for(int d = 0; d < 10; ++d)
1129             buffer[idx++] = INT__g_h_d_f[q*18900+b*900+a*60+d*6+c];
1130 
1131         memcpy(INT__g_h_d_f+q*18900, buffer, 18900*sizeof(double));
1132     }
1133 
1134     return ret;
1135 }
1136 
1137