1 #include "simint/boys/boys.h"
2 #include "simint/ostei/gen/ostei_generated.h"
3 #include "simint/vectorization/vectorization.h"
4 #include <math.h>
5 #include <string.h>
6 
7 
ostei_h_s_k_i(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__h_s_k_i)8 int ostei_h_s_k_i(struct simint_multi_shellpair const P,
9                   struct simint_multi_shellpair const Q,
10                   double screen_tol,
11                   double * const restrict work,
12                   double * const restrict INT__h_s_k_i)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__h_s_k_i);
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__h_s_k_s = work + (SIMINT_NSHELL_SIMD * 0);
29     double * const INT__h_s_l_s = work + (SIMINT_NSHELL_SIMD * 756);
30     double * const INT__h_s_m_s = work + (SIMINT_NSHELL_SIMD * 1701);
31     double * const INT__h_s_n_s = work + (SIMINT_NSHELL_SIMD * 2856);
32     double * const INT__h_s_o_s = work + (SIMINT_NSHELL_SIMD * 4242);
33     double * const INT__h_s_q_s = work + (SIMINT_NSHELL_SIMD * 5880);
34     double * const INT__h_s_r_s = work + (SIMINT_NSHELL_SIMD * 7791);
35     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*9996);
36     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
37     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_p_s = primwork + 19;
38     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_d_s = primwork + 73;
39     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_f_s = primwork + 175;
40     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_g_s = primwork + 335;
41     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_h_s = primwork + 560;
42     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_i_s = primwork + 854;
43     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_k_s = primwork + 1218;
44     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_l_s = primwork + 1650;
45     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_m_s = primwork + 2145;
46     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_n_s = primwork + 2695;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_o_s = primwork + 3289;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_q_s = primwork + 3913;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_r_s = primwork + 4550;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_f_s = primwork + 5180;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_g_s = primwork + 5330;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_h_s = primwork + 5555;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_i_s = primwork + 5870;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_k_s = primwork + 6290;
55     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_l_s = primwork + 6830;
56     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_m_s = primwork + 7505;
57     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_n_s = primwork + 8330;
58     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_o_s = primwork + 9320;
59     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_q_s = primwork + 10490;
60     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_r_s = primwork + 11855;
61     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_g_s = primwork + 13430;
62     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_h_s = primwork + 13790;
63     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_i_s = primwork + 14294;
64     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_k_s = primwork + 14966;
65     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_l_s = primwork + 15830;
66     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_m_s = primwork + 16910;
67     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_n_s = primwork + 18230;
68     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_o_s = primwork + 19814;
69     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_q_s = primwork + 21686;
70     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_r_s = primwork + 23870;
71     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_h_s = primwork + 26390;
72     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_i_s = primwork + 27020;
73     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_k_s = primwork + 27860;
74     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_l_s = primwork + 28940;
75     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_m_s = primwork + 30290;
76     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_n_s = primwork + 31940;
77     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_o_s = primwork + 33920;
78     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_q_s = primwork + 36260;
79     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_r_s = primwork + 38990;
80     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_i_s = primwork + 42140;
81     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_k_s = primwork + 42980;
82     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_l_s = primwork + 44060;
83     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_m_s = primwork + 45410;
84     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_n_s = primwork + 47060;
85     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_o_s = primwork + 49040;
86     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_q_s = primwork + 51380;
87     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_r_s = primwork + 54110;
88     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_k_s = primwork + 57260;
89     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_l_s = primwork + 58016;
90     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_m_s = primwork + 58961;
91     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_n_s = primwork + 60116;
92     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_o_s = primwork + 61502;
93     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_q_s = primwork + 63140;
94     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_r_s = primwork + 65051;
95     double * const hrrwork = (double *)(primwork + 67256);
96     double * const HRR_INT__h_s_k_p = hrrwork + 0;
97     double * const HRR_INT__h_s_k_d = hrrwork + 2268;
98     double * const HRR_INT__h_s_k_f = hrrwork + 6804;
99     double * const HRR_INT__h_s_k_g = hrrwork + 14364;
100     double * const HRR_INT__h_s_k_h = hrrwork + 25704;
101     double * const HRR_INT__h_s_l_p = hrrwork + 41580;
102     double * const HRR_INT__h_s_l_d = hrrwork + 44415;
103     double * const HRR_INT__h_s_l_f = hrrwork + 50085;
104     double * const HRR_INT__h_s_l_g = hrrwork + 59535;
105     double * const HRR_INT__h_s_l_h = hrrwork + 73710;
106     double * const HRR_INT__h_s_m_p = hrrwork + 93555;
107     double * const HRR_INT__h_s_m_d = hrrwork + 97020;
108     double * const HRR_INT__h_s_m_f = hrrwork + 103950;
109     double * const HRR_INT__h_s_m_g = hrrwork + 115500;
110     double * const HRR_INT__h_s_n_p = hrrwork + 132825;
111     double * const HRR_INT__h_s_n_d = hrrwork + 136983;
112     double * const HRR_INT__h_s_n_f = hrrwork + 145299;
113     double * const HRR_INT__h_s_o_p = hrrwork + 159159;
114     double * const HRR_INT__h_s_o_d = hrrwork + 164073;
115     double * const HRR_INT__h_s_q_p = hrrwork + 173901;
116 
117 
118     // Create constants
119     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
120     const SIMINT_DBLTYPE const_10 = SIMINT_DBLSET1(10);
121     const SIMINT_DBLTYPE const_11 = SIMINT_DBLSET1(11);
122     const SIMINT_DBLTYPE const_12 = SIMINT_DBLSET1(12);
123     const SIMINT_DBLTYPE const_13 = SIMINT_DBLSET1(13);
124     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
125     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
126     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
127     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
128     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
129     const SIMINT_DBLTYPE const_7 = SIMINT_DBLSET1(7);
130     const SIMINT_DBLTYPE const_8 = SIMINT_DBLSET1(8);
131     const SIMINT_DBLTYPE const_9 = SIMINT_DBLSET1(9);
132     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
133 
134 
135     ////////////////////////////////////////
136     // Loop over shells and primitives
137     ////////////////////////////////////////
138 
139     real_abcd = 0;
140     istart = 0;
141     for(ab = 0; ab < P.nshell12_clip; ++ab)
142     {
143         const int iend = istart + P.nprim12[ab];
144 
145         cd = 0;
146         jstart = 0;
147 
148         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
149         {
150             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
151             int jend = jstart;
152             for(i = 0; i < nshellbatch; i++)
153                 jend += Q.nprim12[cd+i];
154 
155             // Clear the beginning of the workspace (where we are accumulating integrals)
156             memset(work, 0, SIMINT_NSHELL_SIMD * 9996 * sizeof(double));
157             abcd = 0;
158 
159 
160             for(i = istart; i < iend; ++i)
161             {
162                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
163 
164                 if(check_screen)
165                 {
166                     // Skip this whole thing if always insignificant
167                     if((P.screen[i] * Q.screen_max) < screen_tol)
168                         continue;
169                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
170                 }
171 
172                 icd = 0;
173                 iprimcd = 0;
174                 nprim_icd = Q.nprim12[cd];
175                 double * restrict PRIM_PTR_INT__h_s_k_s = INT__h_s_k_s + abcd * 756;
176                 double * restrict PRIM_PTR_INT__h_s_l_s = INT__h_s_l_s + abcd * 945;
177                 double * restrict PRIM_PTR_INT__h_s_m_s = INT__h_s_m_s + abcd * 1155;
178                 double * restrict PRIM_PTR_INT__h_s_n_s = INT__h_s_n_s + abcd * 1386;
179                 double * restrict PRIM_PTR_INT__h_s_o_s = INT__h_s_o_s + abcd * 1638;
180                 double * restrict PRIM_PTR_INT__h_s_q_s = INT__h_s_q_s + abcd * 1911;
181                 double * restrict PRIM_PTR_INT__h_s_r_s = INT__h_s_r_s + abcd * 2205;
182 
183 
184 
185                 // Load these one per loop over i
186                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
187                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
188                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
189 
190                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
191 
192                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
193                 {
194                     // calculate the shell offsets
195                     // these are the offset from the shell pointed to by cd
196                     // for each element
197                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
198                     int lastoffset = 0;
199                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
200 
201                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
202                     {
203                         // Handle if the first element of the vector is a new shell
204                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
205                         {
206                             nprim_icd += Q.nprim12[cd + (++icd)];
207                             PRIM_PTR_INT__h_s_k_s += 756;
208                             PRIM_PTR_INT__h_s_l_s += 945;
209                             PRIM_PTR_INT__h_s_m_s += 1155;
210                             PRIM_PTR_INT__h_s_n_s += 1386;
211                             PRIM_PTR_INT__h_s_o_s += 1638;
212                             PRIM_PTR_INT__h_s_q_s += 1911;
213                             PRIM_PTR_INT__h_s_r_s += 2205;
214                         }
215                         iprimcd++;
216                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
217                         {
218                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
219                             {
220                                 shelloffsets[n] = shelloffsets[n-1] + 1;
221                                 lastoffset++;
222                                 nprim_icd += Q.nprim12[cd + (++icd)];
223                             }
224                             else
225                                 shelloffsets[n] = shelloffsets[n-1];
226                             iprimcd++;
227                         }
228                     }
229                     else
230                         iprimcd += SIMINT_SIMD_LEN;
231 
232                     // Do we have to compute this vector (or has it been screened out)?
233                     // (not_screened != 0 means we have to do this vector)
234                     if(check_screen)
235                     {
236                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
237                         if(vmax < screen_tol)
238                         {
239                             PRIM_PTR_INT__h_s_k_s += lastoffset*756;
240                             PRIM_PTR_INT__h_s_l_s += lastoffset*945;
241                             PRIM_PTR_INT__h_s_m_s += lastoffset*1155;
242                             PRIM_PTR_INT__h_s_n_s += lastoffset*1386;
243                             PRIM_PTR_INT__h_s_o_s += lastoffset*1638;
244                             PRIM_PTR_INT__h_s_q_s += lastoffset*1911;
245                             PRIM_PTR_INT__h_s_r_s += lastoffset*2205;
246                             continue;
247                         }
248                     }
249 
250                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
251                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
252                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
253                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
254 
255 
256                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
257                     SIMINT_DBLTYPE PQ[3];
258                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
259                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
260                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
261                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
262                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
263                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
264 
265                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
266                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
267                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
268                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
269                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
270                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
271                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
272 
273                     // NOTE: Minus sign!
274                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
275                     SIMINT_DBLTYPE aop_PQ[3];
276                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
277                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
278                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
279 
280                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
281                     SIMINT_DBLTYPE aoq_PQ[3];
282                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
283                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
284                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
285                     // Put a minus sign here so we don't have to in RR routines
286                     a_over_q = SIMINT_NEG(a_over_q);
287 
288 
289                     //////////////////////////////////////////////
290                     // Fjt function section
291                     // Maximum v value: 18
292                     //////////////////////////////////////////////
293                     // The parameter to the Fjt function
294                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
295 
296 
297                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
298 
299 
300                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 18);
301                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
302                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
303                     for(n = 0; n <= 18; n++)
304                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
305 
306                     //////////////////////////////////////////////
307                     // Primitive integrals: Vertical recurrance
308                     //////////////////////////////////////////////
309 
310                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
311                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
312                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
313                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
314                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
315                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
316                     const SIMINT_DBLTYPE vrr_const_3_over_2q = SIMINT_MUL(const_3, one_over_2q);
317                     const SIMINT_DBLTYPE vrr_const_4_over_2q = SIMINT_MUL(const_4, one_over_2q);
318                     const SIMINT_DBLTYPE vrr_const_5_over_2q = SIMINT_MUL(const_5, one_over_2q);
319                     const SIMINT_DBLTYPE vrr_const_6_over_2q = SIMINT_MUL(const_6, one_over_2q);
320                     const SIMINT_DBLTYPE vrr_const_7_over_2q = SIMINT_MUL(const_7, one_over_2q);
321                     const SIMINT_DBLTYPE vrr_const_8_over_2q = SIMINT_MUL(const_8, one_over_2q);
322                     const SIMINT_DBLTYPE vrr_const_9_over_2q = SIMINT_MUL(const_9, one_over_2q);
323                     const SIMINT_DBLTYPE vrr_const_10_over_2q = SIMINT_MUL(const_10, one_over_2q);
324                     const SIMINT_DBLTYPE vrr_const_11_over_2q = SIMINT_MUL(const_11, one_over_2q);
325                     const SIMINT_DBLTYPE vrr_const_12_over_2q = SIMINT_MUL(const_12, one_over_2q);
326                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
327                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
328                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
329                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
330                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
331                     const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
332                     const SIMINT_DBLTYPE vrr_const_7_over_2pq = SIMINT_MUL(const_7, one_over_2pq);
333                     const SIMINT_DBLTYPE vrr_const_8_over_2pq = SIMINT_MUL(const_8, one_over_2pq);
334                     const SIMINT_DBLTYPE vrr_const_9_over_2pq = SIMINT_MUL(const_9, one_over_2pq);
335                     const SIMINT_DBLTYPE vrr_const_10_over_2pq = SIMINT_MUL(const_10, one_over_2pq);
336                     const SIMINT_DBLTYPE vrr_const_11_over_2pq = SIMINT_MUL(const_11, one_over_2pq);
337                     const SIMINT_DBLTYPE vrr_const_12_over_2pq = SIMINT_MUL(const_12, one_over_2pq);
338                     const SIMINT_DBLTYPE vrr_const_13_over_2pq = SIMINT_MUL(const_13, one_over_2pq);
339 
340 
341 
342                     // Forming PRIM_INT__s_s_p_s[18 * 3];
343                     for(n = 0; n < 18; ++n)  // loop over orders of auxiliary function
344                     {
345 
346                         PRIM_INT__s_s_p_s[n * 3 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
347                         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]);
348 
349                         PRIM_INT__s_s_p_s[n * 3 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
350                         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]);
351 
352                         PRIM_INT__s_s_p_s[n * 3 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
353                         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]);
354 
355                     }
356 
357 
358 
359                     // Forming PRIM_INT__s_s_d_s[17 * 6];
360                     for(n = 0; n < 17; ++n)  // loop over orders of auxiliary function
361                     {
362 
363                         PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_p_s[n * 3 + 0]);
364                         PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_d_s[n * 6 + 0]);
365                         PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__s_s_d_s[n * 6 + 0]);
366 
367                         PRIM_INT__s_s_d_s[n * 6 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_p_s[n * 3 + 0]);
368                         PRIM_INT__s_s_d_s[n * 6 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_d_s[n * 6 + 1]);
369 
370                         PRIM_INT__s_s_d_s[n * 6 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 0]);
371                         PRIM_INT__s_s_d_s[n * 6 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_d_s[n * 6 + 2]);
372 
373                         PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_p_s[n * 3 + 1]);
374                         PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_d_s[n * 6 + 3]);
375                         PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__s_s_d_s[n * 6 + 3]);
376 
377                         PRIM_INT__s_s_d_s[n * 6 + 4] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 1]);
378                         PRIM_INT__s_s_d_s[n * 6 + 4] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_d_s[n * 6 + 4]);
379 
380                         PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 2]);
381                         PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__s_s_d_s[n * 6 + 5]);
382                         PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__s_s_d_s[n * 6 + 5]);
383 
384                     }
385 
386 
387 
388                     // Forming PRIM_INT__s_s_f_s[16 * 10];
389                     for(n = 0; n < 16; ++n)  // loop over orders of auxiliary function
390                     {
391 
392                         PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 0]);
393                         PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__s_s_f_s[n * 10 + 0]);
394                         PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_p_s[n * 3 + 0]), PRIM_INT__s_s_f_s[n * 10 + 0]);
395 
396                         PRIM_INT__s_s_f_s[n * 10 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 0]);
397                         PRIM_INT__s_s_f_s[n * 10 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__s_s_f_s[n * 10 + 1]);
398 
399                         PRIM_INT__s_s_f_s[n * 10 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 0]);
400                         PRIM_INT__s_s_f_s[n * 10 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__s_s_f_s[n * 10 + 2]);
401 
402                         PRIM_INT__s_s_f_s[n * 10 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 3]);
403                         PRIM_INT__s_s_f_s[n * 10 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__s_s_f_s[n * 10 + 3]);
404 
405                         PRIM_INT__s_s_f_s[n * 10 + 4] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 1]);
406                         PRIM_INT__s_s_f_s[n * 10 + 4] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__s_s_f_s[n * 10 + 4]);
407 
408                         PRIM_INT__s_s_f_s[n * 10 + 5] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 5]);
409                         PRIM_INT__s_s_f_s[n * 10 + 5] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__s_s_f_s[n * 10 + 5]);
410 
411                         PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 3]);
412                         PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__s_s_f_s[n * 10 + 6]);
413                         PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_p_s[n * 3 + 1]), PRIM_INT__s_s_f_s[n * 10 + 6]);
414 
415                         PRIM_INT__s_s_f_s[n * 10 + 7] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 3]);
416                         PRIM_INT__s_s_f_s[n * 10 + 7] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__s_s_f_s[n * 10 + 7]);
417 
418                         PRIM_INT__s_s_f_s[n * 10 + 8] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 5]);
419                         PRIM_INT__s_s_f_s[n * 10 + 8] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__s_s_f_s[n * 10 + 8]);
420 
421                         PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 5]);
422                         PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__s_s_f_s[n * 10 + 9]);
423                         PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__s_s_p_s[n * 3 + 2]), PRIM_INT__s_s_f_s[n * 10 + 9]);
424 
425                     }
426 
427 
428                     VRR_K_s_s_g_s(
429                             PRIM_INT__s_s_g_s,
430                             PRIM_INT__s_s_f_s,
431                             PRIM_INT__s_s_d_s,
432                             Q_PA,
433                             a_over_q,
434                             aoq_PQ,
435                             one_over_2q,
436                             15);
437 
438 
439                     VRR_K_s_s_h_s(
440                             PRIM_INT__s_s_h_s,
441                             PRIM_INT__s_s_g_s,
442                             PRIM_INT__s_s_f_s,
443                             Q_PA,
444                             a_over_q,
445                             aoq_PQ,
446                             one_over_2q,
447                             14);
448 
449 
450                     ostei_general_vrr1_K(6, 13,
451                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
452                             PRIM_INT__s_s_h_s, PRIM_INT__s_s_g_s, PRIM_INT__s_s_i_s);
453 
454 
455                     ostei_general_vrr1_K(7, 12,
456                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
457                             PRIM_INT__s_s_i_s, PRIM_INT__s_s_h_s, PRIM_INT__s_s_k_s);
458 
459 
460                     ostei_general_vrr_I(1, 0, 7, 0, 5,
461                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
462                             PRIM_INT__s_s_k_s, NULL, NULL, PRIM_INT__s_s_i_s, NULL, PRIM_INT__p_s_k_s);
463 
464 
465                     ostei_general_vrr_I(1, 0, 6, 0, 5,
466                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
467                             PRIM_INT__s_s_i_s, NULL, NULL, PRIM_INT__s_s_h_s, NULL, PRIM_INT__p_s_i_s);
468 
469 
470                     ostei_general_vrr_I(2, 0, 7, 0, 4,
471                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
472                             PRIM_INT__p_s_k_s, PRIM_INT__s_s_k_s, NULL, PRIM_INT__p_s_i_s, NULL, PRIM_INT__d_s_k_s);
473 
474 
475                     ostei_general_vrr_I(1, 0, 5, 0, 5,
476                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
477                             PRIM_INT__s_s_h_s, NULL, NULL, PRIM_INT__s_s_g_s, NULL, PRIM_INT__p_s_h_s);
478 
479 
480                     ostei_general_vrr_I(2, 0, 6, 0, 4,
481                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
482                             PRIM_INT__p_s_i_s, PRIM_INT__s_s_i_s, NULL, PRIM_INT__p_s_h_s, NULL, PRIM_INT__d_s_i_s);
483 
484 
485                     ostei_general_vrr_I(3, 0, 7, 0, 3,
486                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
487                             PRIM_INT__d_s_k_s, PRIM_INT__p_s_k_s, NULL, PRIM_INT__d_s_i_s, NULL, PRIM_INT__f_s_k_s);
488 
489 
490                     VRR_I_p_s_g_s(
491                             PRIM_INT__p_s_g_s,
492                             PRIM_INT__s_s_g_s,
493                             PRIM_INT__s_s_f_s,
494                             P_PA,
495                             aop_PQ,
496                             one_over_2pq,
497                             5);
498 
499 
500                     ostei_general_vrr_I(2, 0, 5, 0, 4,
501                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
502                             PRIM_INT__p_s_h_s, PRIM_INT__s_s_h_s, NULL, PRIM_INT__p_s_g_s, NULL, PRIM_INT__d_s_h_s);
503 
504 
505                     ostei_general_vrr_I(3, 0, 6, 0, 3,
506                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
507                             PRIM_INT__d_s_i_s, PRIM_INT__p_s_i_s, NULL, PRIM_INT__d_s_h_s, NULL, PRIM_INT__f_s_i_s);
508 
509 
510                     ostei_general_vrr_I(4, 0, 7, 0, 2,
511                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
512                             PRIM_INT__f_s_k_s, PRIM_INT__d_s_k_s, NULL, PRIM_INT__f_s_i_s, NULL, PRIM_INT__g_s_k_s);
513 
514 
515                     VRR_I_p_s_f_s(
516                             PRIM_INT__p_s_f_s,
517                             PRIM_INT__s_s_f_s,
518                             PRIM_INT__s_s_d_s,
519                             P_PA,
520                             aop_PQ,
521                             one_over_2pq,
522                             5);
523 
524 
525                     ostei_general_vrr_I(2, 0, 4, 0, 4,
526                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
527                             PRIM_INT__p_s_g_s, PRIM_INT__s_s_g_s, NULL, PRIM_INT__p_s_f_s, NULL, PRIM_INT__d_s_g_s);
528 
529 
530                     ostei_general_vrr_I(3, 0, 5, 0, 3,
531                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
532                             PRIM_INT__d_s_h_s, PRIM_INT__p_s_h_s, NULL, PRIM_INT__d_s_g_s, NULL, PRIM_INT__f_s_h_s);
533 
534 
535                     ostei_general_vrr_I(4, 0, 6, 0, 2,
536                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
537                             PRIM_INT__f_s_i_s, PRIM_INT__d_s_i_s, NULL, PRIM_INT__f_s_h_s, NULL, PRIM_INT__g_s_i_s);
538 
539 
540                     ostei_general_vrr_I(5, 0, 7, 0, 1,
541                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
542                             PRIM_INT__g_s_k_s, PRIM_INT__f_s_k_s, NULL, PRIM_INT__g_s_i_s, NULL, PRIM_INT__h_s_k_s);
543 
544 
545                     ostei_general_vrr1_K(8, 11,
546                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
547                             PRIM_INT__s_s_k_s, PRIM_INT__s_s_i_s, PRIM_INT__s_s_l_s);
548 
549 
550                     ostei_general_vrr_I(1, 0, 8, 0, 5,
551                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
552                             PRIM_INT__s_s_l_s, NULL, NULL, PRIM_INT__s_s_k_s, NULL, PRIM_INT__p_s_l_s);
553 
554 
555                     ostei_general_vrr_I(2, 0, 8, 0, 4,
556                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
557                             PRIM_INT__p_s_l_s, PRIM_INT__s_s_l_s, NULL, PRIM_INT__p_s_k_s, NULL, PRIM_INT__d_s_l_s);
558 
559 
560                     ostei_general_vrr_I(3, 0, 8, 0, 3,
561                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
562                             PRIM_INT__d_s_l_s, PRIM_INT__p_s_l_s, NULL, PRIM_INT__d_s_k_s, NULL, PRIM_INT__f_s_l_s);
563 
564 
565                     ostei_general_vrr_I(4, 0, 8, 0, 2,
566                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
567                             PRIM_INT__f_s_l_s, PRIM_INT__d_s_l_s, NULL, PRIM_INT__f_s_k_s, NULL, PRIM_INT__g_s_l_s);
568 
569 
570                     ostei_general_vrr_I(5, 0, 8, 0, 1,
571                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
572                             PRIM_INT__g_s_l_s, PRIM_INT__f_s_l_s, NULL, PRIM_INT__g_s_k_s, NULL, PRIM_INT__h_s_l_s);
573 
574 
575                     ostei_general_vrr1_K(9, 10,
576                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
577                             PRIM_INT__s_s_l_s, PRIM_INT__s_s_k_s, PRIM_INT__s_s_m_s);
578 
579 
580                     ostei_general_vrr_I(1, 0, 9, 0, 5,
581                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
582                             PRIM_INT__s_s_m_s, NULL, NULL, PRIM_INT__s_s_l_s, NULL, PRIM_INT__p_s_m_s);
583 
584 
585                     ostei_general_vrr_I(2, 0, 9, 0, 4,
586                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
587                             PRIM_INT__p_s_m_s, PRIM_INT__s_s_m_s, NULL, PRIM_INT__p_s_l_s, NULL, PRIM_INT__d_s_m_s);
588 
589 
590                     ostei_general_vrr_I(3, 0, 9, 0, 3,
591                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
592                             PRIM_INT__d_s_m_s, PRIM_INT__p_s_m_s, NULL, PRIM_INT__d_s_l_s, NULL, PRIM_INT__f_s_m_s);
593 
594 
595                     ostei_general_vrr_I(4, 0, 9, 0, 2,
596                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
597                             PRIM_INT__f_s_m_s, PRIM_INT__d_s_m_s, NULL, PRIM_INT__f_s_l_s, NULL, PRIM_INT__g_s_m_s);
598 
599 
600                     ostei_general_vrr_I(5, 0, 9, 0, 1,
601                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
602                             PRIM_INT__g_s_m_s, PRIM_INT__f_s_m_s, NULL, PRIM_INT__g_s_l_s, NULL, PRIM_INT__h_s_m_s);
603 
604 
605                     ostei_general_vrr1_K(10, 9,
606                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
607                             PRIM_INT__s_s_m_s, PRIM_INT__s_s_l_s, PRIM_INT__s_s_n_s);
608 
609 
610                     ostei_general_vrr_I(1, 0, 10, 0, 5,
611                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
612                             PRIM_INT__s_s_n_s, NULL, NULL, PRIM_INT__s_s_m_s, NULL, PRIM_INT__p_s_n_s);
613 
614 
615                     ostei_general_vrr_I(2, 0, 10, 0, 4,
616                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
617                             PRIM_INT__p_s_n_s, PRIM_INT__s_s_n_s, NULL, PRIM_INT__p_s_m_s, NULL, PRIM_INT__d_s_n_s);
618 
619 
620                     ostei_general_vrr_I(3, 0, 10, 0, 3,
621                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
622                             PRIM_INT__d_s_n_s, PRIM_INT__p_s_n_s, NULL, PRIM_INT__d_s_m_s, NULL, PRIM_INT__f_s_n_s);
623 
624 
625                     ostei_general_vrr_I(4, 0, 10, 0, 2,
626                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
627                             PRIM_INT__f_s_n_s, PRIM_INT__d_s_n_s, NULL, PRIM_INT__f_s_m_s, NULL, PRIM_INT__g_s_n_s);
628 
629 
630                     ostei_general_vrr_I(5, 0, 10, 0, 1,
631                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
632                             PRIM_INT__g_s_n_s, PRIM_INT__f_s_n_s, NULL, PRIM_INT__g_s_m_s, NULL, PRIM_INT__h_s_n_s);
633 
634 
635                     ostei_general_vrr1_K(11, 8,
636                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
637                             PRIM_INT__s_s_n_s, PRIM_INT__s_s_m_s, PRIM_INT__s_s_o_s);
638 
639 
640                     ostei_general_vrr_I(1, 0, 11, 0, 5,
641                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
642                             PRIM_INT__s_s_o_s, NULL, NULL, PRIM_INT__s_s_n_s, NULL, PRIM_INT__p_s_o_s);
643 
644 
645                     ostei_general_vrr_I(2, 0, 11, 0, 4,
646                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
647                             PRIM_INT__p_s_o_s, PRIM_INT__s_s_o_s, NULL, PRIM_INT__p_s_n_s, NULL, PRIM_INT__d_s_o_s);
648 
649 
650                     ostei_general_vrr_I(3, 0, 11, 0, 3,
651                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
652                             PRIM_INT__d_s_o_s, PRIM_INT__p_s_o_s, NULL, PRIM_INT__d_s_n_s, NULL, PRIM_INT__f_s_o_s);
653 
654 
655                     ostei_general_vrr_I(4, 0, 11, 0, 2,
656                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
657                             PRIM_INT__f_s_o_s, PRIM_INT__d_s_o_s, NULL, PRIM_INT__f_s_n_s, NULL, PRIM_INT__g_s_o_s);
658 
659 
660                     ostei_general_vrr_I(5, 0, 11, 0, 1,
661                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
662                             PRIM_INT__g_s_o_s, PRIM_INT__f_s_o_s, NULL, PRIM_INT__g_s_n_s, NULL, PRIM_INT__h_s_o_s);
663 
664 
665                     ostei_general_vrr1_K(12, 7,
666                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
667                             PRIM_INT__s_s_o_s, PRIM_INT__s_s_n_s, PRIM_INT__s_s_q_s);
668 
669 
670                     ostei_general_vrr_I(1, 0, 12, 0, 5,
671                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
672                             PRIM_INT__s_s_q_s, NULL, NULL, PRIM_INT__s_s_o_s, NULL, PRIM_INT__p_s_q_s);
673 
674 
675                     ostei_general_vrr_I(2, 0, 12, 0, 4,
676                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
677                             PRIM_INT__p_s_q_s, PRIM_INT__s_s_q_s, NULL, PRIM_INT__p_s_o_s, NULL, PRIM_INT__d_s_q_s);
678 
679 
680                     ostei_general_vrr_I(3, 0, 12, 0, 3,
681                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
682                             PRIM_INT__d_s_q_s, PRIM_INT__p_s_q_s, NULL, PRIM_INT__d_s_o_s, NULL, PRIM_INT__f_s_q_s);
683 
684 
685                     ostei_general_vrr_I(4, 0, 12, 0, 2,
686                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
687                             PRIM_INT__f_s_q_s, PRIM_INT__d_s_q_s, NULL, PRIM_INT__f_s_o_s, NULL, PRIM_INT__g_s_q_s);
688 
689 
690                     ostei_general_vrr_I(5, 0, 12, 0, 1,
691                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
692                             PRIM_INT__g_s_q_s, PRIM_INT__f_s_q_s, NULL, PRIM_INT__g_s_o_s, NULL, PRIM_INT__h_s_q_s);
693 
694 
695                     ostei_general_vrr1_K(13, 6,
696                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
697                             PRIM_INT__s_s_q_s, PRIM_INT__s_s_o_s, PRIM_INT__s_s_r_s);
698 
699 
700                     ostei_general_vrr_I(1, 0, 13, 0, 5,
701                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
702                             PRIM_INT__s_s_r_s, NULL, NULL, PRIM_INT__s_s_q_s, NULL, PRIM_INT__p_s_r_s);
703 
704 
705                     ostei_general_vrr_I(2, 0, 13, 0, 4,
706                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
707                             PRIM_INT__p_s_r_s, PRIM_INT__s_s_r_s, NULL, PRIM_INT__p_s_q_s, NULL, PRIM_INT__d_s_r_s);
708 
709 
710                     ostei_general_vrr_I(3, 0, 13, 0, 3,
711                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
712                             PRIM_INT__d_s_r_s, PRIM_INT__p_s_r_s, NULL, PRIM_INT__d_s_q_s, NULL, PRIM_INT__f_s_r_s);
713 
714 
715                     ostei_general_vrr_I(4, 0, 13, 0, 2,
716                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
717                             PRIM_INT__f_s_r_s, PRIM_INT__d_s_r_s, NULL, PRIM_INT__f_s_q_s, NULL, PRIM_INT__g_s_r_s);
718 
719 
720                     ostei_general_vrr_I(5, 0, 13, 0, 1,
721                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
722                             PRIM_INT__g_s_r_s, PRIM_INT__f_s_r_s, NULL, PRIM_INT__g_s_q_s, NULL, PRIM_INT__h_s_r_s);
723 
724 
725 
726 
727                     ////////////////////////////////////
728                     // Accumulate contracted integrals
729                     ////////////////////////////////////
730                     if(lastoffset == 0)
731                     {
732                         contract_all(756, PRIM_INT__h_s_k_s, PRIM_PTR_INT__h_s_k_s);
733                         contract_all(945, PRIM_INT__h_s_l_s, PRIM_PTR_INT__h_s_l_s);
734                         contract_all(1155, PRIM_INT__h_s_m_s, PRIM_PTR_INT__h_s_m_s);
735                         contract_all(1386, PRIM_INT__h_s_n_s, PRIM_PTR_INT__h_s_n_s);
736                         contract_all(1638, PRIM_INT__h_s_o_s, PRIM_PTR_INT__h_s_o_s);
737                         contract_all(1911, PRIM_INT__h_s_q_s, PRIM_PTR_INT__h_s_q_s);
738                         contract_all(2205, PRIM_INT__h_s_r_s, PRIM_PTR_INT__h_s_r_s);
739                     }
740                     else
741                     {
742                         contract(756, shelloffsets, PRIM_INT__h_s_k_s, PRIM_PTR_INT__h_s_k_s);
743                         contract(945, shelloffsets, PRIM_INT__h_s_l_s, PRIM_PTR_INT__h_s_l_s);
744                         contract(1155, shelloffsets, PRIM_INT__h_s_m_s, PRIM_PTR_INT__h_s_m_s);
745                         contract(1386, shelloffsets, PRIM_INT__h_s_n_s, PRIM_PTR_INT__h_s_n_s);
746                         contract(1638, shelloffsets, PRIM_INT__h_s_o_s, PRIM_PTR_INT__h_s_o_s);
747                         contract(1911, shelloffsets, PRIM_INT__h_s_q_s, PRIM_PTR_INT__h_s_q_s);
748                         contract(2205, shelloffsets, PRIM_INT__h_s_r_s, PRIM_PTR_INT__h_s_r_s);
749                         PRIM_PTR_INT__h_s_k_s += lastoffset*756;
750                         PRIM_PTR_INT__h_s_l_s += lastoffset*945;
751                         PRIM_PTR_INT__h_s_m_s += lastoffset*1155;
752                         PRIM_PTR_INT__h_s_n_s += lastoffset*1386;
753                         PRIM_PTR_INT__h_s_o_s += lastoffset*1638;
754                         PRIM_PTR_INT__h_s_q_s += lastoffset*1911;
755                         PRIM_PTR_INT__h_s_r_s += lastoffset*2205;
756                     }
757 
758                 }  // close loop over j
759             }  // close loop over i
760 
761             //Advance to the next batch
762             jstart = SIMINT_SIMD_ROUND(jend);
763 
764             //////////////////////////////////////////////
765             // Contracted integrals: Horizontal recurrance
766             //////////////////////////////////////////////
767 
768 
769 
770 
771             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
772             {
773                 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
774 
775                 // set up HRR pointers
776                 double const * restrict HRR_INT__h_s_k_s = INT__h_s_k_s + abcd * 756;
777                 double const * restrict HRR_INT__h_s_l_s = INT__h_s_l_s + abcd * 945;
778                 double const * restrict HRR_INT__h_s_m_s = INT__h_s_m_s + abcd * 1155;
779                 double const * restrict HRR_INT__h_s_n_s = INT__h_s_n_s + abcd * 1386;
780                 double const * restrict HRR_INT__h_s_o_s = INT__h_s_o_s + abcd * 1638;
781                 double const * restrict HRR_INT__h_s_q_s = INT__h_s_q_s + abcd * 1911;
782                 double const * restrict HRR_INT__h_s_r_s = INT__h_s_r_s + abcd * 2205;
783                 double * restrict HRR_INT__h_s_k_i = INT__h_s_k_i + real_abcd * 21168;
784 
785                 // form INT__h_s_k_p
786                 ostei_general_hrr_L(5, 0, 7, 1, hCD, HRR_INT__h_s_l_s, HRR_INT__h_s_k_s, HRR_INT__h_s_k_p);
787 
788                 // form INT__h_s_l_p
789                 ostei_general_hrr_L(5, 0, 8, 1, hCD, HRR_INT__h_s_m_s, HRR_INT__h_s_l_s, HRR_INT__h_s_l_p);
790 
791                 // form INT__h_s_m_p
792                 ostei_general_hrr_L(5, 0, 9, 1, hCD, HRR_INT__h_s_n_s, HRR_INT__h_s_m_s, HRR_INT__h_s_m_p);
793 
794                 // form INT__h_s_n_p
795                 ostei_general_hrr_L(5, 0, 10, 1, hCD, HRR_INT__h_s_o_s, HRR_INT__h_s_n_s, HRR_INT__h_s_n_p);
796 
797                 // form INT__h_s_o_p
798                 ostei_general_hrr_L(5, 0, 11, 1, hCD, HRR_INT__h_s_q_s, HRR_INT__h_s_o_s, HRR_INT__h_s_o_p);
799 
800                 // form INT__h_s_q_p
801                 ostei_general_hrr_L(5, 0, 12, 1, hCD, HRR_INT__h_s_r_s, HRR_INT__h_s_q_s, HRR_INT__h_s_q_p);
802 
803                 // form INT__h_s_k_d
804                 ostei_general_hrr_L(5, 0, 7, 2, hCD, HRR_INT__h_s_l_p, HRR_INT__h_s_k_p, HRR_INT__h_s_k_d);
805 
806                 // form INT__h_s_l_d
807                 ostei_general_hrr_L(5, 0, 8, 2, hCD, HRR_INT__h_s_m_p, HRR_INT__h_s_l_p, HRR_INT__h_s_l_d);
808 
809                 // form INT__h_s_m_d
810                 ostei_general_hrr_L(5, 0, 9, 2, hCD, HRR_INT__h_s_n_p, HRR_INT__h_s_m_p, HRR_INT__h_s_m_d);
811 
812                 // form INT__h_s_n_d
813                 ostei_general_hrr_L(5, 0, 10, 2, hCD, HRR_INT__h_s_o_p, HRR_INT__h_s_n_p, HRR_INT__h_s_n_d);
814 
815                 // form INT__h_s_o_d
816                 ostei_general_hrr_L(5, 0, 11, 2, hCD, HRR_INT__h_s_q_p, HRR_INT__h_s_o_p, HRR_INT__h_s_o_d);
817 
818                 // form INT__h_s_k_f
819                 ostei_general_hrr_L(5, 0, 7, 3, hCD, HRR_INT__h_s_l_d, HRR_INT__h_s_k_d, HRR_INT__h_s_k_f);
820 
821                 // form INT__h_s_l_f
822                 ostei_general_hrr_L(5, 0, 8, 3, hCD, HRR_INT__h_s_m_d, HRR_INT__h_s_l_d, HRR_INT__h_s_l_f);
823 
824                 // form INT__h_s_m_f
825                 ostei_general_hrr_L(5, 0, 9, 3, hCD, HRR_INT__h_s_n_d, HRR_INT__h_s_m_d, HRR_INT__h_s_m_f);
826 
827                 // form INT__h_s_n_f
828                 ostei_general_hrr_L(5, 0, 10, 3, hCD, HRR_INT__h_s_o_d, HRR_INT__h_s_n_d, HRR_INT__h_s_n_f);
829 
830                 // form INT__h_s_k_g
831                 ostei_general_hrr_L(5, 0, 7, 4, hCD, HRR_INT__h_s_l_f, HRR_INT__h_s_k_f, HRR_INT__h_s_k_g);
832 
833                 // form INT__h_s_l_g
834                 ostei_general_hrr_L(5, 0, 8, 4, hCD, HRR_INT__h_s_m_f, HRR_INT__h_s_l_f, HRR_INT__h_s_l_g);
835 
836                 // form INT__h_s_m_g
837                 ostei_general_hrr_L(5, 0, 9, 4, hCD, HRR_INT__h_s_n_f, HRR_INT__h_s_m_f, HRR_INT__h_s_m_g);
838 
839                 // form INT__h_s_k_h
840                 ostei_general_hrr_L(5, 0, 7, 5, hCD, HRR_INT__h_s_l_g, HRR_INT__h_s_k_g, HRR_INT__h_s_k_h);
841 
842                 // form INT__h_s_l_h
843                 ostei_general_hrr_L(5, 0, 8, 5, hCD, HRR_INT__h_s_m_g, HRR_INT__h_s_l_g, HRR_INT__h_s_l_h);
844 
845                 // form INT__h_s_k_i
846                 ostei_general_hrr_L(5, 0, 7, 6, hCD, HRR_INT__h_s_l_h, HRR_INT__h_s_k_h, HRR_INT__h_s_k_i);
847 
848 
849             }  // close HRR loop
850 
851 
852         }   // close loop cdbatch
853 
854         istart = iend;
855     }  // close loop over ab
856 
857     return P.nshell12_clip * Q.nshell12_clip;
858 }
859 
ostei_s_h_k_i(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__s_h_k_i)860 int ostei_s_h_k_i(struct simint_multi_shellpair const P,
861                   struct simint_multi_shellpair const Q,
862                   double screen_tol,
863                   double * const restrict work,
864                   double * const restrict INT__s_h_k_i)
865 {
866     double P_AB[3*P.nshell12];
867     struct simint_multi_shellpair P_tmp = P;
868     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
869     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
870     P_tmp.AB_x = P_AB;
871     P_tmp.AB_y = P_AB + P.nshell12;
872     P_tmp.AB_z = P_AB + 2*P.nshell12;
873 
874     for(int i = 0; i < P.nshell12; i++)
875     {
876         P_tmp.AB_x[i] = -P.AB_x[i];
877         P_tmp.AB_y[i] = -P.AB_y[i];
878         P_tmp.AB_z[i] = -P.AB_z[i];
879     }
880 
881     int ret = ostei_h_s_k_i(P_tmp, Q, screen_tol, work, INT__s_h_k_i);
882     double buffer[21168] SIMINT_ALIGN_ARRAY_DBL;
883 
884     for(int q = 0; q < ret; q++)
885     {
886         int idx = 0;
887         for(int a = 0; a < 1; ++a)
888         for(int b = 0; b < 21; ++b)
889         for(int c = 0; c < 36; ++c)
890         for(int d = 0; d < 28; ++d)
891             buffer[idx++] = INT__s_h_k_i[q*21168+b*1008+a*1008+c*28+d];
892 
893         memcpy(INT__s_h_k_i+q*21168, buffer, 21168*sizeof(double));
894     }
895 
896     return ret;
897 }
898 
ostei_h_s_i_k(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__h_s_i_k)899 int ostei_h_s_i_k(struct simint_multi_shellpair const P,
900                   struct simint_multi_shellpair const Q,
901                   double screen_tol,
902                   double * const restrict work,
903                   double * const restrict INT__h_s_i_k)
904 {
905     double Q_AB[3*Q.nshell12];
906     struct simint_multi_shellpair Q_tmp = Q;
907     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
908     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
909     Q_tmp.AB_x = Q_AB;
910     Q_tmp.AB_y = Q_AB + Q.nshell12;
911     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
912 
913     for(int i = 0; i < Q.nshell12; i++)
914     {
915         Q_tmp.AB_x[i] = -Q.AB_x[i];
916         Q_tmp.AB_y[i] = -Q.AB_y[i];
917         Q_tmp.AB_z[i] = -Q.AB_z[i];
918     }
919 
920     int ret = ostei_h_s_k_i(P, Q_tmp, screen_tol, work, INT__h_s_i_k);
921     double buffer[21168] SIMINT_ALIGN_ARRAY_DBL;
922 
923     for(int q = 0; q < ret; q++)
924     {
925         int idx = 0;
926         for(int a = 0; a < 21; ++a)
927         for(int b = 0; b < 1; ++b)
928         for(int c = 0; c < 28; ++c)
929         for(int d = 0; d < 36; ++d)
930             buffer[idx++] = INT__h_s_i_k[q*21168+a*1008+b*1008+d*28+c];
931 
932         memcpy(INT__h_s_i_k+q*21168, buffer, 21168*sizeof(double));
933     }
934 
935     return ret;
936 }
937 
ostei_s_h_i_k(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__s_h_i_k)938 int ostei_s_h_i_k(struct simint_multi_shellpair const P,
939                   struct simint_multi_shellpair const Q,
940                   double screen_tol,
941                   double * const restrict work,
942                   double * const restrict INT__s_h_i_k)
943 {
944     double P_AB[3*P.nshell12];
945     struct simint_multi_shellpair P_tmp = P;
946     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
947     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
948     P_tmp.AB_x = P_AB;
949     P_tmp.AB_y = P_AB + P.nshell12;
950     P_tmp.AB_z = P_AB + 2*P.nshell12;
951 
952     for(int i = 0; i < P.nshell12; i++)
953     {
954         P_tmp.AB_x[i] = -P.AB_x[i];
955         P_tmp.AB_y[i] = -P.AB_y[i];
956         P_tmp.AB_z[i] = -P.AB_z[i];
957     }
958 
959     double Q_AB[3*Q.nshell12];
960     struct simint_multi_shellpair Q_tmp = Q;
961     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
962     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
963     Q_tmp.AB_x = Q_AB;
964     Q_tmp.AB_y = Q_AB + Q.nshell12;
965     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
966 
967     for(int i = 0; i < Q.nshell12; i++)
968     {
969         Q_tmp.AB_x[i] = -Q.AB_x[i];
970         Q_tmp.AB_y[i] = -Q.AB_y[i];
971         Q_tmp.AB_z[i] = -Q.AB_z[i];
972     }
973 
974     int ret = ostei_h_s_k_i(P_tmp, Q_tmp, screen_tol, work, INT__s_h_i_k);
975     double buffer[21168] SIMINT_ALIGN_ARRAY_DBL;
976 
977     for(int q = 0; q < ret; q++)
978     {
979         int idx = 0;
980         for(int a = 0; a < 1; ++a)
981         for(int b = 0; b < 21; ++b)
982         for(int c = 0; c < 28; ++c)
983         for(int d = 0; d < 36; ++d)
984             buffer[idx++] = INT__s_h_i_k[q*21168+b*1008+a*1008+d*28+c];
985 
986         memcpy(INT__s_h_i_k+q*21168, buffer, 21168*sizeof(double));
987     }
988 
989     return ret;
990 }
991 
992