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_i_d_g_p(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__i_d_g_p)8 int ostei_i_d_g_p(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__i_d_g_p)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__i_d_g_p);
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__i_s_g_s = work + (SIMINT_NSHELL_SIMD * 0);
30     double * const INT__i_s_h_s = work + (SIMINT_NSHELL_SIMD * 420);
31     double * const INT__k_s_g_s = work + (SIMINT_NSHELL_SIMD * 1008);
32     double * const INT__k_s_h_s = work + (SIMINT_NSHELL_SIMD * 1548);
33     double * const INT__l_s_g_s = work + (SIMINT_NSHELL_SIMD * 2304);
34     double * const INT__l_s_h_s = work + (SIMINT_NSHELL_SIMD * 2979);
35     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*3924);
36     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
37     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_s = primwork + 14;
38     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 53;
39     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_p_s = primwork + 125;
40     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_s = primwork + 215;
41     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_p_s = primwork + 325;
42     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_d_s = primwork + 475;
43     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_s_s = primwork + 715;
44     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_p_s = primwork + 865;
45     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_d_s = primwork + 1090;
46     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_f_s = primwork + 1450;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_s_s = primwork + 1900;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_p_s = primwork + 2089;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_d_s = primwork + 2404;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_f_s = primwork + 2908;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_g_s = primwork + 3538;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_s_s = primwork + 4168;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_p_s = primwork + 4392;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_d_s = primwork + 4812;
55     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_f_s = primwork + 5484;
56     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_g_s = primwork + 6324;
57     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_h_s = primwork + 7164;
58     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_s_s = primwork + 7752;
59     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_p_s = primwork + 8004;
60     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_d_s = primwork + 8544;
61     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_f_s = primwork + 9408;
62     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_g_s = primwork + 10488;
63     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_h_s = primwork + 11568;
64     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_s_s = primwork + 12324;
65     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_p_s = primwork + 12594;
66     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_d_s = primwork + 13269;
67     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_f_s = primwork + 14349;
68     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_g_s = primwork + 15699;
69     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_h_s = primwork + 17049;
70     double * const hrrwork = (double *)(primwork + 17994);
71     double * const HRR_INT__i_p_g_s = hrrwork + 0;
72     double * const HRR_INT__i_p_h_s = hrrwork + 1260;
73     double * const HRR_INT__i_d_g_s = hrrwork + 3024;
74     double * const HRR_INT__i_d_h_s = hrrwork + 5544;
75     double * const HRR_INT__k_p_g_s = hrrwork + 9072;
76     double * const HRR_INT__k_p_h_s = hrrwork + 10692;
77 
78 
79     // Create constants
80     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
81     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
82     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
83     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
84     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
85     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
86     const SIMINT_DBLTYPE const_7 = SIMINT_DBLSET1(7);
87     const SIMINT_DBLTYPE const_8 = SIMINT_DBLSET1(8);
88     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
89 
90 
91     ////////////////////////////////////////
92     // Loop over shells and primitives
93     ////////////////////////////////////////
94 
95     real_abcd = 0;
96     istart = 0;
97     for(ab = 0; ab < P.nshell12_clip; ++ab)
98     {
99         const int iend = istart + P.nprim12[ab];
100 
101         cd = 0;
102         jstart = 0;
103 
104         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
105         {
106             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
107             int jend = jstart;
108             for(i = 0; i < nshellbatch; i++)
109                 jend += Q.nprim12[cd+i];
110 
111             // Clear the beginning of the workspace (where we are accumulating integrals)
112             memset(work, 0, SIMINT_NSHELL_SIMD * 3924 * sizeof(double));
113             abcd = 0;
114 
115 
116             for(i = istart; i < iend; ++i)
117             {
118                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
119 
120                 if(check_screen)
121                 {
122                     // Skip this whole thing if always insignificant
123                     if((P.screen[i] * Q.screen_max) < screen_tol)
124                         continue;
125                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
126                 }
127 
128                 icd = 0;
129                 iprimcd = 0;
130                 nprim_icd = Q.nprim12[cd];
131                 double * restrict PRIM_PTR_INT__i_s_g_s = INT__i_s_g_s + abcd * 420;
132                 double * restrict PRIM_PTR_INT__i_s_h_s = INT__i_s_h_s + abcd * 588;
133                 double * restrict PRIM_PTR_INT__k_s_g_s = INT__k_s_g_s + abcd * 540;
134                 double * restrict PRIM_PTR_INT__k_s_h_s = INT__k_s_h_s + abcd * 756;
135                 double * restrict PRIM_PTR_INT__l_s_g_s = INT__l_s_g_s + abcd * 675;
136                 double * restrict PRIM_PTR_INT__l_s_h_s = INT__l_s_h_s + abcd * 945;
137 
138 
139 
140                 // Load these one per loop over i
141                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
142                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
143                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
144 
145                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
146 
147                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
148                 {
149                     // calculate the shell offsets
150                     // these are the offset from the shell pointed to by cd
151                     // for each element
152                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
153                     int lastoffset = 0;
154                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
155 
156                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
157                     {
158                         // Handle if the first element of the vector is a new shell
159                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
160                         {
161                             nprim_icd += Q.nprim12[cd + (++icd)];
162                             PRIM_PTR_INT__i_s_g_s += 420;
163                             PRIM_PTR_INT__i_s_h_s += 588;
164                             PRIM_PTR_INT__k_s_g_s += 540;
165                             PRIM_PTR_INT__k_s_h_s += 756;
166                             PRIM_PTR_INT__l_s_g_s += 675;
167                             PRIM_PTR_INT__l_s_h_s += 945;
168                         }
169                         iprimcd++;
170                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
171                         {
172                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
173                             {
174                                 shelloffsets[n] = shelloffsets[n-1] + 1;
175                                 lastoffset++;
176                                 nprim_icd += Q.nprim12[cd + (++icd)];
177                             }
178                             else
179                                 shelloffsets[n] = shelloffsets[n-1];
180                             iprimcd++;
181                         }
182                     }
183                     else
184                         iprimcd += SIMINT_SIMD_LEN;
185 
186                     // Do we have to compute this vector (or has it been screened out)?
187                     // (not_screened != 0 means we have to do this vector)
188                     if(check_screen)
189                     {
190                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
191                         if(vmax < screen_tol)
192                         {
193                             PRIM_PTR_INT__i_s_g_s += lastoffset*420;
194                             PRIM_PTR_INT__i_s_h_s += lastoffset*588;
195                             PRIM_PTR_INT__k_s_g_s += lastoffset*540;
196                             PRIM_PTR_INT__k_s_h_s += lastoffset*756;
197                             PRIM_PTR_INT__l_s_g_s += lastoffset*675;
198                             PRIM_PTR_INT__l_s_h_s += lastoffset*945;
199                             continue;
200                         }
201                     }
202 
203                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
204                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
205                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
206                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
207 
208 
209                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
210                     SIMINT_DBLTYPE PQ[3];
211                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
212                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
213                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
214                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
215                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
216                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
217 
218                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
219                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
220                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
221                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
222                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
223                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
224                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
225 
226                     // NOTE: Minus sign!
227                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
228                     SIMINT_DBLTYPE aop_PQ[3];
229                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
230                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
231                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
232 
233                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
234                     SIMINT_DBLTYPE aoq_PQ[3];
235                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
236                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
237                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
238                     // Put a minus sign here so we don't have to in RR routines
239                     a_over_q = SIMINT_NEG(a_over_q);
240 
241 
242                     //////////////////////////////////////////////
243                     // Fjt function section
244                     // Maximum v value: 13
245                     //////////////////////////////////////////////
246                     // The parameter to the Fjt function
247                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
248 
249 
250                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
251 
252 
253                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 13);
254                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
255                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
256                     for(n = 0; n <= 13; n++)
257                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
258 
259                     //////////////////////////////////////////////
260                     // Primitive integrals: Vertical recurrance
261                     //////////////////////////////////////////////
262 
263                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
264                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
265                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
266                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
267                     const SIMINT_DBLTYPE vrr_const_5_over_2p = SIMINT_MUL(const_5, one_over_2p);
268                     const SIMINT_DBLTYPE vrr_const_6_over_2p = SIMINT_MUL(const_6, one_over_2p);
269                     const SIMINT_DBLTYPE vrr_const_7_over_2p = SIMINT_MUL(const_7, one_over_2p);
270                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
271                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
272                     const SIMINT_DBLTYPE vrr_const_3_over_2q = SIMINT_MUL(const_3, one_over_2q);
273                     const SIMINT_DBLTYPE vrr_const_4_over_2q = SIMINT_MUL(const_4, one_over_2q);
274                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
275                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
276                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
277                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
278                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
279                     const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
280                     const SIMINT_DBLTYPE vrr_const_7_over_2pq = SIMINT_MUL(const_7, one_over_2pq);
281                     const SIMINT_DBLTYPE vrr_const_8_over_2pq = SIMINT_MUL(const_8, one_over_2pq);
282 
283 
284 
285                     // Forming PRIM_INT__p_s_s_s[13 * 3];
286                     for(n = 0; n < 13; ++n)  // loop over orders of auxiliary function
287                     {
288 
289                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
290                         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]);
291 
292                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
293                         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]);
294 
295                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
296                         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]);
297 
298                     }
299 
300 
301 
302                     // Forming PRIM_INT__d_s_s_s[12 * 6];
303                     for(n = 0; n < 12; ++n)  // loop over orders of auxiliary function
304                     {
305 
306                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
307                         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]);
308                         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]);
309 
310                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
311                         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]);
312 
313                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
314                         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]);
315 
316                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
317                         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]);
318                         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]);
319 
320                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
321                         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]);
322 
323                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
324                         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]);
325                         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]);
326 
327                     }
328 
329 
330 
331                     // Forming PRIM_INT__f_s_s_s[11 * 10];
332                     for(n = 0; n < 11; ++n)  // loop over orders of auxiliary function
333                     {
334 
335                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
336                         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]);
337                         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]);
338 
339                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
340                         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]);
341 
342                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
343                         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]);
344 
345                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
346                         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]);
347 
348                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
349                         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]);
350 
351                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
352                         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]);
353 
354                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
355                         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]);
356                         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]);
357 
358                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
359                         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]);
360 
361                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
362                         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]);
363 
364                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
365                         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]);
366                         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]);
367 
368                     }
369 
370 
371                     VRR_I_g_s_s_s(
372                             PRIM_INT__g_s_s_s,
373                             PRIM_INT__f_s_s_s,
374                             PRIM_INT__d_s_s_s,
375                             P_PA,
376                             a_over_p,
377                             aop_PQ,
378                             one_over_2p,
379                             10);
380 
381 
382                     VRR_I_h_s_s_s(
383                             PRIM_INT__h_s_s_s,
384                             PRIM_INT__g_s_s_s,
385                             PRIM_INT__f_s_s_s,
386                             P_PA,
387                             a_over_p,
388                             aop_PQ,
389                             one_over_2p,
390                             9);
391 
392 
393                     ostei_general_vrr1_I(6, 8,
394                             one_over_2p, a_over_p, aop_PQ, P_PA,
395                             PRIM_INT__h_s_s_s, PRIM_INT__g_s_s_s, PRIM_INT__i_s_s_s);
396 
397 
398                     ostei_general_vrr_K(6, 0, 1, 0, 5,
399                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
400                             PRIM_INT__i_s_s_s, NULL, NULL, PRIM_INT__h_s_s_s, NULL, PRIM_INT__i_s_p_s);
401 
402 
403                     ostei_general_vrr_K(5, 0, 1, 0, 5,
404                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
405                             PRIM_INT__h_s_s_s, NULL, NULL, PRIM_INT__g_s_s_s, NULL, PRIM_INT__h_s_p_s);
406 
407 
408                     ostei_general_vrr_K(6, 0, 2, 0, 4,
409                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
410                             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);
411 
412 
413                     VRR_K_g_s_p_s(
414                             PRIM_INT__g_s_p_s,
415                             PRIM_INT__g_s_s_s,
416                             PRIM_INT__f_s_s_s,
417                             Q_PA,
418                             aoq_PQ,
419                             one_over_2pq,
420                             5);
421 
422 
423                     ostei_general_vrr_K(5, 0, 2, 0, 4,
424                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
425                             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);
426 
427 
428                     ostei_general_vrr_K(6, 0, 3, 0, 3,
429                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
430                             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);
431 
432 
433                     VRR_K_f_s_p_s(
434                             PRIM_INT__f_s_p_s,
435                             PRIM_INT__f_s_s_s,
436                             PRIM_INT__d_s_s_s,
437                             Q_PA,
438                             aoq_PQ,
439                             one_over_2pq,
440                             5);
441 
442 
443                     ostei_general_vrr_K(4, 0, 2, 0, 4,
444                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
445                             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);
446 
447 
448                     ostei_general_vrr_K(5, 0, 3, 0, 3,
449                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
450                             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);
451 
452 
453                     ostei_general_vrr_K(6, 0, 4, 0, 2,
454                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
455                             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);
456 
457 
458 
459                     // Forming PRIM_INT__d_s_p_s[5 * 18];
460                     for(n = 0; n < 5; ++n)  // loop over orders of auxiliary function
461                     {
462 
463                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
464                         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]);
465                         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]);
466 
467                         PRIM_INT__d_s_p_s[n * 18 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
468                         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]);
469 
470                         PRIM_INT__d_s_p_s[n * 18 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
471                         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]);
472 
473                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 1]);
474                         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]);
475                         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]);
476 
477                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 1]);
478                         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]);
479                         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]);
480 
481                         PRIM_INT__d_s_p_s[n * 18 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
482                         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]);
483 
484                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 2]);
485                         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]);
486                         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]);
487 
488                         PRIM_INT__d_s_p_s[n * 18 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 2]);
489                         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]);
490 
491                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 2]);
492                         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]);
493                         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]);
494 
495                         PRIM_INT__d_s_p_s[n * 18 + 9] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
496                         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]);
497 
498                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
499                         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]);
500                         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]);
501 
502                         PRIM_INT__d_s_p_s[n * 18 + 11] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
503                         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]);
504 
505                         PRIM_INT__d_s_p_s[n * 18 + 12] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 4]);
506                         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]);
507 
508                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 4]);
509                         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]);
510                         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]);
511 
512                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 4]);
513                         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]);
514                         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]);
515 
516                         PRIM_INT__d_s_p_s[n * 18 + 15] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
517                         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]);
518 
519                         PRIM_INT__d_s_p_s[n * 18 + 16] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
520                         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]);
521 
522                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
523                         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]);
524                         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]);
525 
526                     }
527 
528 
529                     VRR_K_f_s_d_s(
530                             PRIM_INT__f_s_d_s,
531                             PRIM_INT__f_s_p_s,
532                             PRIM_INT__f_s_s_s,
533                             PRIM_INT__d_s_p_s,
534                             Q_PA,
535                             a_over_q,
536                             aoq_PQ,
537                             one_over_2pq,
538                             one_over_2q,
539                             4);
540 
541 
542                     ostei_general_vrr_K(4, 0, 3, 0, 3,
543                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
544                             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);
545 
546 
547                     ostei_general_vrr_K(5, 0, 4, 0, 2,
548                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
549                             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);
550 
551 
552                     ostei_general_vrr_K(6, 0, 5, 0, 1,
553                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
554                             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);
555 
556 
557                     ostei_general_vrr1_I(7, 7,
558                             one_over_2p, a_over_p, aop_PQ, P_PA,
559                             PRIM_INT__i_s_s_s, PRIM_INT__h_s_s_s, PRIM_INT__k_s_s_s);
560 
561 
562                     ostei_general_vrr_K(7, 0, 1, 0, 5,
563                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
564                             PRIM_INT__k_s_s_s, NULL, NULL, PRIM_INT__i_s_s_s, NULL, PRIM_INT__k_s_p_s);
565 
566 
567                     ostei_general_vrr_K(7, 0, 2, 0, 4,
568                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
569                             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);
570 
571 
572                     ostei_general_vrr_K(7, 0, 3, 0, 3,
573                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
574                             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);
575 
576 
577                     ostei_general_vrr_K(7, 0, 4, 0, 2,
578                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
579                             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);
580 
581 
582                     ostei_general_vrr_K(7, 0, 5, 0, 1,
583                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
584                             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);
585 
586 
587                     ostei_general_vrr1_I(8, 6,
588                             one_over_2p, a_over_p, aop_PQ, P_PA,
589                             PRIM_INT__k_s_s_s, PRIM_INT__i_s_s_s, PRIM_INT__l_s_s_s);
590 
591 
592                     ostei_general_vrr_K(8, 0, 1, 0, 5,
593                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
594                             PRIM_INT__l_s_s_s, NULL, NULL, PRIM_INT__k_s_s_s, NULL, PRIM_INT__l_s_p_s);
595 
596 
597                     ostei_general_vrr_K(8, 0, 2, 0, 4,
598                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
599                             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);
600 
601 
602                     ostei_general_vrr_K(8, 0, 3, 0, 3,
603                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
604                             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);
605 
606 
607                     ostei_general_vrr_K(8, 0, 4, 0, 2,
608                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
609                             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);
610 
611 
612                     ostei_general_vrr_K(8, 0, 5, 0, 1,
613                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
614                             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);
615 
616 
617 
618 
619                     ////////////////////////////////////
620                     // Accumulate contracted integrals
621                     ////////////////////////////////////
622                     if(lastoffset == 0)
623                     {
624                         contract_all(420, PRIM_INT__i_s_g_s, PRIM_PTR_INT__i_s_g_s);
625                         contract_all(588, PRIM_INT__i_s_h_s, PRIM_PTR_INT__i_s_h_s);
626                         contract_all(540, PRIM_INT__k_s_g_s, PRIM_PTR_INT__k_s_g_s);
627                         contract_all(756, PRIM_INT__k_s_h_s, PRIM_PTR_INT__k_s_h_s);
628                         contract_all(675, PRIM_INT__l_s_g_s, PRIM_PTR_INT__l_s_g_s);
629                         contract_all(945, PRIM_INT__l_s_h_s, PRIM_PTR_INT__l_s_h_s);
630                     }
631                     else
632                     {
633                         contract(420, shelloffsets, PRIM_INT__i_s_g_s, PRIM_PTR_INT__i_s_g_s);
634                         contract(588, shelloffsets, PRIM_INT__i_s_h_s, PRIM_PTR_INT__i_s_h_s);
635                         contract(540, shelloffsets, PRIM_INT__k_s_g_s, PRIM_PTR_INT__k_s_g_s);
636                         contract(756, shelloffsets, PRIM_INT__k_s_h_s, PRIM_PTR_INT__k_s_h_s);
637                         contract(675, shelloffsets, PRIM_INT__l_s_g_s, PRIM_PTR_INT__l_s_g_s);
638                         contract(945, shelloffsets, PRIM_INT__l_s_h_s, PRIM_PTR_INT__l_s_h_s);
639                         PRIM_PTR_INT__i_s_g_s += lastoffset*420;
640                         PRIM_PTR_INT__i_s_h_s += lastoffset*588;
641                         PRIM_PTR_INT__k_s_g_s += lastoffset*540;
642                         PRIM_PTR_INT__k_s_h_s += lastoffset*756;
643                         PRIM_PTR_INT__l_s_g_s += lastoffset*675;
644                         PRIM_PTR_INT__l_s_h_s += lastoffset*945;
645                     }
646 
647                 }  // close loop over j
648             }  // close loop over i
649 
650             //Advance to the next batch
651             jstart = SIMINT_SIMD_ROUND(jend);
652 
653             //////////////////////////////////////////////
654             // Contracted integrals: Horizontal recurrance
655             //////////////////////////////////////////////
656 
657 
658             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
659 
660 
661             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
662             {
663                 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
664 
665                 // set up HRR pointers
666                 double const * restrict HRR_INT__i_s_g_s = INT__i_s_g_s + abcd * 420;
667                 double const * restrict HRR_INT__i_s_h_s = INT__i_s_h_s + abcd * 588;
668                 double const * restrict HRR_INT__k_s_g_s = INT__k_s_g_s + abcd * 540;
669                 double const * restrict HRR_INT__k_s_h_s = INT__k_s_h_s + abcd * 756;
670                 double const * restrict HRR_INT__l_s_g_s = INT__l_s_g_s + abcd * 675;
671                 double const * restrict HRR_INT__l_s_h_s = INT__l_s_h_s + abcd * 945;
672                 double * restrict HRR_INT__i_d_g_p = INT__i_d_g_p + real_abcd * 7560;
673 
674                 // form INT__i_p_g_s
675                 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);
676 
677                 // form INT__i_p_h_s
678                 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);
679 
680                 // form INT__k_p_g_s
681                 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);
682 
683                 // form INT__k_p_h_s
684                 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);
685 
686                 // form INT__i_d_g_s
687                 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);
688 
689                 // form INT__i_d_h_s
690                 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);
691 
692                 // form INT__i_d_g_p
693                 HRR_L_g_p(
694                     HRR_INT__i_d_g_p,
695                     HRR_INT__i_d_g_s,
696                     HRR_INT__i_d_h_s,
697                     hCD, 168);
698 
699 
700             }  // close HRR loop
701 
702 
703         }   // close loop cdbatch
704 
705         istart = iend;
706     }  // close loop over ab
707 
708     return P.nshell12_clip * Q.nshell12_clip;
709 }
710 
ostei_d_i_g_p(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__d_i_g_p)711 int ostei_d_i_g_p(struct simint_multi_shellpair const P,
712                   struct simint_multi_shellpair const Q,
713                   double screen_tol,
714                   double * const restrict work,
715                   double * const restrict INT__d_i_g_p)
716 {
717     double P_AB[3*P.nshell12];
718     struct simint_multi_shellpair P_tmp = P;
719     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
720     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
721     P_tmp.AB_x = P_AB;
722     P_tmp.AB_y = P_AB + P.nshell12;
723     P_tmp.AB_z = P_AB + 2*P.nshell12;
724 
725     for(int i = 0; i < P.nshell12; i++)
726     {
727         P_tmp.AB_x[i] = -P.AB_x[i];
728         P_tmp.AB_y[i] = -P.AB_y[i];
729         P_tmp.AB_z[i] = -P.AB_z[i];
730     }
731 
732     int ret = ostei_i_d_g_p(P_tmp, Q, screen_tol, work, INT__d_i_g_p);
733     double buffer[7560] SIMINT_ALIGN_ARRAY_DBL;
734 
735     for(int q = 0; q < ret; q++)
736     {
737         int idx = 0;
738         for(int a = 0; a < 6; ++a)
739         for(int b = 0; b < 28; ++b)
740         for(int c = 0; c < 15; ++c)
741         for(int d = 0; d < 3; ++d)
742             buffer[idx++] = INT__d_i_g_p[q*7560+b*270+a*45+c*3+d];
743 
744         memcpy(INT__d_i_g_p+q*7560, buffer, 7560*sizeof(double));
745     }
746 
747     return ret;
748 }
749 
ostei_i_d_p_g(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__i_d_p_g)750 int ostei_i_d_p_g(struct simint_multi_shellpair const P,
751                   struct simint_multi_shellpair const Q,
752                   double screen_tol,
753                   double * const restrict work,
754                   double * const restrict INT__i_d_p_g)
755 {
756     double Q_AB[3*Q.nshell12];
757     struct simint_multi_shellpair Q_tmp = Q;
758     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
759     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
760     Q_tmp.AB_x = Q_AB;
761     Q_tmp.AB_y = Q_AB + Q.nshell12;
762     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
763 
764     for(int i = 0; i < Q.nshell12; i++)
765     {
766         Q_tmp.AB_x[i] = -Q.AB_x[i];
767         Q_tmp.AB_y[i] = -Q.AB_y[i];
768         Q_tmp.AB_z[i] = -Q.AB_z[i];
769     }
770 
771     int ret = ostei_i_d_g_p(P, Q_tmp, screen_tol, work, INT__i_d_p_g);
772     double buffer[7560] SIMINT_ALIGN_ARRAY_DBL;
773 
774     for(int q = 0; q < ret; q++)
775     {
776         int idx = 0;
777         for(int a = 0; a < 28; ++a)
778         for(int b = 0; b < 6; ++b)
779         for(int c = 0; c < 3; ++c)
780         for(int d = 0; d < 15; ++d)
781             buffer[idx++] = INT__i_d_p_g[q*7560+a*270+b*45+d*3+c];
782 
783         memcpy(INT__i_d_p_g+q*7560, buffer, 7560*sizeof(double));
784     }
785 
786     return ret;
787 }
788 
ostei_d_i_p_g(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__d_i_p_g)789 int ostei_d_i_p_g(struct simint_multi_shellpair const P,
790                   struct simint_multi_shellpair const Q,
791                   double screen_tol,
792                   double * const restrict work,
793                   double * const restrict INT__d_i_p_g)
794 {
795     double P_AB[3*P.nshell12];
796     struct simint_multi_shellpair P_tmp = P;
797     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
798     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
799     P_tmp.AB_x = P_AB;
800     P_tmp.AB_y = P_AB + P.nshell12;
801     P_tmp.AB_z = P_AB + 2*P.nshell12;
802 
803     for(int i = 0; i < P.nshell12; i++)
804     {
805         P_tmp.AB_x[i] = -P.AB_x[i];
806         P_tmp.AB_y[i] = -P.AB_y[i];
807         P_tmp.AB_z[i] = -P.AB_z[i];
808     }
809 
810     double Q_AB[3*Q.nshell12];
811     struct simint_multi_shellpair Q_tmp = Q;
812     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
813     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
814     Q_tmp.AB_x = Q_AB;
815     Q_tmp.AB_y = Q_AB + Q.nshell12;
816     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
817 
818     for(int i = 0; i < Q.nshell12; i++)
819     {
820         Q_tmp.AB_x[i] = -Q.AB_x[i];
821         Q_tmp.AB_y[i] = -Q.AB_y[i];
822         Q_tmp.AB_z[i] = -Q.AB_z[i];
823     }
824 
825     int ret = ostei_i_d_g_p(P_tmp, Q_tmp, screen_tol, work, INT__d_i_p_g);
826     double buffer[7560] SIMINT_ALIGN_ARRAY_DBL;
827 
828     for(int q = 0; q < ret; q++)
829     {
830         int idx = 0;
831         for(int a = 0; a < 6; ++a)
832         for(int b = 0; b < 28; ++b)
833         for(int c = 0; c < 3; ++c)
834         for(int d = 0; d < 15; ++d)
835             buffer[idx++] = INT__d_i_p_g[q*7560+b*270+a*45+d*3+c];
836 
837         memcpy(INT__d_i_p_g+q*7560, buffer, 7560*sizeof(double));
838     }
839 
840     return ret;
841 }
842 
843