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