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