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_g_f_k_f(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__g_f_k_f)8 int ostei_g_f_k_f(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__g_f_k_f)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__g_f_k_f);
17     int ab, cd, abcd;
18     int istart, jstart;
19     int iprimcd, nprim_icd, icd;
20     const int check_screen = (screen_tol > 0.0);
21     int i, j;
22     int n;
23     int not_screened;
24     int real_abcd;
25     int iket;
26     int ibra;
27 
28     // partition workspace
29     double * const INT__g_s_k_s = work + (SIMINT_NSHELL_SIMD * 0);
30     double * const INT__g_s_l_s = work + (SIMINT_NSHELL_SIMD * 540);
31     double * const INT__g_s_m_s = work + (SIMINT_NSHELL_SIMD * 1215);
32     double * const INT__g_s_n_s = work + (SIMINT_NSHELL_SIMD * 2040);
33     double * const INT__h_s_k_s = work + (SIMINT_NSHELL_SIMD * 3030);
34     double * const INT__h_s_l_s = work + (SIMINT_NSHELL_SIMD * 3786);
35     double * const INT__h_s_m_s = work + (SIMINT_NSHELL_SIMD * 4731);
36     double * const INT__h_s_n_s = work + (SIMINT_NSHELL_SIMD * 5886);
37     double * const INT__i_s_k_s = work + (SIMINT_NSHELL_SIMD * 7272);
38     double * const INT__i_s_l_s = work + (SIMINT_NSHELL_SIMD * 8280);
39     double * const INT__i_s_m_s = work + (SIMINT_NSHELL_SIMD * 9540);
40     double * const INT__i_s_n_s = work + (SIMINT_NSHELL_SIMD * 11080);
41     double * const INT__k_s_k_s = work + (SIMINT_NSHELL_SIMD * 12928);
42     double * const INT__k_s_l_s = work + (SIMINT_NSHELL_SIMD * 14224);
43     double * const INT__k_s_m_s = work + (SIMINT_NSHELL_SIMD * 15844);
44     double * const INT__k_s_n_s = work + (SIMINT_NSHELL_SIMD * 17824);
45     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*20200);
46     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_p_s = primwork + 18;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_d_s = primwork + 69;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_f_s = primwork + 165;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_g_s = primwork + 315;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_h_s = primwork + 525;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_i_s = primwork + 798;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_k_s = primwork + 1134;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_l_s = primwork + 1530;
55     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_m_s = primwork + 1980;
56     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_n_s = primwork + 2475;
57     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_p_s = primwork + 3003;
58     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_d_s = primwork + 3066;
59     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_f_s = primwork + 3192;
60     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_g_s = primwork + 3402;
61     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_h_s = primwork + 3717;
62     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_i_s = primwork + 4158;
63     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_k_s = primwork + 4746;
64     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_l_s = primwork + 5502;
65     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_m_s = primwork + 6447;
66     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_n_s = primwork + 7602;
67     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_d_s = primwork + 8988;
68     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_f_s = primwork + 9204;
69     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_g_s = primwork + 9564;
70     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_h_s = primwork + 10104;
71     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_i_s = primwork + 10860;
72     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_k_s = primwork + 11868;
73     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_l_s = primwork + 13164;
74     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_m_s = primwork + 14784;
75     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_n_s = primwork + 16764;
76     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_f_s = primwork + 19140;
77     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_g_s = primwork + 19640;
78     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_h_s = primwork + 20390;
79     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_i_s = primwork + 21440;
80     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_k_s = primwork + 22840;
81     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_l_s = primwork + 24640;
82     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_m_s = primwork + 26890;
83     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_n_s = primwork + 29640;
84     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_g_s = primwork + 32940;
85     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_h_s = primwork + 33840;
86     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_i_s = primwork + 35100;
87     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_k_s = primwork + 36780;
88     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_l_s = primwork + 38940;
89     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_m_s = primwork + 41640;
90     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_n_s = primwork + 44940;
91     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_h_s = primwork + 48900;
92     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_i_s = primwork + 50223;
93     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_k_s = primwork + 51987;
94     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_l_s = primwork + 54255;
95     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_m_s = primwork + 57090;
96     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_n_s = primwork + 60555;
97     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_i_s = primwork + 64713;
98     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_k_s = primwork + 66281;
99     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_l_s = primwork + 68297;
100     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_m_s = primwork + 70817;
101     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_n_s = primwork + 73897;
102     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_k_s = primwork + 77593;
103     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_l_s = primwork + 78889;
104     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_m_s = primwork + 80509;
105     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_n_s = primwork + 82489;
106     double * const hrrwork = (double *)(primwork + 84865);
107     double * const HRR_INT__g_p_k_s = hrrwork + 0;
108     double * const HRR_INT__g_p_l_s = hrrwork + 1620;
109     double * const HRR_INT__g_p_m_s = hrrwork + 3645;
110     double * const HRR_INT__g_p_n_s = hrrwork + 6120;
111     double * const HRR_INT__g_d_k_s = hrrwork + 9090;
112     double * const HRR_INT__g_d_l_s = hrrwork + 12330;
113     double * const HRR_INT__g_d_m_s = hrrwork + 16380;
114     double * const HRR_INT__g_d_n_s = hrrwork + 21330;
115     double * const HRR_INT__g_f_k_s = hrrwork + 27270;
116     double * const HRR_INT__g_f_k_p = hrrwork + 32670;
117     double * const HRR_INT__g_f_k_d = hrrwork + 48870;
118     double * const HRR_INT__g_f_l_s = hrrwork + 81270;
119     double * const HRR_INT__g_f_l_p = hrrwork + 88020;
120     double * const HRR_INT__g_f_l_d = hrrwork + 108270;
121     double * const HRR_INT__g_f_m_s = hrrwork + 148770;
122     double * const HRR_INT__g_f_m_p = hrrwork + 157020;
123     double * const HRR_INT__g_f_n_s = hrrwork + 181770;
124     double * const HRR_INT__h_p_k_s = hrrwork + 191670;
125     double * const HRR_INT__h_p_l_s = hrrwork + 193938;
126     double * const HRR_INT__h_p_m_s = hrrwork + 196773;
127     double * const HRR_INT__h_p_n_s = hrrwork + 200238;
128     double * const HRR_INT__h_d_k_s = hrrwork + 204396;
129     double * const HRR_INT__h_d_l_s = hrrwork + 208932;
130     double * const HRR_INT__h_d_m_s = hrrwork + 214602;
131     double * const HRR_INT__h_d_n_s = hrrwork + 221532;
132     double * const HRR_INT__i_p_k_s = hrrwork + 229848;
133     double * const HRR_INT__i_p_l_s = hrrwork + 232872;
134     double * const HRR_INT__i_p_m_s = hrrwork + 236652;
135     double * const HRR_INT__i_p_n_s = hrrwork + 241272;
136 
137 
138     // Create constants
139     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
140     const SIMINT_DBLTYPE const_10 = SIMINT_DBLSET1(10);
141     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
142     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
143     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
144     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
145     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
146     const SIMINT_DBLTYPE const_7 = SIMINT_DBLSET1(7);
147     const SIMINT_DBLTYPE const_8 = SIMINT_DBLSET1(8);
148     const SIMINT_DBLTYPE const_9 = SIMINT_DBLSET1(9);
149     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
150 
151 
152     ////////////////////////////////////////
153     // Loop over shells and primitives
154     ////////////////////////////////////////
155 
156     real_abcd = 0;
157     istart = 0;
158     for(ab = 0; ab < P.nshell12_clip; ++ab)
159     {
160         const int iend = istart + P.nprim12[ab];
161 
162         cd = 0;
163         jstart = 0;
164 
165         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
166         {
167             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
168             int jend = jstart;
169             for(i = 0; i < nshellbatch; i++)
170                 jend += Q.nprim12[cd+i];
171 
172             // Clear the beginning of the workspace (where we are accumulating integrals)
173             memset(work, 0, SIMINT_NSHELL_SIMD * 20200 * sizeof(double));
174             abcd = 0;
175 
176 
177             for(i = istart; i < iend; ++i)
178             {
179                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
180 
181                 if(check_screen)
182                 {
183                     // Skip this whole thing if always insignificant
184                     if((P.screen[i] * Q.screen_max) < screen_tol)
185                         continue;
186                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
187                 }
188 
189                 icd = 0;
190                 iprimcd = 0;
191                 nprim_icd = Q.nprim12[cd];
192                 double * restrict PRIM_PTR_INT__g_s_k_s = INT__g_s_k_s + abcd * 540;
193                 double * restrict PRIM_PTR_INT__g_s_l_s = INT__g_s_l_s + abcd * 675;
194                 double * restrict PRIM_PTR_INT__g_s_m_s = INT__g_s_m_s + abcd * 825;
195                 double * restrict PRIM_PTR_INT__g_s_n_s = INT__g_s_n_s + abcd * 990;
196                 double * restrict PRIM_PTR_INT__h_s_k_s = INT__h_s_k_s + abcd * 756;
197                 double * restrict PRIM_PTR_INT__h_s_l_s = INT__h_s_l_s + abcd * 945;
198                 double * restrict PRIM_PTR_INT__h_s_m_s = INT__h_s_m_s + abcd * 1155;
199                 double * restrict PRIM_PTR_INT__h_s_n_s = INT__h_s_n_s + abcd * 1386;
200                 double * restrict PRIM_PTR_INT__i_s_k_s = INT__i_s_k_s + abcd * 1008;
201                 double * restrict PRIM_PTR_INT__i_s_l_s = INT__i_s_l_s + abcd * 1260;
202                 double * restrict PRIM_PTR_INT__i_s_m_s = INT__i_s_m_s + abcd * 1540;
203                 double * restrict PRIM_PTR_INT__i_s_n_s = INT__i_s_n_s + abcd * 1848;
204                 double * restrict PRIM_PTR_INT__k_s_k_s = INT__k_s_k_s + abcd * 1296;
205                 double * restrict PRIM_PTR_INT__k_s_l_s = INT__k_s_l_s + abcd * 1620;
206                 double * restrict PRIM_PTR_INT__k_s_m_s = INT__k_s_m_s + abcd * 1980;
207                 double * restrict PRIM_PTR_INT__k_s_n_s = INT__k_s_n_s + abcd * 2376;
208 
209 
210 
211                 // Load these one per loop over i
212                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
213                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
214                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
215 
216                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
217 
218                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
219                 {
220                     // calculate the shell offsets
221                     // these are the offset from the shell pointed to by cd
222                     // for each element
223                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
224                     int lastoffset = 0;
225                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
226 
227                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
228                     {
229                         // Handle if the first element of the vector is a new shell
230                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
231                         {
232                             nprim_icd += Q.nprim12[cd + (++icd)];
233                             PRIM_PTR_INT__g_s_k_s += 540;
234                             PRIM_PTR_INT__g_s_l_s += 675;
235                             PRIM_PTR_INT__g_s_m_s += 825;
236                             PRIM_PTR_INT__g_s_n_s += 990;
237                             PRIM_PTR_INT__h_s_k_s += 756;
238                             PRIM_PTR_INT__h_s_l_s += 945;
239                             PRIM_PTR_INT__h_s_m_s += 1155;
240                             PRIM_PTR_INT__h_s_n_s += 1386;
241                             PRIM_PTR_INT__i_s_k_s += 1008;
242                             PRIM_PTR_INT__i_s_l_s += 1260;
243                             PRIM_PTR_INT__i_s_m_s += 1540;
244                             PRIM_PTR_INT__i_s_n_s += 1848;
245                             PRIM_PTR_INT__k_s_k_s += 1296;
246                             PRIM_PTR_INT__k_s_l_s += 1620;
247                             PRIM_PTR_INT__k_s_m_s += 1980;
248                             PRIM_PTR_INT__k_s_n_s += 2376;
249                         }
250                         iprimcd++;
251                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
252                         {
253                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
254                             {
255                                 shelloffsets[n] = shelloffsets[n-1] + 1;
256                                 lastoffset++;
257                                 nprim_icd += Q.nprim12[cd + (++icd)];
258                             }
259                             else
260                                 shelloffsets[n] = shelloffsets[n-1];
261                             iprimcd++;
262                         }
263                     }
264                     else
265                         iprimcd += SIMINT_SIMD_LEN;
266 
267                     // Do we have to compute this vector (or has it been screened out)?
268                     // (not_screened != 0 means we have to do this vector)
269                     if(check_screen)
270                     {
271                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
272                         if(vmax < screen_tol)
273                         {
274                             PRIM_PTR_INT__g_s_k_s += lastoffset*540;
275                             PRIM_PTR_INT__g_s_l_s += lastoffset*675;
276                             PRIM_PTR_INT__g_s_m_s += lastoffset*825;
277                             PRIM_PTR_INT__g_s_n_s += lastoffset*990;
278                             PRIM_PTR_INT__h_s_k_s += lastoffset*756;
279                             PRIM_PTR_INT__h_s_l_s += lastoffset*945;
280                             PRIM_PTR_INT__h_s_m_s += lastoffset*1155;
281                             PRIM_PTR_INT__h_s_n_s += lastoffset*1386;
282                             PRIM_PTR_INT__i_s_k_s += lastoffset*1008;
283                             PRIM_PTR_INT__i_s_l_s += lastoffset*1260;
284                             PRIM_PTR_INT__i_s_m_s += lastoffset*1540;
285                             PRIM_PTR_INT__i_s_n_s += lastoffset*1848;
286                             PRIM_PTR_INT__k_s_k_s += lastoffset*1296;
287                             PRIM_PTR_INT__k_s_l_s += lastoffset*1620;
288                             PRIM_PTR_INT__k_s_m_s += lastoffset*1980;
289                             PRIM_PTR_INT__k_s_n_s += lastoffset*2376;
290                             continue;
291                         }
292                     }
293 
294                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
295                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
296                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
297                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
298 
299 
300                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
301                     SIMINT_DBLTYPE PQ[3];
302                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
303                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
304                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
305                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
306                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
307                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
308 
309                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
310                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
311                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
312                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
313                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
314                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
315                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
316 
317                     // NOTE: Minus sign!
318                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
319                     SIMINT_DBLTYPE aop_PQ[3];
320                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
321                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
322                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
323 
324                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
325                     SIMINT_DBLTYPE aoq_PQ[3];
326                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
327                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
328                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
329                     // Put a minus sign here so we don't have to in RR routines
330                     a_over_q = SIMINT_NEG(a_over_q);
331 
332 
333                     //////////////////////////////////////////////
334                     // Fjt function section
335                     // Maximum v value: 17
336                     //////////////////////////////////////////////
337                     // The parameter to the Fjt function
338                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
339 
340 
341                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
342 
343 
344                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 17);
345                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
346                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
347                     for(n = 0; n <= 17; n++)
348                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
349 
350                     //////////////////////////////////////////////
351                     // Primitive integrals: Vertical recurrance
352                     //////////////////////////////////////////////
353 
354                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
355                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
356                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
357                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
358                     const SIMINT_DBLTYPE vrr_const_5_over_2p = SIMINT_MUL(const_5, one_over_2p);
359                     const SIMINT_DBLTYPE vrr_const_6_over_2p = SIMINT_MUL(const_6, one_over_2p);
360                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
361                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
362                     const SIMINT_DBLTYPE vrr_const_3_over_2q = SIMINT_MUL(const_3, one_over_2q);
363                     const SIMINT_DBLTYPE vrr_const_4_over_2q = SIMINT_MUL(const_4, one_over_2q);
364                     const SIMINT_DBLTYPE vrr_const_5_over_2q = SIMINT_MUL(const_5, one_over_2q);
365                     const SIMINT_DBLTYPE vrr_const_6_over_2q = SIMINT_MUL(const_6, one_over_2q);
366                     const SIMINT_DBLTYPE vrr_const_7_over_2q = SIMINT_MUL(const_7, one_over_2q);
367                     const SIMINT_DBLTYPE vrr_const_8_over_2q = SIMINT_MUL(const_8, one_over_2q);
368                     const SIMINT_DBLTYPE vrr_const_9_over_2q = SIMINT_MUL(const_9, one_over_2q);
369                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
370                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
371                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
372                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
373                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
374                     const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
375                     const SIMINT_DBLTYPE vrr_const_7_over_2pq = SIMINT_MUL(const_7, one_over_2pq);
376                     const SIMINT_DBLTYPE vrr_const_8_over_2pq = SIMINT_MUL(const_8, one_over_2pq);
377                     const SIMINT_DBLTYPE vrr_const_9_over_2pq = SIMINT_MUL(const_9, one_over_2pq);
378                     const SIMINT_DBLTYPE vrr_const_10_over_2pq = SIMINT_MUL(const_10, one_over_2pq);
379 
380 
381 
382                     // Forming PRIM_INT__s_s_p_s[17 * 3];
383                     for(n = 0; n < 17; ++n)  // loop over orders of auxiliary function
384                     {
385 
386                         PRIM_INT__s_s_p_s[n * 3 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
387                         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]);
388 
389                         PRIM_INT__s_s_p_s[n * 3 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
390                         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]);
391 
392                         PRIM_INT__s_s_p_s[n * 3 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
393                         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]);
394 
395                     }
396 
397 
398 
399                     // Forming PRIM_INT__s_s_d_s[16 * 6];
400                     for(n = 0; n < 16; ++n)  // loop over orders of auxiliary function
401                     {
402 
403                         PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_p_s[n * 3 + 0]);
404                         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]);
405                         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]);
406 
407                         PRIM_INT__s_s_d_s[n * 6 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_p_s[n * 3 + 0]);
408                         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]);
409 
410                         PRIM_INT__s_s_d_s[n * 6 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 0]);
411                         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]);
412 
413                         PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_p_s[n * 3 + 1]);
414                         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]);
415                         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]);
416 
417                         PRIM_INT__s_s_d_s[n * 6 + 4] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 1]);
418                         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]);
419 
420                         PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 2]);
421                         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]);
422                         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]);
423 
424                     }
425 
426 
427 
428                     // Forming PRIM_INT__s_s_f_s[15 * 10];
429                     for(n = 0; n < 15; ++n)  // loop over orders of auxiliary function
430                     {
431 
432                         PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 0]);
433                         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]);
434                         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]);
435 
436                         PRIM_INT__s_s_f_s[n * 10 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 0]);
437                         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]);
438 
439                         PRIM_INT__s_s_f_s[n * 10 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 0]);
440                         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]);
441 
442                         PRIM_INT__s_s_f_s[n * 10 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 3]);
443                         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]);
444 
445                         PRIM_INT__s_s_f_s[n * 10 + 4] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 1]);
446                         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]);
447 
448                         PRIM_INT__s_s_f_s[n * 10 + 5] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 5]);
449                         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]);
450 
451                         PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 3]);
452                         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]);
453                         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]);
454 
455                         PRIM_INT__s_s_f_s[n * 10 + 7] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 3]);
456                         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]);
457 
458                         PRIM_INT__s_s_f_s[n * 10 + 8] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 5]);
459                         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]);
460 
461                         PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 5]);
462                         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]);
463                         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]);
464 
465                     }
466 
467 
468                     VRR_K_s_s_g_s(
469                             PRIM_INT__s_s_g_s,
470                             PRIM_INT__s_s_f_s,
471                             PRIM_INT__s_s_d_s,
472                             Q_PA,
473                             a_over_q,
474                             aoq_PQ,
475                             one_over_2q,
476                             14);
477 
478 
479                     VRR_K_s_s_h_s(
480                             PRIM_INT__s_s_h_s,
481                             PRIM_INT__s_s_g_s,
482                             PRIM_INT__s_s_f_s,
483                             Q_PA,
484                             a_over_q,
485                             aoq_PQ,
486                             one_over_2q,
487                             13);
488 
489 
490                     ostei_general_vrr1_K(6, 12,
491                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
492                             PRIM_INT__s_s_h_s, PRIM_INT__s_s_g_s, PRIM_INT__s_s_i_s);
493 
494 
495                     ostei_general_vrr1_K(7, 11,
496                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
497                             PRIM_INT__s_s_i_s, PRIM_INT__s_s_h_s, PRIM_INT__s_s_k_s);
498 
499 
500                     ostei_general_vrr_I(1, 0, 7, 0, 7,
501                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
502                             PRIM_INT__s_s_k_s, NULL, NULL, PRIM_INT__s_s_i_s, NULL, PRIM_INT__p_s_k_s);
503 
504 
505                     ostei_general_vrr_I(1, 0, 6, 0, 7,
506                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
507                             PRIM_INT__s_s_i_s, NULL, NULL, PRIM_INT__s_s_h_s, NULL, PRIM_INT__p_s_i_s);
508 
509 
510                     ostei_general_vrr_I(2, 0, 7, 0, 6,
511                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
512                             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);
513 
514 
515                     ostei_general_vrr_I(1, 0, 5, 0, 7,
516                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
517                             PRIM_INT__s_s_h_s, NULL, NULL, PRIM_INT__s_s_g_s, NULL, PRIM_INT__p_s_h_s);
518 
519 
520                     ostei_general_vrr_I(2, 0, 6, 0, 6,
521                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
522                             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);
523 
524 
525                     ostei_general_vrr_I(3, 0, 7, 0, 5,
526                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
527                             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);
528 
529 
530                     VRR_I_p_s_g_s(
531                             PRIM_INT__p_s_g_s,
532                             PRIM_INT__s_s_g_s,
533                             PRIM_INT__s_s_f_s,
534                             P_PA,
535                             aop_PQ,
536                             one_over_2pq,
537                             7);
538 
539 
540                     ostei_general_vrr_I(2, 0, 5, 0, 6,
541                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
542                             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);
543 
544 
545                     ostei_general_vrr_I(3, 0, 6, 0, 5,
546                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
547                             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);
548 
549 
550                     ostei_general_vrr_I(4, 0, 7, 0, 4,
551                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
552                             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);
553 
554 
555                     ostei_general_vrr1_K(8, 10,
556                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
557                             PRIM_INT__s_s_k_s, PRIM_INT__s_s_i_s, PRIM_INT__s_s_l_s);
558 
559 
560                     ostei_general_vrr_I(1, 0, 8, 0, 7,
561                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
562                             PRIM_INT__s_s_l_s, NULL, NULL, PRIM_INT__s_s_k_s, NULL, PRIM_INT__p_s_l_s);
563 
564 
565                     ostei_general_vrr_I(2, 0, 8, 0, 6,
566                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
567                             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);
568 
569 
570                     ostei_general_vrr_I(3, 0, 8, 0, 5,
571                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
572                             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);
573 
574 
575                     ostei_general_vrr_I(4, 0, 8, 0, 4,
576                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
577                             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);
578 
579 
580                     VRR_I_p_s_f_s(
581                             PRIM_INT__p_s_f_s,
582                             PRIM_INT__s_s_f_s,
583                             PRIM_INT__s_s_d_s,
584                             P_PA,
585                             aop_PQ,
586                             one_over_2pq,
587                             7);
588 
589 
590                     ostei_general_vrr_I(2, 0, 4, 0, 6,
591                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
592                             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);
593 
594 
595                     ostei_general_vrr_I(3, 0, 5, 0, 5,
596                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
597                             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);
598 
599 
600                     ostei_general_vrr_I(4, 0, 6, 0, 4,
601                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
602                             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);
603 
604 
605                     ostei_general_vrr_I(5, 0, 7, 0, 3,
606                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
607                             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);
608 
609 
610                     ostei_general_vrr1_K(9, 9,
611                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
612                             PRIM_INT__s_s_l_s, PRIM_INT__s_s_k_s, PRIM_INT__s_s_m_s);
613 
614 
615                     ostei_general_vrr_I(1, 0, 9, 0, 7,
616                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
617                             PRIM_INT__s_s_m_s, NULL, NULL, PRIM_INT__s_s_l_s, NULL, PRIM_INT__p_s_m_s);
618 
619 
620                     ostei_general_vrr_I(2, 0, 9, 0, 6,
621                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
622                             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);
623 
624 
625                     ostei_general_vrr_I(3, 0, 9, 0, 5,
626                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
627                             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);
628 
629 
630                     ostei_general_vrr_I(4, 0, 9, 0, 4,
631                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
632                             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);
633 
634 
635                     ostei_general_vrr_I(5, 0, 8, 0, 3,
636                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
637                             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);
638 
639 
640 
641                     // Forming PRIM_INT__p_s_d_s[7 * 18];
642                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
643                     {
644 
645                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 0]);
646                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__p_s_d_s[n * 18 + 0]);
647                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_d_s[n * 18 + 0]);
648 
649                         PRIM_INT__p_s_d_s[n * 18 + 1] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 1]);
650                         PRIM_INT__p_s_d_s[n * 18 + 1] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__p_s_d_s[n * 18 + 1]);
651                         PRIM_INT__p_s_d_s[n * 18 + 1] = 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 + 1]);
652 
653                         PRIM_INT__p_s_d_s[n * 18 + 2] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 2]);
654                         PRIM_INT__p_s_d_s[n * 18 + 2] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 2], PRIM_INT__p_s_d_s[n * 18 + 2]);
655                         PRIM_INT__p_s_d_s[n * 18 + 2] = 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 + 2]);
656 
657                         PRIM_INT__p_s_d_s[n * 18 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 3]);
658                         PRIM_INT__p_s_d_s[n * 18 + 3] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__p_s_d_s[n * 18 + 3]);
659 
660                         PRIM_INT__p_s_d_s[n * 18 + 4] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 4]);
661                         PRIM_INT__p_s_d_s[n * 18 + 4] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 4], PRIM_INT__p_s_d_s[n * 18 + 4]);
662 
663                         PRIM_INT__p_s_d_s[n * 18 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 5]);
664                         PRIM_INT__p_s_d_s[n * 18 + 5] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__p_s_d_s[n * 18 + 5]);
665 
666                         PRIM_INT__p_s_d_s[n * 18 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 0]);
667                         PRIM_INT__p_s_d_s[n * 18 + 6] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__p_s_d_s[n * 18 + 6]);
668 
669                         PRIM_INT__p_s_d_s[n * 18 + 7] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 1]);
670                         PRIM_INT__p_s_d_s[n * 18 + 7] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__p_s_d_s[n * 18 + 7]);
671                         PRIM_INT__p_s_d_s[n * 18 + 7] = 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 + 7]);
672 
673                         PRIM_INT__p_s_d_s[n * 18 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 2]);
674                         PRIM_INT__p_s_d_s[n * 18 + 8] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 2], PRIM_INT__p_s_d_s[n * 18 + 8]);
675 
676                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 3]);
677                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__p_s_d_s[n * 18 + 9]);
678                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_d_s[n * 18 + 9]);
679 
680                         PRIM_INT__p_s_d_s[n * 18 + 10] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 4]);
681                         PRIM_INT__p_s_d_s[n * 18 + 10] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 4], PRIM_INT__p_s_d_s[n * 18 + 10]);
682                         PRIM_INT__p_s_d_s[n * 18 + 10] = 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 + 10]);
683 
684                         PRIM_INT__p_s_d_s[n * 18 + 11] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 5]);
685                         PRIM_INT__p_s_d_s[n * 18 + 11] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__p_s_d_s[n * 18 + 11]);
686 
687                         PRIM_INT__p_s_d_s[n * 18 + 12] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 0]);
688                         PRIM_INT__p_s_d_s[n * 18 + 12] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__p_s_d_s[n * 18 + 12]);
689 
690                         PRIM_INT__p_s_d_s[n * 18 + 13] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 1]);
691                         PRIM_INT__p_s_d_s[n * 18 + 13] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__p_s_d_s[n * 18 + 13]);
692 
693                         PRIM_INT__p_s_d_s[n * 18 + 14] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 2]);
694                         PRIM_INT__p_s_d_s[n * 18 + 14] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 2], PRIM_INT__p_s_d_s[n * 18 + 14]);
695                         PRIM_INT__p_s_d_s[n * 18 + 14] = 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 + 14]);
696 
697                         PRIM_INT__p_s_d_s[n * 18 + 15] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 3]);
698                         PRIM_INT__p_s_d_s[n * 18 + 15] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__p_s_d_s[n * 18 + 15]);
699 
700                         PRIM_INT__p_s_d_s[n * 18 + 16] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 4]);
701                         PRIM_INT__p_s_d_s[n * 18 + 16] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 4], PRIM_INT__p_s_d_s[n * 18 + 16]);
702                         PRIM_INT__p_s_d_s[n * 18 + 16] = 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 + 16]);
703 
704                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 5]);
705                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__p_s_d_s[n * 18 + 17]);
706                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_d_s[n * 18 + 17]);
707 
708                     }
709 
710 
711                     VRR_I_d_s_f_s(
712                             PRIM_INT__d_s_f_s,
713                             PRIM_INT__p_s_f_s,
714                             PRIM_INT__s_s_f_s,
715                             PRIM_INT__p_s_d_s,
716                             P_PA,
717                             a_over_p,
718                             aop_PQ,
719                             one_over_2p,
720                             one_over_2pq,
721                             6);
722 
723 
724                     ostei_general_vrr_I(3, 0, 4, 0, 5,
725                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
726                             PRIM_INT__d_s_g_s, PRIM_INT__p_s_g_s, NULL, PRIM_INT__d_s_f_s, NULL, PRIM_INT__f_s_g_s);
727 
728 
729                     ostei_general_vrr_I(4, 0, 5, 0, 4,
730                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
731                             PRIM_INT__f_s_h_s, PRIM_INT__d_s_h_s, NULL, PRIM_INT__f_s_g_s, NULL, PRIM_INT__g_s_h_s);
732 
733 
734                     ostei_general_vrr_I(5, 0, 6, 0, 3,
735                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
736                             PRIM_INT__g_s_i_s, PRIM_INT__f_s_i_s, NULL, PRIM_INT__g_s_h_s, NULL, PRIM_INT__h_s_i_s);
737 
738 
739                     ostei_general_vrr_I(6, 0, 7, 0, 2,
740                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
741                             PRIM_INT__h_s_k_s, PRIM_INT__g_s_k_s, NULL, PRIM_INT__h_s_i_s, NULL, PRIM_INT__i_s_k_s);
742 
743 
744                     ostei_general_vrr1_K(10, 8,
745                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
746                             PRIM_INT__s_s_m_s, PRIM_INT__s_s_l_s, PRIM_INT__s_s_n_s);
747 
748 
749                     ostei_general_vrr_I(1, 0, 10, 0, 7,
750                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
751                             PRIM_INT__s_s_n_s, NULL, NULL, PRIM_INT__s_s_m_s, NULL, PRIM_INT__p_s_n_s);
752 
753 
754                     ostei_general_vrr_I(2, 0, 10, 0, 6,
755                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
756                             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);
757 
758 
759                     ostei_general_vrr_I(3, 0, 10, 0, 5,
760                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
761                             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);
762 
763 
764                     ostei_general_vrr_I(4, 0, 10, 0, 4,
765                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
766                             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);
767 
768 
769                     ostei_general_vrr_I(5, 0, 9, 0, 3,
770                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
771                             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);
772 
773 
774                     ostei_general_vrr_I(6, 0, 8, 0, 2,
775                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
776                             PRIM_INT__h_s_l_s, PRIM_INT__g_s_l_s, NULL, PRIM_INT__h_s_k_s, NULL, PRIM_INT__i_s_l_s);
777 
778 
779 
780                     // Forming PRIM_INT__p_s_p_s[7 * 9];
781                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
782                     {
783 
784                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_p_s[n * 3 + 0]);
785                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 0]);
786                         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]);
787 
788                         PRIM_INT__p_s_p_s[n * 9 + 1] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_p_s[n * 3 + 1]);
789                         PRIM_INT__p_s_p_s[n * 9 + 1] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 1]);
790 
791                         PRIM_INT__p_s_p_s[n * 9 + 2] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_p_s[n * 3 + 2]);
792                         PRIM_INT__p_s_p_s[n * 9 + 2] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 2]);
793 
794                         PRIM_INT__p_s_p_s[n * 9 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_p_s[n * 3 + 0]);
795                         PRIM_INT__p_s_p_s[n * 9 + 3] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 3]);
796 
797                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_p_s[n * 3 + 1]);
798                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 4]);
799                         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]);
800 
801                         PRIM_INT__p_s_p_s[n * 9 + 5] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_p_s[n * 3 + 2]);
802                         PRIM_INT__p_s_p_s[n * 9 + 5] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 5]);
803 
804                         PRIM_INT__p_s_p_s[n * 9 + 6] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_p_s[n * 3 + 0]);
805                         PRIM_INT__p_s_p_s[n * 9 + 6] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 6]);
806 
807                         PRIM_INT__p_s_p_s[n * 9 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_p_s[n * 3 + 1]);
808                         PRIM_INT__p_s_p_s[n * 9 + 7] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 7]);
809 
810                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_p_s[n * 3 + 2]);
811                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 8]);
812                         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]);
813 
814                     }
815 
816 
817                     VRR_I_d_s_d_s(
818                             PRIM_INT__d_s_d_s,
819                             PRIM_INT__p_s_d_s,
820                             PRIM_INT__s_s_d_s,
821                             PRIM_INT__p_s_p_s,
822                             P_PA,
823                             a_over_p,
824                             aop_PQ,
825                             one_over_2p,
826                             one_over_2pq,
827                             6);
828 
829 
830                     ostei_general_vrr_I(3, 0, 3, 0, 5,
831                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
832                             PRIM_INT__d_s_f_s, PRIM_INT__p_s_f_s, NULL, PRIM_INT__d_s_d_s, NULL, PRIM_INT__f_s_f_s);
833 
834 
835                     ostei_general_vrr_I(4, 0, 4, 0, 4,
836                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
837                             PRIM_INT__f_s_g_s, PRIM_INT__d_s_g_s, NULL, PRIM_INT__f_s_f_s, NULL, PRIM_INT__g_s_g_s);
838 
839 
840                     ostei_general_vrr_I(5, 0, 5, 0, 3,
841                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
842                             PRIM_INT__g_s_h_s, PRIM_INT__f_s_h_s, NULL, PRIM_INT__g_s_g_s, NULL, PRIM_INT__h_s_h_s);
843 
844 
845                     ostei_general_vrr_I(6, 0, 6, 0, 2,
846                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
847                             PRIM_INT__h_s_i_s, PRIM_INT__g_s_i_s, NULL, PRIM_INT__h_s_h_s, NULL, PRIM_INT__i_s_i_s);
848 
849 
850                     ostei_general_vrr_I(7, 0, 7, 0, 1,
851                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
852                             PRIM_INT__i_s_k_s, PRIM_INT__h_s_k_s, NULL, PRIM_INT__i_s_i_s, NULL, PRIM_INT__k_s_k_s);
853 
854 
855                     ostei_general_vrr_I(5, 0, 10, 0, 3,
856                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
857                             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);
858 
859 
860                     ostei_general_vrr_I(6, 0, 9, 0, 2,
861                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
862                             PRIM_INT__h_s_m_s, PRIM_INT__g_s_m_s, NULL, PRIM_INT__h_s_l_s, NULL, PRIM_INT__i_s_m_s);
863 
864 
865                     ostei_general_vrr_I(7, 0, 8, 0, 1,
866                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
867                             PRIM_INT__i_s_l_s, PRIM_INT__h_s_l_s, NULL, PRIM_INT__i_s_k_s, NULL, PRIM_INT__k_s_l_s);
868 
869 
870                     ostei_general_vrr_I(6, 0, 10, 0, 2,
871                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
872                             PRIM_INT__h_s_n_s, PRIM_INT__g_s_n_s, NULL, PRIM_INT__h_s_m_s, NULL, PRIM_INT__i_s_n_s);
873 
874 
875                     ostei_general_vrr_I(7, 0, 9, 0, 1,
876                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
877                             PRIM_INT__i_s_m_s, PRIM_INT__h_s_m_s, NULL, PRIM_INT__i_s_l_s, NULL, PRIM_INT__k_s_m_s);
878 
879 
880                     ostei_general_vrr_I(7, 0, 10, 0, 1,
881                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
882                             PRIM_INT__i_s_n_s, PRIM_INT__h_s_n_s, NULL, PRIM_INT__i_s_m_s, NULL, PRIM_INT__k_s_n_s);
883 
884 
885 
886 
887                     ////////////////////////////////////
888                     // Accumulate contracted integrals
889                     ////////////////////////////////////
890                     if(lastoffset == 0)
891                     {
892                         contract_all(540, PRIM_INT__g_s_k_s, PRIM_PTR_INT__g_s_k_s);
893                         contract_all(675, PRIM_INT__g_s_l_s, PRIM_PTR_INT__g_s_l_s);
894                         contract_all(825, PRIM_INT__g_s_m_s, PRIM_PTR_INT__g_s_m_s);
895                         contract_all(990, PRIM_INT__g_s_n_s, PRIM_PTR_INT__g_s_n_s);
896                         contract_all(756, PRIM_INT__h_s_k_s, PRIM_PTR_INT__h_s_k_s);
897                         contract_all(945, PRIM_INT__h_s_l_s, PRIM_PTR_INT__h_s_l_s);
898                         contract_all(1155, PRIM_INT__h_s_m_s, PRIM_PTR_INT__h_s_m_s);
899                         contract_all(1386, PRIM_INT__h_s_n_s, PRIM_PTR_INT__h_s_n_s);
900                         contract_all(1008, PRIM_INT__i_s_k_s, PRIM_PTR_INT__i_s_k_s);
901                         contract_all(1260, PRIM_INT__i_s_l_s, PRIM_PTR_INT__i_s_l_s);
902                         contract_all(1540, PRIM_INT__i_s_m_s, PRIM_PTR_INT__i_s_m_s);
903                         contract_all(1848, PRIM_INT__i_s_n_s, PRIM_PTR_INT__i_s_n_s);
904                         contract_all(1296, PRIM_INT__k_s_k_s, PRIM_PTR_INT__k_s_k_s);
905                         contract_all(1620, PRIM_INT__k_s_l_s, PRIM_PTR_INT__k_s_l_s);
906                         contract_all(1980, PRIM_INT__k_s_m_s, PRIM_PTR_INT__k_s_m_s);
907                         contract_all(2376, PRIM_INT__k_s_n_s, PRIM_PTR_INT__k_s_n_s);
908                     }
909                     else
910                     {
911                         contract(540, shelloffsets, PRIM_INT__g_s_k_s, PRIM_PTR_INT__g_s_k_s);
912                         contract(675, shelloffsets, PRIM_INT__g_s_l_s, PRIM_PTR_INT__g_s_l_s);
913                         contract(825, shelloffsets, PRIM_INT__g_s_m_s, PRIM_PTR_INT__g_s_m_s);
914                         contract(990, shelloffsets, PRIM_INT__g_s_n_s, PRIM_PTR_INT__g_s_n_s);
915                         contract(756, shelloffsets, PRIM_INT__h_s_k_s, PRIM_PTR_INT__h_s_k_s);
916                         contract(945, shelloffsets, PRIM_INT__h_s_l_s, PRIM_PTR_INT__h_s_l_s);
917                         contract(1155, shelloffsets, PRIM_INT__h_s_m_s, PRIM_PTR_INT__h_s_m_s);
918                         contract(1386, shelloffsets, PRIM_INT__h_s_n_s, PRIM_PTR_INT__h_s_n_s);
919                         contract(1008, shelloffsets, PRIM_INT__i_s_k_s, PRIM_PTR_INT__i_s_k_s);
920                         contract(1260, shelloffsets, PRIM_INT__i_s_l_s, PRIM_PTR_INT__i_s_l_s);
921                         contract(1540, shelloffsets, PRIM_INT__i_s_m_s, PRIM_PTR_INT__i_s_m_s);
922                         contract(1848, shelloffsets, PRIM_INT__i_s_n_s, PRIM_PTR_INT__i_s_n_s);
923                         contract(1296, shelloffsets, PRIM_INT__k_s_k_s, PRIM_PTR_INT__k_s_k_s);
924                         contract(1620, shelloffsets, PRIM_INT__k_s_l_s, PRIM_PTR_INT__k_s_l_s);
925                         contract(1980, shelloffsets, PRIM_INT__k_s_m_s, PRIM_PTR_INT__k_s_m_s);
926                         contract(2376, shelloffsets, PRIM_INT__k_s_n_s, PRIM_PTR_INT__k_s_n_s);
927                         PRIM_PTR_INT__g_s_k_s += lastoffset*540;
928                         PRIM_PTR_INT__g_s_l_s += lastoffset*675;
929                         PRIM_PTR_INT__g_s_m_s += lastoffset*825;
930                         PRIM_PTR_INT__g_s_n_s += lastoffset*990;
931                         PRIM_PTR_INT__h_s_k_s += lastoffset*756;
932                         PRIM_PTR_INT__h_s_l_s += lastoffset*945;
933                         PRIM_PTR_INT__h_s_m_s += lastoffset*1155;
934                         PRIM_PTR_INT__h_s_n_s += lastoffset*1386;
935                         PRIM_PTR_INT__i_s_k_s += lastoffset*1008;
936                         PRIM_PTR_INT__i_s_l_s += lastoffset*1260;
937                         PRIM_PTR_INT__i_s_m_s += lastoffset*1540;
938                         PRIM_PTR_INT__i_s_n_s += lastoffset*1848;
939                         PRIM_PTR_INT__k_s_k_s += lastoffset*1296;
940                         PRIM_PTR_INT__k_s_l_s += lastoffset*1620;
941                         PRIM_PTR_INT__k_s_m_s += lastoffset*1980;
942                         PRIM_PTR_INT__k_s_n_s += lastoffset*2376;
943                     }
944 
945                 }  // close loop over j
946             }  // close loop over i
947 
948             //Advance to the next batch
949             jstart = SIMINT_SIMD_ROUND(jend);
950 
951             //////////////////////////////////////////////
952             // Contracted integrals: Horizontal recurrance
953             //////////////////////////////////////////////
954 
955 
956             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
957 
958 
959             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
960             {
961                 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
962 
963                 // set up HRR pointers
964                 double const * restrict HRR_INT__g_s_k_s = INT__g_s_k_s + abcd * 540;
965                 double const * restrict HRR_INT__g_s_l_s = INT__g_s_l_s + abcd * 675;
966                 double const * restrict HRR_INT__g_s_m_s = INT__g_s_m_s + abcd * 825;
967                 double const * restrict HRR_INT__g_s_n_s = INT__g_s_n_s + abcd * 990;
968                 double const * restrict HRR_INT__h_s_k_s = INT__h_s_k_s + abcd * 756;
969                 double const * restrict HRR_INT__h_s_l_s = INT__h_s_l_s + abcd * 945;
970                 double const * restrict HRR_INT__h_s_m_s = INT__h_s_m_s + abcd * 1155;
971                 double const * restrict HRR_INT__h_s_n_s = INT__h_s_n_s + abcd * 1386;
972                 double const * restrict HRR_INT__i_s_k_s = INT__i_s_k_s + abcd * 1008;
973                 double const * restrict HRR_INT__i_s_l_s = INT__i_s_l_s + abcd * 1260;
974                 double const * restrict HRR_INT__i_s_m_s = INT__i_s_m_s + abcd * 1540;
975                 double const * restrict HRR_INT__i_s_n_s = INT__i_s_n_s + abcd * 1848;
976                 double const * restrict HRR_INT__k_s_k_s = INT__k_s_k_s + abcd * 1296;
977                 double const * restrict HRR_INT__k_s_l_s = INT__k_s_l_s + abcd * 1620;
978                 double const * restrict HRR_INT__k_s_m_s = INT__k_s_m_s + abcd * 1980;
979                 double const * restrict HRR_INT__k_s_n_s = INT__k_s_n_s + abcd * 2376;
980                 double * restrict HRR_INT__g_f_k_f = INT__g_f_k_f + real_abcd * 54000;
981 
982                 // form INT__g_p_k_s
983                 HRR_J_g_p(
984                     HRR_INT__g_p_k_s,
985                     HRR_INT__g_s_k_s,
986                     HRR_INT__h_s_k_s,
987                     hAB, 36);
988 
989                 // form INT__g_p_l_s
990                 HRR_J_g_p(
991                     HRR_INT__g_p_l_s,
992                     HRR_INT__g_s_l_s,
993                     HRR_INT__h_s_l_s,
994                     hAB, 45);
995 
996                 // form INT__g_p_m_s
997                 HRR_J_g_p(
998                     HRR_INT__g_p_m_s,
999                     HRR_INT__g_s_m_s,
1000                     HRR_INT__h_s_m_s,
1001                     hAB, 55);
1002 
1003                 // form INT__g_p_n_s
1004                 HRR_J_g_p(
1005                     HRR_INT__g_p_n_s,
1006                     HRR_INT__g_s_n_s,
1007                     HRR_INT__h_s_n_s,
1008                     hAB, 66);
1009 
1010                 // form INT__h_p_k_s
1011                 ostei_general_hrr_J(5, 1, 7, 0, hAB, HRR_INT__i_s_k_s, HRR_INT__h_s_k_s, HRR_INT__h_p_k_s);
1012 
1013                 // form INT__h_p_l_s
1014                 ostei_general_hrr_J(5, 1, 8, 0, hAB, HRR_INT__i_s_l_s, HRR_INT__h_s_l_s, HRR_INT__h_p_l_s);
1015 
1016                 // form INT__h_p_m_s
1017                 ostei_general_hrr_J(5, 1, 9, 0, hAB, HRR_INT__i_s_m_s, HRR_INT__h_s_m_s, HRR_INT__h_p_m_s);
1018 
1019                 // form INT__h_p_n_s
1020                 ostei_general_hrr_J(5, 1, 10, 0, hAB, HRR_INT__i_s_n_s, HRR_INT__h_s_n_s, HRR_INT__h_p_n_s);
1021 
1022                 // form INT__i_p_k_s
1023                 ostei_general_hrr_J(6, 1, 7, 0, hAB, HRR_INT__k_s_k_s, HRR_INT__i_s_k_s, HRR_INT__i_p_k_s);
1024 
1025                 // form INT__i_p_l_s
1026                 ostei_general_hrr_J(6, 1, 8, 0, hAB, HRR_INT__k_s_l_s, HRR_INT__i_s_l_s, HRR_INT__i_p_l_s);
1027 
1028                 // form INT__i_p_m_s
1029                 ostei_general_hrr_J(6, 1, 9, 0, hAB, HRR_INT__k_s_m_s, HRR_INT__i_s_m_s, HRR_INT__i_p_m_s);
1030 
1031                 // form INT__i_p_n_s
1032                 ostei_general_hrr_J(6, 1, 10, 0, hAB, HRR_INT__k_s_n_s, HRR_INT__i_s_n_s, HRR_INT__i_p_n_s);
1033 
1034                 // form INT__g_d_k_s
1035                 ostei_general_hrr_J(4, 2, 7, 0, hAB, HRR_INT__h_p_k_s, HRR_INT__g_p_k_s, HRR_INT__g_d_k_s);
1036 
1037                 // form INT__g_d_l_s
1038                 ostei_general_hrr_J(4, 2, 8, 0, hAB, HRR_INT__h_p_l_s, HRR_INT__g_p_l_s, HRR_INT__g_d_l_s);
1039 
1040                 // form INT__g_d_m_s
1041                 ostei_general_hrr_J(4, 2, 9, 0, hAB, HRR_INT__h_p_m_s, HRR_INT__g_p_m_s, HRR_INT__g_d_m_s);
1042 
1043                 // form INT__g_d_n_s
1044                 ostei_general_hrr_J(4, 2, 10, 0, hAB, HRR_INT__h_p_n_s, HRR_INT__g_p_n_s, HRR_INT__g_d_n_s);
1045 
1046                 // form INT__h_d_k_s
1047                 ostei_general_hrr_J(5, 2, 7, 0, hAB, HRR_INT__i_p_k_s, HRR_INT__h_p_k_s, HRR_INT__h_d_k_s);
1048 
1049                 // form INT__h_d_l_s
1050                 ostei_general_hrr_J(5, 2, 8, 0, hAB, HRR_INT__i_p_l_s, HRR_INT__h_p_l_s, HRR_INT__h_d_l_s);
1051 
1052                 // form INT__h_d_m_s
1053                 ostei_general_hrr_J(5, 2, 9, 0, hAB, HRR_INT__i_p_m_s, HRR_INT__h_p_m_s, HRR_INT__h_d_m_s);
1054 
1055                 // form INT__h_d_n_s
1056                 ostei_general_hrr_J(5, 2, 10, 0, hAB, HRR_INT__i_p_n_s, HRR_INT__h_p_n_s, HRR_INT__h_d_n_s);
1057 
1058                 // form INT__g_f_k_s
1059                 ostei_general_hrr_J(4, 3, 7, 0, hAB, HRR_INT__h_d_k_s, HRR_INT__g_d_k_s, HRR_INT__g_f_k_s);
1060 
1061                 // form INT__g_f_l_s
1062                 ostei_general_hrr_J(4, 3, 8, 0, hAB, HRR_INT__h_d_l_s, HRR_INT__g_d_l_s, HRR_INT__g_f_l_s);
1063 
1064                 // form INT__g_f_m_s
1065                 ostei_general_hrr_J(4, 3, 9, 0, hAB, HRR_INT__h_d_m_s, HRR_INT__g_d_m_s, HRR_INT__g_f_m_s);
1066 
1067                 // form INT__g_f_n_s
1068                 ostei_general_hrr_J(4, 3, 10, 0, hAB, HRR_INT__h_d_n_s, HRR_INT__g_d_n_s, HRR_INT__g_f_n_s);
1069 
1070                 // form INT__g_f_k_p
1071                 ostei_general_hrr_L(4, 3, 7, 1, hCD, HRR_INT__g_f_l_s, HRR_INT__g_f_k_s, HRR_INT__g_f_k_p);
1072 
1073                 // form INT__g_f_l_p
1074                 ostei_general_hrr_L(4, 3, 8, 1, hCD, HRR_INT__g_f_m_s, HRR_INT__g_f_l_s, HRR_INT__g_f_l_p);
1075 
1076                 // form INT__g_f_m_p
1077                 ostei_general_hrr_L(4, 3, 9, 1, hCD, HRR_INT__g_f_n_s, HRR_INT__g_f_m_s, HRR_INT__g_f_m_p);
1078 
1079                 // form INT__g_f_k_d
1080                 ostei_general_hrr_L(4, 3, 7, 2, hCD, HRR_INT__g_f_l_p, HRR_INT__g_f_k_p, HRR_INT__g_f_k_d);
1081 
1082                 // form INT__g_f_l_d
1083                 ostei_general_hrr_L(4, 3, 8, 2, hCD, HRR_INT__g_f_m_p, HRR_INT__g_f_l_p, HRR_INT__g_f_l_d);
1084 
1085                 // form INT__g_f_k_f
1086                 ostei_general_hrr_L(4, 3, 7, 3, hCD, HRR_INT__g_f_l_d, HRR_INT__g_f_k_d, HRR_INT__g_f_k_f);
1087 
1088 
1089             }  // close HRR loop
1090 
1091 
1092         }   // close loop cdbatch
1093 
1094         istart = iend;
1095     }  // close loop over ab
1096 
1097     return P.nshell12_clip * Q.nshell12_clip;
1098 }
1099 
ostei_f_g_k_f(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__f_g_k_f)1100 int ostei_f_g_k_f(struct simint_multi_shellpair const P,
1101                   struct simint_multi_shellpair const Q,
1102                   double screen_tol,
1103                   double * const restrict work,
1104                   double * const restrict INT__f_g_k_f)
1105 {
1106     double P_AB[3*P.nshell12];
1107     struct simint_multi_shellpair P_tmp = P;
1108     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
1109     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
1110     P_tmp.AB_x = P_AB;
1111     P_tmp.AB_y = P_AB + P.nshell12;
1112     P_tmp.AB_z = P_AB + 2*P.nshell12;
1113 
1114     for(int i = 0; i < P.nshell12; i++)
1115     {
1116         P_tmp.AB_x[i] = -P.AB_x[i];
1117         P_tmp.AB_y[i] = -P.AB_y[i];
1118         P_tmp.AB_z[i] = -P.AB_z[i];
1119     }
1120 
1121     int ret = ostei_g_f_k_f(P_tmp, Q, screen_tol, work, INT__f_g_k_f);
1122     double buffer[54000] SIMINT_ALIGN_ARRAY_DBL;
1123 
1124     for(int q = 0; q < ret; q++)
1125     {
1126         int idx = 0;
1127         for(int a = 0; a < 10; ++a)
1128         for(int b = 0; b < 15; ++b)
1129         for(int c = 0; c < 36; ++c)
1130         for(int d = 0; d < 10; ++d)
1131             buffer[idx++] = INT__f_g_k_f[q*54000+b*3600+a*360+c*10+d];
1132 
1133         memcpy(INT__f_g_k_f+q*54000, buffer, 54000*sizeof(double));
1134     }
1135 
1136     return ret;
1137 }
1138 
ostei_g_f_f_k(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__g_f_f_k)1139 int ostei_g_f_f_k(struct simint_multi_shellpair const P,
1140                   struct simint_multi_shellpair const Q,
1141                   double screen_tol,
1142                   double * const restrict work,
1143                   double * const restrict INT__g_f_f_k)
1144 {
1145     double Q_AB[3*Q.nshell12];
1146     struct simint_multi_shellpair Q_tmp = Q;
1147     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
1148     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
1149     Q_tmp.AB_x = Q_AB;
1150     Q_tmp.AB_y = Q_AB + Q.nshell12;
1151     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
1152 
1153     for(int i = 0; i < Q.nshell12; i++)
1154     {
1155         Q_tmp.AB_x[i] = -Q.AB_x[i];
1156         Q_tmp.AB_y[i] = -Q.AB_y[i];
1157         Q_tmp.AB_z[i] = -Q.AB_z[i];
1158     }
1159 
1160     int ret = ostei_g_f_k_f(P, Q_tmp, screen_tol, work, INT__g_f_f_k);
1161     double buffer[54000] SIMINT_ALIGN_ARRAY_DBL;
1162 
1163     for(int q = 0; q < ret; q++)
1164     {
1165         int idx = 0;
1166         for(int a = 0; a < 15; ++a)
1167         for(int b = 0; b < 10; ++b)
1168         for(int c = 0; c < 10; ++c)
1169         for(int d = 0; d < 36; ++d)
1170             buffer[idx++] = INT__g_f_f_k[q*54000+a*3600+b*360+d*10+c];
1171 
1172         memcpy(INT__g_f_f_k+q*54000, buffer, 54000*sizeof(double));
1173     }
1174 
1175     return ret;
1176 }
1177 
ostei_f_g_f_k(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__f_g_f_k)1178 int ostei_f_g_f_k(struct simint_multi_shellpair const P,
1179                   struct simint_multi_shellpair const Q,
1180                   double screen_tol,
1181                   double * const restrict work,
1182                   double * const restrict INT__f_g_f_k)
1183 {
1184     double P_AB[3*P.nshell12];
1185     struct simint_multi_shellpair P_tmp = P;
1186     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
1187     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
1188     P_tmp.AB_x = P_AB;
1189     P_tmp.AB_y = P_AB + P.nshell12;
1190     P_tmp.AB_z = P_AB + 2*P.nshell12;
1191 
1192     for(int i = 0; i < P.nshell12; i++)
1193     {
1194         P_tmp.AB_x[i] = -P.AB_x[i];
1195         P_tmp.AB_y[i] = -P.AB_y[i];
1196         P_tmp.AB_z[i] = -P.AB_z[i];
1197     }
1198 
1199     double Q_AB[3*Q.nshell12];
1200     struct simint_multi_shellpair Q_tmp = Q;
1201     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
1202     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
1203     Q_tmp.AB_x = Q_AB;
1204     Q_tmp.AB_y = Q_AB + Q.nshell12;
1205     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
1206 
1207     for(int i = 0; i < Q.nshell12; i++)
1208     {
1209         Q_tmp.AB_x[i] = -Q.AB_x[i];
1210         Q_tmp.AB_y[i] = -Q.AB_y[i];
1211         Q_tmp.AB_z[i] = -Q.AB_z[i];
1212     }
1213 
1214     int ret = ostei_g_f_k_f(P_tmp, Q_tmp, screen_tol, work, INT__f_g_f_k);
1215     double buffer[54000] SIMINT_ALIGN_ARRAY_DBL;
1216 
1217     for(int q = 0; q < ret; q++)
1218     {
1219         int idx = 0;
1220         for(int a = 0; a < 10; ++a)
1221         for(int b = 0; b < 15; ++b)
1222         for(int c = 0; c < 10; ++c)
1223         for(int d = 0; d < 36; ++d)
1224             buffer[idx++] = INT__f_g_f_k[q*54000+b*3600+a*360+d*10+c];
1225 
1226         memcpy(INT__f_g_f_k+q*54000, buffer, 54000*sizeof(double));
1227     }
1228 
1229     return ret;
1230 }
1231 
1232