1 #include "simint/boys/boys.h"
2 #include "simint/ostei/gen/ostei_generated.h"
3 #include "simint/vectorization/vectorization.h"
4 #include <math.h>
5 #include <string.h>
6 
7 
ostei_k_g_g_p(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__k_g_g_p)8 int ostei_k_g_g_p(struct simint_multi_shellpair const P,
9                   struct simint_multi_shellpair const Q,
10                   double screen_tol,
11                   double * const restrict work,
12                   double * const restrict INT__k_g_g_p)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__k_g_g_p);
17     int ab, cd, abcd;
18     int istart, jstart;
19     int iprimcd, nprim_icd, icd;
20     const int check_screen = (screen_tol > 0.0);
21     int i, j;
22     int n;
23     int not_screened;
24     int real_abcd;
25     int iket;
26     int ibra;
27 
28     // partition workspace
29     double * const INT__k_s_g_s = work + (SIMINT_NSHELL_SIMD * 0);
30     double * const INT__k_s_h_s = work + (SIMINT_NSHELL_SIMD * 540);
31     double * const INT__l_s_g_s = work + (SIMINT_NSHELL_SIMD * 1296);
32     double * const INT__l_s_h_s = work + (SIMINT_NSHELL_SIMD * 1971);
33     double * const INT__m_s_g_s = work + (SIMINT_NSHELL_SIMD * 2916);
34     double * const INT__m_s_h_s = work + (SIMINT_NSHELL_SIMD * 3741);
35     double * const INT__n_s_g_s = work + (SIMINT_NSHELL_SIMD * 4896);
36     double * const INT__n_s_h_s = work + (SIMINT_NSHELL_SIMD * 5886);
37     double * const INT__o_s_g_s = work + (SIMINT_NSHELL_SIMD * 7272);
38     double * const INT__o_s_h_s = work + (SIMINT_NSHELL_SIMD * 8442);
39     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*10080);
40     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
41     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_s = primwork + 17;
42     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 65;
43     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_s = primwork + 155;
44     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_p_s = primwork + 295;
45     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_s_s = primwork + 445;
46     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_p_s = primwork + 640;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_d_s = primwork + 865;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_s_s = primwork + 1225;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_p_s = primwork + 1477;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_d_s = primwork + 1792;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_f_s = primwork + 2296;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_s_s = primwork + 2926;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_p_s = primwork + 3234;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_d_s = primwork + 3654;
55     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_f_s = primwork + 4326;
56     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_g_s = primwork + 5166;
57     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_s_s = primwork + 6006;
58     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_p_s = primwork + 6366;
59     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_d_s = primwork + 6906;
60     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_f_s = primwork + 7770;
61     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_g_s = primwork + 8850;
62     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_h_s = primwork + 9930;
63     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_s_s = primwork + 10686;
64     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_p_s = primwork + 11091;
65     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_d_s = primwork + 11766;
66     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_f_s = primwork + 12846;
67     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_g_s = primwork + 14196;
68     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_h_s = primwork + 15546;
69     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_s_s = primwork + 16491;
70     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_p_s = primwork + 16931;
71     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_d_s = primwork + 17756;
72     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_f_s = primwork + 19076;
73     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_g_s = primwork + 20726;
74     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_h_s = primwork + 22376;
75     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_s_s = primwork + 23531;
76     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_p_s = primwork + 23993;
77     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_d_s = primwork + 24983;
78     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_f_s = primwork + 26567;
79     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_g_s = primwork + 28547;
80     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_h_s = primwork + 30527;
81     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_s_s = primwork + 31913;
82     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_p_s = primwork + 32381;
83     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_d_s = primwork + 33551;
84     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_f_s = primwork + 35423;
85     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_g_s = primwork + 37763;
86     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_h_s = primwork + 40103;
87     double * const hrrwork = (double *)(primwork + 41741);
88     double * const HRR_INT__k_p_g_s = hrrwork + 0;
89     double * const HRR_INT__k_p_h_s = hrrwork + 1620;
90     double * const HRR_INT__k_d_g_s = hrrwork + 3888;
91     double * const HRR_INT__k_d_h_s = hrrwork + 7128;
92     double * const HRR_INT__k_f_g_s = hrrwork + 11664;
93     double * const HRR_INT__k_f_h_s = hrrwork + 17064;
94     double * const HRR_INT__k_g_g_s = hrrwork + 24624;
95     double * const HRR_INT__k_g_h_s = hrrwork + 32724;
96     double * const HRR_INT__l_p_g_s = hrrwork + 44064;
97     double * const HRR_INT__l_p_h_s = hrrwork + 46089;
98     double * const HRR_INT__l_d_g_s = hrrwork + 48924;
99     double * const HRR_INT__l_d_h_s = hrrwork + 52974;
100     double * const HRR_INT__l_f_g_s = hrrwork + 58644;
101     double * const HRR_INT__l_f_h_s = hrrwork + 65394;
102     double * const HRR_INT__m_p_g_s = hrrwork + 74844;
103     double * const HRR_INT__m_p_h_s = hrrwork + 77319;
104     double * const HRR_INT__m_d_g_s = hrrwork + 80784;
105     double * const HRR_INT__m_d_h_s = hrrwork + 85734;
106     double * const HRR_INT__n_p_g_s = hrrwork + 92664;
107     double * const HRR_INT__n_p_h_s = hrrwork + 95634;
108 
109 
110     // Create constants
111     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
112     const SIMINT_DBLTYPE const_10 = SIMINT_DBLSET1(10);
113     const SIMINT_DBLTYPE const_11 = SIMINT_DBLSET1(11);
114     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
115     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
116     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
117     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
118     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
119     const SIMINT_DBLTYPE const_7 = SIMINT_DBLSET1(7);
120     const SIMINT_DBLTYPE const_8 = SIMINT_DBLSET1(8);
121     const SIMINT_DBLTYPE const_9 = SIMINT_DBLSET1(9);
122     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
123 
124 
125     ////////////////////////////////////////
126     // Loop over shells and primitives
127     ////////////////////////////////////////
128 
129     real_abcd = 0;
130     istart = 0;
131     for(ab = 0; ab < P.nshell12_clip; ++ab)
132     {
133         const int iend = istart + P.nprim12[ab];
134 
135         cd = 0;
136         jstart = 0;
137 
138         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
139         {
140             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
141             int jend = jstart;
142             for(i = 0; i < nshellbatch; i++)
143                 jend += Q.nprim12[cd+i];
144 
145             // Clear the beginning of the workspace (where we are accumulating integrals)
146             memset(work, 0, SIMINT_NSHELL_SIMD * 10080 * sizeof(double));
147             abcd = 0;
148 
149 
150             for(i = istart; i < iend; ++i)
151             {
152                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
153 
154                 if(check_screen)
155                 {
156                     // Skip this whole thing if always insignificant
157                     if((P.screen[i] * Q.screen_max) < screen_tol)
158                         continue;
159                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
160                 }
161 
162                 icd = 0;
163                 iprimcd = 0;
164                 nprim_icd = Q.nprim12[cd];
165                 double * restrict PRIM_PTR_INT__k_s_g_s = INT__k_s_g_s + abcd * 540;
166                 double * restrict PRIM_PTR_INT__k_s_h_s = INT__k_s_h_s + abcd * 756;
167                 double * restrict PRIM_PTR_INT__l_s_g_s = INT__l_s_g_s + abcd * 675;
168                 double * restrict PRIM_PTR_INT__l_s_h_s = INT__l_s_h_s + abcd * 945;
169                 double * restrict PRIM_PTR_INT__m_s_g_s = INT__m_s_g_s + abcd * 825;
170                 double * restrict PRIM_PTR_INT__m_s_h_s = INT__m_s_h_s + abcd * 1155;
171                 double * restrict PRIM_PTR_INT__n_s_g_s = INT__n_s_g_s + abcd * 990;
172                 double * restrict PRIM_PTR_INT__n_s_h_s = INT__n_s_h_s + abcd * 1386;
173                 double * restrict PRIM_PTR_INT__o_s_g_s = INT__o_s_g_s + abcd * 1170;
174                 double * restrict PRIM_PTR_INT__o_s_h_s = INT__o_s_h_s + abcd * 1638;
175 
176 
177 
178                 // Load these one per loop over i
179                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
180                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
181                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
182 
183                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
184 
185                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
186                 {
187                     // calculate the shell offsets
188                     // these are the offset from the shell pointed to by cd
189                     // for each element
190                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
191                     int lastoffset = 0;
192                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
193 
194                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
195                     {
196                         // Handle if the first element of the vector is a new shell
197                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
198                         {
199                             nprim_icd += Q.nprim12[cd + (++icd)];
200                             PRIM_PTR_INT__k_s_g_s += 540;
201                             PRIM_PTR_INT__k_s_h_s += 756;
202                             PRIM_PTR_INT__l_s_g_s += 675;
203                             PRIM_PTR_INT__l_s_h_s += 945;
204                             PRIM_PTR_INT__m_s_g_s += 825;
205                             PRIM_PTR_INT__m_s_h_s += 1155;
206                             PRIM_PTR_INT__n_s_g_s += 990;
207                             PRIM_PTR_INT__n_s_h_s += 1386;
208                             PRIM_PTR_INT__o_s_g_s += 1170;
209                             PRIM_PTR_INT__o_s_h_s += 1638;
210                         }
211                         iprimcd++;
212                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
213                         {
214                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
215                             {
216                                 shelloffsets[n] = shelloffsets[n-1] + 1;
217                                 lastoffset++;
218                                 nprim_icd += Q.nprim12[cd + (++icd)];
219                             }
220                             else
221                                 shelloffsets[n] = shelloffsets[n-1];
222                             iprimcd++;
223                         }
224                     }
225                     else
226                         iprimcd += SIMINT_SIMD_LEN;
227 
228                     // Do we have to compute this vector (or has it been screened out)?
229                     // (not_screened != 0 means we have to do this vector)
230                     if(check_screen)
231                     {
232                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
233                         if(vmax < screen_tol)
234                         {
235                             PRIM_PTR_INT__k_s_g_s += lastoffset*540;
236                             PRIM_PTR_INT__k_s_h_s += lastoffset*756;
237                             PRIM_PTR_INT__l_s_g_s += lastoffset*675;
238                             PRIM_PTR_INT__l_s_h_s += lastoffset*945;
239                             PRIM_PTR_INT__m_s_g_s += lastoffset*825;
240                             PRIM_PTR_INT__m_s_h_s += lastoffset*1155;
241                             PRIM_PTR_INT__n_s_g_s += lastoffset*990;
242                             PRIM_PTR_INT__n_s_h_s += lastoffset*1386;
243                             PRIM_PTR_INT__o_s_g_s += lastoffset*1170;
244                             PRIM_PTR_INT__o_s_h_s += lastoffset*1638;
245                             continue;
246                         }
247                     }
248 
249                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
250                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
251                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
252                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
253 
254 
255                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
256                     SIMINT_DBLTYPE PQ[3];
257                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
258                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
259                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
260                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
261                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
262                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
263 
264                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
265                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
266                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
267                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
268                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
269                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
270                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
271 
272                     // NOTE: Minus sign!
273                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
274                     SIMINT_DBLTYPE aop_PQ[3];
275                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
276                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
277                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
278 
279                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
280                     SIMINT_DBLTYPE aoq_PQ[3];
281                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
282                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
283                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
284                     // Put a minus sign here so we don't have to in RR routines
285                     a_over_q = SIMINT_NEG(a_over_q);
286 
287 
288                     //////////////////////////////////////////////
289                     // Fjt function section
290                     // Maximum v value: 16
291                     //////////////////////////////////////////////
292                     // The parameter to the Fjt function
293                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
294 
295 
296                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
297 
298 
299                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 16);
300                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
301                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
302                     for(n = 0; n <= 16; n++)
303                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
304 
305                     //////////////////////////////////////////////
306                     // Primitive integrals: Vertical recurrance
307                     //////////////////////////////////////////////
308 
309                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
310                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
311                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
312                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
313                     const SIMINT_DBLTYPE vrr_const_5_over_2p = SIMINT_MUL(const_5, one_over_2p);
314                     const SIMINT_DBLTYPE vrr_const_6_over_2p = SIMINT_MUL(const_6, one_over_2p);
315                     const SIMINT_DBLTYPE vrr_const_7_over_2p = SIMINT_MUL(const_7, one_over_2p);
316                     const SIMINT_DBLTYPE vrr_const_8_over_2p = SIMINT_MUL(const_8, one_over_2p);
317                     const SIMINT_DBLTYPE vrr_const_9_over_2p = SIMINT_MUL(const_9, one_over_2p);
318                     const SIMINT_DBLTYPE vrr_const_10_over_2p = SIMINT_MUL(const_10, one_over_2p);
319                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
320                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
321                     const SIMINT_DBLTYPE vrr_const_3_over_2q = SIMINT_MUL(const_3, one_over_2q);
322                     const SIMINT_DBLTYPE vrr_const_4_over_2q = SIMINT_MUL(const_4, one_over_2q);
323                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
324                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
325                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
326                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
327                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
328                     const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
329                     const SIMINT_DBLTYPE vrr_const_7_over_2pq = SIMINT_MUL(const_7, one_over_2pq);
330                     const SIMINT_DBLTYPE vrr_const_8_over_2pq = SIMINT_MUL(const_8, one_over_2pq);
331                     const SIMINT_DBLTYPE vrr_const_9_over_2pq = SIMINT_MUL(const_9, one_over_2pq);
332                     const SIMINT_DBLTYPE vrr_const_10_over_2pq = SIMINT_MUL(const_10, one_over_2pq);
333                     const SIMINT_DBLTYPE vrr_const_11_over_2pq = SIMINT_MUL(const_11, one_over_2pq);
334 
335 
336 
337                     // Forming PRIM_INT__p_s_s_s[16 * 3];
338                     for(n = 0; n < 16; ++n)  // loop over orders of auxiliary function
339                     {
340 
341                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
342                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]);
343 
344                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
345                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_s[n * 3 + 1]);
346 
347                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
348                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_s[n * 3 + 2]);
349 
350                     }
351 
352 
353 
354                     // Forming PRIM_INT__d_s_s_s[15 * 6];
355                     for(n = 0; n < 15; ++n)  // loop over orders of auxiliary function
356                     {
357 
358                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
359                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_s_s[n * 6 + 0]);
360                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__d_s_s_s[n * 6 + 0]);
361 
362                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
363                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_s_s[n * 6 + 1]);
364 
365                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
366                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_s_s[n * 6 + 2]);
367 
368                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
369                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_s_s[n * 6 + 3]);
370                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__d_s_s_s[n * 6 + 3]);
371 
372                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
373                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_s_s[n * 6 + 4]);
374 
375                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
376                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_s_s[n * 6 + 5]);
377                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__d_s_s_s[n * 6 + 5]);
378 
379                     }
380 
381 
382 
383                     // Forming PRIM_INT__f_s_s_s[14 * 10];
384                     for(n = 0; n < 14; ++n)  // loop over orders of auxiliary function
385                     {
386 
387                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
388                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__f_s_s_s[n * 10 + 0]);
389                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]), PRIM_INT__f_s_s_s[n * 10 + 0]);
390 
391                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
392                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__f_s_s_s[n * 10 + 1]);
393 
394                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
395                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__f_s_s_s[n * 10 + 2]);
396 
397                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
398                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__f_s_s_s[n * 10 + 3]);
399 
400                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
401                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__f_s_s_s[n * 10 + 4]);
402 
403                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
404                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__f_s_s_s[n * 10 + 5]);
405 
406                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
407                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__f_s_s_s[n * 10 + 6]);
408                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_s[n * 3 + 1]), PRIM_INT__f_s_s_s[n * 10 + 6]);
409 
410                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
411                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__f_s_s_s[n * 10 + 7]);
412 
413                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
414                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__f_s_s_s[n * 10 + 8]);
415 
416                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
417                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__f_s_s_s[n * 10 + 9]);
418                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_s[n * 3 + 2]), PRIM_INT__f_s_s_s[n * 10 + 9]);
419 
420                     }
421 
422 
423                     VRR_I_g_s_s_s(
424                             PRIM_INT__g_s_s_s,
425                             PRIM_INT__f_s_s_s,
426                             PRIM_INT__d_s_s_s,
427                             P_PA,
428                             a_over_p,
429                             aop_PQ,
430                             one_over_2p,
431                             13);
432 
433 
434                     VRR_I_h_s_s_s(
435                             PRIM_INT__h_s_s_s,
436                             PRIM_INT__g_s_s_s,
437                             PRIM_INT__f_s_s_s,
438                             P_PA,
439                             a_over_p,
440                             aop_PQ,
441                             one_over_2p,
442                             12);
443 
444 
445                     ostei_general_vrr1_I(6, 11,
446                             one_over_2p, a_over_p, aop_PQ, P_PA,
447                             PRIM_INT__h_s_s_s, PRIM_INT__g_s_s_s, PRIM_INT__i_s_s_s);
448 
449 
450                     ostei_general_vrr1_I(7, 10,
451                             one_over_2p, a_over_p, aop_PQ, P_PA,
452                             PRIM_INT__i_s_s_s, PRIM_INT__h_s_s_s, PRIM_INT__k_s_s_s);
453 
454 
455                     ostei_general_vrr_K(7, 0, 1, 0, 5,
456                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
457                             PRIM_INT__k_s_s_s, NULL, NULL, PRIM_INT__i_s_s_s, NULL, PRIM_INT__k_s_p_s);
458 
459 
460                     ostei_general_vrr_K(6, 0, 1, 0, 5,
461                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
462                             PRIM_INT__i_s_s_s, NULL, NULL, PRIM_INT__h_s_s_s, NULL, PRIM_INT__i_s_p_s);
463 
464 
465                     ostei_general_vrr_K(7, 0, 2, 0, 4,
466                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
467                             PRIM_INT__k_s_p_s, PRIM_INT__k_s_s_s, NULL, PRIM_INT__i_s_p_s, NULL, PRIM_INT__k_s_d_s);
468 
469 
470                     ostei_general_vrr_K(5, 0, 1, 0, 5,
471                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
472                             PRIM_INT__h_s_s_s, NULL, NULL, PRIM_INT__g_s_s_s, NULL, PRIM_INT__h_s_p_s);
473 
474 
475                     ostei_general_vrr_K(6, 0, 2, 0, 4,
476                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
477                             PRIM_INT__i_s_p_s, PRIM_INT__i_s_s_s, NULL, PRIM_INT__h_s_p_s, NULL, PRIM_INT__i_s_d_s);
478 
479 
480                     ostei_general_vrr_K(7, 0, 3, 0, 3,
481                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
482                             PRIM_INT__k_s_d_s, PRIM_INT__k_s_p_s, NULL, PRIM_INT__i_s_d_s, NULL, PRIM_INT__k_s_f_s);
483 
484 
485                     VRR_K_g_s_p_s(
486                             PRIM_INT__g_s_p_s,
487                             PRIM_INT__g_s_s_s,
488                             PRIM_INT__f_s_s_s,
489                             Q_PA,
490                             aoq_PQ,
491                             one_over_2pq,
492                             5);
493 
494 
495                     ostei_general_vrr_K(5, 0, 2, 0, 4,
496                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
497                             PRIM_INT__h_s_p_s, PRIM_INT__h_s_s_s, NULL, PRIM_INT__g_s_p_s, NULL, PRIM_INT__h_s_d_s);
498 
499 
500                     ostei_general_vrr_K(6, 0, 3, 0, 3,
501                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
502                             PRIM_INT__i_s_d_s, PRIM_INT__i_s_p_s, NULL, PRIM_INT__h_s_d_s, NULL, PRIM_INT__i_s_f_s);
503 
504 
505                     ostei_general_vrr_K(7, 0, 4, 0, 2,
506                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
507                             PRIM_INT__k_s_f_s, PRIM_INT__k_s_d_s, NULL, PRIM_INT__i_s_f_s, NULL, PRIM_INT__k_s_g_s);
508 
509 
510                     VRR_K_f_s_p_s(
511                             PRIM_INT__f_s_p_s,
512                             PRIM_INT__f_s_s_s,
513                             PRIM_INT__d_s_s_s,
514                             Q_PA,
515                             aoq_PQ,
516                             one_over_2pq,
517                             5);
518 
519 
520                     ostei_general_vrr_K(4, 0, 2, 0, 4,
521                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
522                             PRIM_INT__g_s_p_s, PRIM_INT__g_s_s_s, NULL, PRIM_INT__f_s_p_s, NULL, PRIM_INT__g_s_d_s);
523 
524 
525                     ostei_general_vrr_K(5, 0, 3, 0, 3,
526                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
527                             PRIM_INT__h_s_d_s, PRIM_INT__h_s_p_s, NULL, PRIM_INT__g_s_d_s, NULL, PRIM_INT__h_s_f_s);
528 
529 
530                     ostei_general_vrr_K(6, 0, 4, 0, 2,
531                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
532                             PRIM_INT__i_s_f_s, PRIM_INT__i_s_d_s, NULL, PRIM_INT__h_s_f_s, NULL, PRIM_INT__i_s_g_s);
533 
534 
535                     ostei_general_vrr_K(7, 0, 5, 0, 1,
536                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
537                             PRIM_INT__k_s_g_s, PRIM_INT__k_s_f_s, NULL, PRIM_INT__i_s_g_s, NULL, PRIM_INT__k_s_h_s);
538 
539 
540                     ostei_general_vrr1_I(8, 9,
541                             one_over_2p, a_over_p, aop_PQ, P_PA,
542                             PRIM_INT__k_s_s_s, PRIM_INT__i_s_s_s, PRIM_INT__l_s_s_s);
543 
544 
545                     ostei_general_vrr_K(8, 0, 1, 0, 5,
546                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
547                             PRIM_INT__l_s_s_s, NULL, NULL, PRIM_INT__k_s_s_s, NULL, PRIM_INT__l_s_p_s);
548 
549 
550                     ostei_general_vrr_K(8, 0, 2, 0, 4,
551                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
552                             PRIM_INT__l_s_p_s, PRIM_INT__l_s_s_s, NULL, PRIM_INT__k_s_p_s, NULL, PRIM_INT__l_s_d_s);
553 
554 
555                     ostei_general_vrr_K(8, 0, 3, 0, 3,
556                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
557                             PRIM_INT__l_s_d_s, PRIM_INT__l_s_p_s, NULL, PRIM_INT__k_s_d_s, NULL, PRIM_INT__l_s_f_s);
558 
559 
560                     ostei_general_vrr_K(8, 0, 4, 0, 2,
561                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
562                             PRIM_INT__l_s_f_s, PRIM_INT__l_s_d_s, NULL, PRIM_INT__k_s_f_s, NULL, PRIM_INT__l_s_g_s);
563 
564 
565                     ostei_general_vrr_K(8, 0, 5, 0, 1,
566                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
567                             PRIM_INT__l_s_g_s, PRIM_INT__l_s_f_s, NULL, PRIM_INT__k_s_g_s, NULL, PRIM_INT__l_s_h_s);
568 
569 
570                     ostei_general_vrr1_I(9, 8,
571                             one_over_2p, a_over_p, aop_PQ, P_PA,
572                             PRIM_INT__l_s_s_s, PRIM_INT__k_s_s_s, PRIM_INT__m_s_s_s);
573 
574 
575                     ostei_general_vrr_K(9, 0, 1, 0, 5,
576                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
577                             PRIM_INT__m_s_s_s, NULL, NULL, PRIM_INT__l_s_s_s, NULL, PRIM_INT__m_s_p_s);
578 
579 
580                     ostei_general_vrr_K(9, 0, 2, 0, 4,
581                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
582                             PRIM_INT__m_s_p_s, PRIM_INT__m_s_s_s, NULL, PRIM_INT__l_s_p_s, NULL, PRIM_INT__m_s_d_s);
583 
584 
585                     ostei_general_vrr_K(9, 0, 3, 0, 3,
586                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
587                             PRIM_INT__m_s_d_s, PRIM_INT__m_s_p_s, NULL, PRIM_INT__l_s_d_s, NULL, PRIM_INT__m_s_f_s);
588 
589 
590                     ostei_general_vrr_K(9, 0, 4, 0, 2,
591                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
592                             PRIM_INT__m_s_f_s, PRIM_INT__m_s_d_s, NULL, PRIM_INT__l_s_f_s, NULL, PRIM_INT__m_s_g_s);
593 
594 
595                     ostei_general_vrr_K(9, 0, 5, 0, 1,
596                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
597                             PRIM_INT__m_s_g_s, PRIM_INT__m_s_f_s, NULL, PRIM_INT__l_s_g_s, NULL, PRIM_INT__m_s_h_s);
598 
599 
600                     ostei_general_vrr1_I(10, 7,
601                             one_over_2p, a_over_p, aop_PQ, P_PA,
602                             PRIM_INT__m_s_s_s, PRIM_INT__l_s_s_s, PRIM_INT__n_s_s_s);
603 
604 
605                     ostei_general_vrr_K(10, 0, 1, 0, 5,
606                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
607                             PRIM_INT__n_s_s_s, NULL, NULL, PRIM_INT__m_s_s_s, NULL, PRIM_INT__n_s_p_s);
608 
609 
610                     ostei_general_vrr_K(10, 0, 2, 0, 4,
611                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
612                             PRIM_INT__n_s_p_s, PRIM_INT__n_s_s_s, NULL, PRIM_INT__m_s_p_s, NULL, PRIM_INT__n_s_d_s);
613 
614 
615                     ostei_general_vrr_K(10, 0, 3, 0, 3,
616                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
617                             PRIM_INT__n_s_d_s, PRIM_INT__n_s_p_s, NULL, PRIM_INT__m_s_d_s, NULL, PRIM_INT__n_s_f_s);
618 
619 
620                     ostei_general_vrr_K(10, 0, 4, 0, 2,
621                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
622                             PRIM_INT__n_s_f_s, PRIM_INT__n_s_d_s, NULL, PRIM_INT__m_s_f_s, NULL, PRIM_INT__n_s_g_s);
623 
624 
625                     ostei_general_vrr_K(10, 0, 5, 0, 1,
626                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
627                             PRIM_INT__n_s_g_s, PRIM_INT__n_s_f_s, NULL, PRIM_INT__m_s_g_s, NULL, PRIM_INT__n_s_h_s);
628 
629 
630                     ostei_general_vrr1_I(11, 6,
631                             one_over_2p, a_over_p, aop_PQ, P_PA,
632                             PRIM_INT__n_s_s_s, PRIM_INT__m_s_s_s, PRIM_INT__o_s_s_s);
633 
634 
635                     ostei_general_vrr_K(11, 0, 1, 0, 5,
636                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
637                             PRIM_INT__o_s_s_s, NULL, NULL, PRIM_INT__n_s_s_s, NULL, PRIM_INT__o_s_p_s);
638 
639 
640                     ostei_general_vrr_K(11, 0, 2, 0, 4,
641                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
642                             PRIM_INT__o_s_p_s, PRIM_INT__o_s_s_s, NULL, PRIM_INT__n_s_p_s, NULL, PRIM_INT__o_s_d_s);
643 
644 
645                     ostei_general_vrr_K(11, 0, 3, 0, 3,
646                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
647                             PRIM_INT__o_s_d_s, PRIM_INT__o_s_p_s, NULL, PRIM_INT__n_s_d_s, NULL, PRIM_INT__o_s_f_s);
648 
649 
650                     ostei_general_vrr_K(11, 0, 4, 0, 2,
651                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
652                             PRIM_INT__o_s_f_s, PRIM_INT__o_s_d_s, NULL, PRIM_INT__n_s_f_s, NULL, PRIM_INT__o_s_g_s);
653 
654 
655                     ostei_general_vrr_K(11, 0, 5, 0, 1,
656                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
657                             PRIM_INT__o_s_g_s, PRIM_INT__o_s_f_s, NULL, PRIM_INT__n_s_g_s, NULL, PRIM_INT__o_s_h_s);
658 
659 
660 
661 
662                     ////////////////////////////////////
663                     // Accumulate contracted integrals
664                     ////////////////////////////////////
665                     if(lastoffset == 0)
666                     {
667                         contract_all(540, PRIM_INT__k_s_g_s, PRIM_PTR_INT__k_s_g_s);
668                         contract_all(756, PRIM_INT__k_s_h_s, PRIM_PTR_INT__k_s_h_s);
669                         contract_all(675, PRIM_INT__l_s_g_s, PRIM_PTR_INT__l_s_g_s);
670                         contract_all(945, PRIM_INT__l_s_h_s, PRIM_PTR_INT__l_s_h_s);
671                         contract_all(825, PRIM_INT__m_s_g_s, PRIM_PTR_INT__m_s_g_s);
672                         contract_all(1155, PRIM_INT__m_s_h_s, PRIM_PTR_INT__m_s_h_s);
673                         contract_all(990, PRIM_INT__n_s_g_s, PRIM_PTR_INT__n_s_g_s);
674                         contract_all(1386, PRIM_INT__n_s_h_s, PRIM_PTR_INT__n_s_h_s);
675                         contract_all(1170, PRIM_INT__o_s_g_s, PRIM_PTR_INT__o_s_g_s);
676                         contract_all(1638, PRIM_INT__o_s_h_s, PRIM_PTR_INT__o_s_h_s);
677                     }
678                     else
679                     {
680                         contract(540, shelloffsets, PRIM_INT__k_s_g_s, PRIM_PTR_INT__k_s_g_s);
681                         contract(756, shelloffsets, PRIM_INT__k_s_h_s, PRIM_PTR_INT__k_s_h_s);
682                         contract(675, shelloffsets, PRIM_INT__l_s_g_s, PRIM_PTR_INT__l_s_g_s);
683                         contract(945, shelloffsets, PRIM_INT__l_s_h_s, PRIM_PTR_INT__l_s_h_s);
684                         contract(825, shelloffsets, PRIM_INT__m_s_g_s, PRIM_PTR_INT__m_s_g_s);
685                         contract(1155, shelloffsets, PRIM_INT__m_s_h_s, PRIM_PTR_INT__m_s_h_s);
686                         contract(990, shelloffsets, PRIM_INT__n_s_g_s, PRIM_PTR_INT__n_s_g_s);
687                         contract(1386, shelloffsets, PRIM_INT__n_s_h_s, PRIM_PTR_INT__n_s_h_s);
688                         contract(1170, shelloffsets, PRIM_INT__o_s_g_s, PRIM_PTR_INT__o_s_g_s);
689                         contract(1638, shelloffsets, PRIM_INT__o_s_h_s, PRIM_PTR_INT__o_s_h_s);
690                         PRIM_PTR_INT__k_s_g_s += lastoffset*540;
691                         PRIM_PTR_INT__k_s_h_s += lastoffset*756;
692                         PRIM_PTR_INT__l_s_g_s += lastoffset*675;
693                         PRIM_PTR_INT__l_s_h_s += lastoffset*945;
694                         PRIM_PTR_INT__m_s_g_s += lastoffset*825;
695                         PRIM_PTR_INT__m_s_h_s += lastoffset*1155;
696                         PRIM_PTR_INT__n_s_g_s += lastoffset*990;
697                         PRIM_PTR_INT__n_s_h_s += lastoffset*1386;
698                         PRIM_PTR_INT__o_s_g_s += lastoffset*1170;
699                         PRIM_PTR_INT__o_s_h_s += lastoffset*1638;
700                     }
701 
702                 }  // close loop over j
703             }  // close loop over i
704 
705             //Advance to the next batch
706             jstart = SIMINT_SIMD_ROUND(jend);
707 
708             //////////////////////////////////////////////
709             // Contracted integrals: Horizontal recurrance
710             //////////////////////////////////////////////
711 
712 
713             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
714 
715 
716             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
717             {
718                 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
719 
720                 // set up HRR pointers
721                 double const * restrict HRR_INT__k_s_g_s = INT__k_s_g_s + abcd * 540;
722                 double const * restrict HRR_INT__k_s_h_s = INT__k_s_h_s + abcd * 756;
723                 double const * restrict HRR_INT__l_s_g_s = INT__l_s_g_s + abcd * 675;
724                 double const * restrict HRR_INT__l_s_h_s = INT__l_s_h_s + abcd * 945;
725                 double const * restrict HRR_INT__m_s_g_s = INT__m_s_g_s + abcd * 825;
726                 double const * restrict HRR_INT__m_s_h_s = INT__m_s_h_s + abcd * 1155;
727                 double const * restrict HRR_INT__n_s_g_s = INT__n_s_g_s + abcd * 990;
728                 double const * restrict HRR_INT__n_s_h_s = INT__n_s_h_s + abcd * 1386;
729                 double const * restrict HRR_INT__o_s_g_s = INT__o_s_g_s + abcd * 1170;
730                 double const * restrict HRR_INT__o_s_h_s = INT__o_s_h_s + abcd * 1638;
731                 double * restrict HRR_INT__k_g_g_p = INT__k_g_g_p + real_abcd * 24300;
732 
733                 // form INT__k_p_g_s
734                 ostei_general_hrr_J(7, 1, 4, 0, hAB, HRR_INT__l_s_g_s, HRR_INT__k_s_g_s, HRR_INT__k_p_g_s);
735 
736                 // form INT__k_p_h_s
737                 ostei_general_hrr_J(7, 1, 5, 0, hAB, HRR_INT__l_s_h_s, HRR_INT__k_s_h_s, HRR_INT__k_p_h_s);
738 
739                 // form INT__l_p_g_s
740                 ostei_general_hrr_J(8, 1, 4, 0, hAB, HRR_INT__m_s_g_s, HRR_INT__l_s_g_s, HRR_INT__l_p_g_s);
741 
742                 // form INT__l_p_h_s
743                 ostei_general_hrr_J(8, 1, 5, 0, hAB, HRR_INT__m_s_h_s, HRR_INT__l_s_h_s, HRR_INT__l_p_h_s);
744 
745                 // form INT__m_p_g_s
746                 ostei_general_hrr_J(9, 1, 4, 0, hAB, HRR_INT__n_s_g_s, HRR_INT__m_s_g_s, HRR_INT__m_p_g_s);
747 
748                 // form INT__m_p_h_s
749                 ostei_general_hrr_J(9, 1, 5, 0, hAB, HRR_INT__n_s_h_s, HRR_INT__m_s_h_s, HRR_INT__m_p_h_s);
750 
751                 // form INT__n_p_g_s
752                 ostei_general_hrr_J(10, 1, 4, 0, hAB, HRR_INT__o_s_g_s, HRR_INT__n_s_g_s, HRR_INT__n_p_g_s);
753 
754                 // form INT__n_p_h_s
755                 ostei_general_hrr_J(10, 1, 5, 0, hAB, HRR_INT__o_s_h_s, HRR_INT__n_s_h_s, HRR_INT__n_p_h_s);
756 
757                 // form INT__k_d_g_s
758                 ostei_general_hrr_J(7, 2, 4, 0, hAB, HRR_INT__l_p_g_s, HRR_INT__k_p_g_s, HRR_INT__k_d_g_s);
759 
760                 // form INT__k_d_h_s
761                 ostei_general_hrr_J(7, 2, 5, 0, hAB, HRR_INT__l_p_h_s, HRR_INT__k_p_h_s, HRR_INT__k_d_h_s);
762 
763                 // form INT__l_d_g_s
764                 ostei_general_hrr_J(8, 2, 4, 0, hAB, HRR_INT__m_p_g_s, HRR_INT__l_p_g_s, HRR_INT__l_d_g_s);
765 
766                 // form INT__l_d_h_s
767                 ostei_general_hrr_J(8, 2, 5, 0, hAB, HRR_INT__m_p_h_s, HRR_INT__l_p_h_s, HRR_INT__l_d_h_s);
768 
769                 // form INT__m_d_g_s
770                 ostei_general_hrr_J(9, 2, 4, 0, hAB, HRR_INT__n_p_g_s, HRR_INT__m_p_g_s, HRR_INT__m_d_g_s);
771 
772                 // form INT__m_d_h_s
773                 ostei_general_hrr_J(9, 2, 5, 0, hAB, HRR_INT__n_p_h_s, HRR_INT__m_p_h_s, HRR_INT__m_d_h_s);
774 
775                 // form INT__k_f_g_s
776                 ostei_general_hrr_J(7, 3, 4, 0, hAB, HRR_INT__l_d_g_s, HRR_INT__k_d_g_s, HRR_INT__k_f_g_s);
777 
778                 // form INT__k_f_h_s
779                 ostei_general_hrr_J(7, 3, 5, 0, hAB, HRR_INT__l_d_h_s, HRR_INT__k_d_h_s, HRR_INT__k_f_h_s);
780 
781                 // form INT__l_f_g_s
782                 ostei_general_hrr_J(8, 3, 4, 0, hAB, HRR_INT__m_d_g_s, HRR_INT__l_d_g_s, HRR_INT__l_f_g_s);
783 
784                 // form INT__l_f_h_s
785                 ostei_general_hrr_J(8, 3, 5, 0, hAB, HRR_INT__m_d_h_s, HRR_INT__l_d_h_s, HRR_INT__l_f_h_s);
786 
787                 // form INT__k_g_g_s
788                 ostei_general_hrr_J(7, 4, 4, 0, hAB, HRR_INT__l_f_g_s, HRR_INT__k_f_g_s, HRR_INT__k_g_g_s);
789 
790                 // form INT__k_g_h_s
791                 ostei_general_hrr_J(7, 4, 5, 0, hAB, HRR_INT__l_f_h_s, HRR_INT__k_f_h_s, HRR_INT__k_g_h_s);
792 
793                 // form INT__k_g_g_p
794                 HRR_L_g_p(
795                     HRR_INT__k_g_g_p,
796                     HRR_INT__k_g_g_s,
797                     HRR_INT__k_g_h_s,
798                     hCD, 540);
799 
800 
801             }  // close HRR loop
802 
803 
804         }   // close loop cdbatch
805 
806         istart = iend;
807     }  // close loop over ab
808 
809     return P.nshell12_clip * Q.nshell12_clip;
810 }
811 
ostei_g_k_g_p(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__g_k_g_p)812 int ostei_g_k_g_p(struct simint_multi_shellpair const P,
813                   struct simint_multi_shellpair const Q,
814                   double screen_tol,
815                   double * const restrict work,
816                   double * const restrict INT__g_k_g_p)
817 {
818     double P_AB[3*P.nshell12];
819     struct simint_multi_shellpair P_tmp = P;
820     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
821     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
822     P_tmp.AB_x = P_AB;
823     P_tmp.AB_y = P_AB + P.nshell12;
824     P_tmp.AB_z = P_AB + 2*P.nshell12;
825 
826     for(int i = 0; i < P.nshell12; i++)
827     {
828         P_tmp.AB_x[i] = -P.AB_x[i];
829         P_tmp.AB_y[i] = -P.AB_y[i];
830         P_tmp.AB_z[i] = -P.AB_z[i];
831     }
832 
833     int ret = ostei_k_g_g_p(P_tmp, Q, screen_tol, work, INT__g_k_g_p);
834     double buffer[24300] SIMINT_ALIGN_ARRAY_DBL;
835 
836     for(int q = 0; q < ret; q++)
837     {
838         int idx = 0;
839         for(int a = 0; a < 15; ++a)
840         for(int b = 0; b < 36; ++b)
841         for(int c = 0; c < 15; ++c)
842         for(int d = 0; d < 3; ++d)
843             buffer[idx++] = INT__g_k_g_p[q*24300+b*675+a*45+c*3+d];
844 
845         memcpy(INT__g_k_g_p+q*24300, buffer, 24300*sizeof(double));
846     }
847 
848     return ret;
849 }
850 
ostei_k_g_p_g(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__k_g_p_g)851 int ostei_k_g_p_g(struct simint_multi_shellpair const P,
852                   struct simint_multi_shellpair const Q,
853                   double screen_tol,
854                   double * const restrict work,
855                   double * const restrict INT__k_g_p_g)
856 {
857     double Q_AB[3*Q.nshell12];
858     struct simint_multi_shellpair Q_tmp = Q;
859     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
860     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
861     Q_tmp.AB_x = Q_AB;
862     Q_tmp.AB_y = Q_AB + Q.nshell12;
863     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
864 
865     for(int i = 0; i < Q.nshell12; i++)
866     {
867         Q_tmp.AB_x[i] = -Q.AB_x[i];
868         Q_tmp.AB_y[i] = -Q.AB_y[i];
869         Q_tmp.AB_z[i] = -Q.AB_z[i];
870     }
871 
872     int ret = ostei_k_g_g_p(P, Q_tmp, screen_tol, work, INT__k_g_p_g);
873     double buffer[24300] SIMINT_ALIGN_ARRAY_DBL;
874 
875     for(int q = 0; q < ret; q++)
876     {
877         int idx = 0;
878         for(int a = 0; a < 36; ++a)
879         for(int b = 0; b < 15; ++b)
880         for(int c = 0; c < 3; ++c)
881         for(int d = 0; d < 15; ++d)
882             buffer[idx++] = INT__k_g_p_g[q*24300+a*675+b*45+d*3+c];
883 
884         memcpy(INT__k_g_p_g+q*24300, buffer, 24300*sizeof(double));
885     }
886 
887     return ret;
888 }
889 
ostei_g_k_p_g(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__g_k_p_g)890 int ostei_g_k_p_g(struct simint_multi_shellpair const P,
891                   struct simint_multi_shellpair const Q,
892                   double screen_tol,
893                   double * const restrict work,
894                   double * const restrict INT__g_k_p_g)
895 {
896     double P_AB[3*P.nshell12];
897     struct simint_multi_shellpair P_tmp = P;
898     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
899     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
900     P_tmp.AB_x = P_AB;
901     P_tmp.AB_y = P_AB + P.nshell12;
902     P_tmp.AB_z = P_AB + 2*P.nshell12;
903 
904     for(int i = 0; i < P.nshell12; i++)
905     {
906         P_tmp.AB_x[i] = -P.AB_x[i];
907         P_tmp.AB_y[i] = -P.AB_y[i];
908         P_tmp.AB_z[i] = -P.AB_z[i];
909     }
910 
911     double Q_AB[3*Q.nshell12];
912     struct simint_multi_shellpair Q_tmp = Q;
913     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
914     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
915     Q_tmp.AB_x = Q_AB;
916     Q_tmp.AB_y = Q_AB + Q.nshell12;
917     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
918 
919     for(int i = 0; i < Q.nshell12; i++)
920     {
921         Q_tmp.AB_x[i] = -Q.AB_x[i];
922         Q_tmp.AB_y[i] = -Q.AB_y[i];
923         Q_tmp.AB_z[i] = -Q.AB_z[i];
924     }
925 
926     int ret = ostei_k_g_g_p(P_tmp, Q_tmp, screen_tol, work, INT__g_k_p_g);
927     double buffer[24300] SIMINT_ALIGN_ARRAY_DBL;
928 
929     for(int q = 0; q < ret; q++)
930     {
931         int idx = 0;
932         for(int a = 0; a < 15; ++a)
933         for(int b = 0; b < 36; ++b)
934         for(int c = 0; c < 3; ++c)
935         for(int d = 0; d < 15; ++d)
936             buffer[idx++] = INT__g_k_p_g[q*24300+b*675+a*45+d*3+c];
937 
938         memcpy(INT__g_k_p_g+q*24300, buffer, 24300*sizeof(double));
939     }
940 
941     return ret;
942 }
943 
944