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_i_h_k_p(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__i_h_k_p)8 int ostei_i_h_k_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__i_h_k_p)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__i_h_k_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__i_s_k_s = work + (SIMINT_NSHELL_SIMD * 0);
30     double * const INT__i_s_l_s = work + (SIMINT_NSHELL_SIMD * 1008);
31     double * const INT__k_s_k_s = work + (SIMINT_NSHELL_SIMD * 2268);
32     double * const INT__k_s_l_s = work + (SIMINT_NSHELL_SIMD * 3564);
33     double * const INT__l_s_k_s = work + (SIMINT_NSHELL_SIMD * 5184);
34     double * const INT__l_s_l_s = work + (SIMINT_NSHELL_SIMD * 6804);
35     double * const INT__m_s_k_s = work + (SIMINT_NSHELL_SIMD * 8829);
36     double * const INT__m_s_l_s = work + (SIMINT_NSHELL_SIMD * 10809);
37     double * const INT__n_s_k_s = work + (SIMINT_NSHELL_SIMD * 13284);
38     double * const INT__n_s_l_s = work + (SIMINT_NSHELL_SIMD * 15660);
39     double * const INT__o_s_k_s = work + (SIMINT_NSHELL_SIMD * 18630);
40     double * const INT__o_s_l_s = work + (SIMINT_NSHELL_SIMD * 21438);
41     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*24948);
42     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
43     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_p_s = primwork + 20;
44     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_d_s = primwork + 44;
45     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_s = primwork + 86;
46     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_p_s = primwork + 143;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_d_s = primwork + 215;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_f_s = primwork + 341;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 521;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_p_s = primwork + 629;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_d_s = primwork + 773;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_f_s = primwork + 1025;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_g_s = primwork + 1385;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_s = primwork + 1835;
55     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_p_s = primwork + 2005;
56     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_d_s = primwork + 2245;
57     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_f_s = primwork + 2665;
58     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_g_s = primwork + 3265;
59     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_h_s = primwork + 4015;
60     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_s_s = primwork + 4855;
61     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_p_s = primwork + 5095;
62     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_d_s = primwork + 5455;
63     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_f_s = primwork + 6085;
64     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_g_s = primwork + 6985;
65     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_h_s = primwork + 8110;
66     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_i_s = primwork + 9370;
67     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_s_s = primwork + 10630;
68     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_p_s = primwork + 10945;
69     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_d_s = primwork + 11449;
70     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_f_s = primwork + 12331;
71     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_g_s = primwork + 13591;
72     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_h_s = primwork + 15166;
73     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_i_s = primwork + 16930;
74     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_k_s = primwork + 18694;
75     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_s_s = primwork + 20206;
76     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_p_s = primwork + 20598;
77     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_d_s = primwork + 21270;
78     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_f_s = primwork + 22446;
79     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_g_s = primwork + 24126;
80     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_h_s = primwork + 26226;
81     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_i_s = primwork + 28578;
82     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_k_s = primwork + 30930;
83     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_l_s = primwork + 32946;
84     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_s_s = primwork + 34206;
85     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_p_s = primwork + 34674;
86     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_d_s = primwork + 35538;
87     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_f_s = primwork + 37050;
88     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_g_s = primwork + 39210;
89     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_h_s = primwork + 41910;
90     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_i_s = primwork + 44934;
91     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_k_s = primwork + 47958;
92     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_l_s = primwork + 50550;
93     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_s_s = primwork + 52170;
94     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_p_s = primwork + 52710;
95     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_d_s = primwork + 53790;
96     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_f_s = primwork + 55680;
97     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_g_s = primwork + 58380;
98     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_h_s = primwork + 61755;
99     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_i_s = primwork + 65535;
100     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_k_s = primwork + 69315;
101     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_l_s = primwork + 72555;
102     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_s_s = primwork + 74580;
103     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_p_s = primwork + 75185;
104     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_d_s = primwork + 76505;
105     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_f_s = primwork + 78815;
106     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_g_s = primwork + 82115;
107     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_h_s = primwork + 86240;
108     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_i_s = primwork + 90860;
109     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_k_s = primwork + 95480;
110     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_l_s = primwork + 99440;
111     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_s_s = primwork + 101915;
112     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_p_s = primwork + 102575;
113     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_d_s = primwork + 104159;
114     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_f_s = primwork + 106931;
115     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_g_s = primwork + 110891;
116     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_h_s = primwork + 115841;
117     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_i_s = primwork + 121385;
118     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_k_s = primwork + 126929;
119     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_l_s = primwork + 131681;
120     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_s_s = primwork + 134651;
121     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_p_s = primwork + 135353;
122     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_d_s = primwork + 137225;
123     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_f_s = primwork + 140501;
124     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_g_s = primwork + 145181;
125     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_h_s = primwork + 151031;
126     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_i_s = primwork + 157583;
127     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_k_s = primwork + 164135;
128     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_l_s = primwork + 169751;
129     double * const hrrwork = (double *)(primwork + 173261);
130     double * const HRR_INT__i_p_k_s = hrrwork + 0;
131     double * const HRR_INT__i_p_l_s = hrrwork + 3024;
132     double * const HRR_INT__i_d_k_s = hrrwork + 6804;
133     double * const HRR_INT__i_d_l_s = hrrwork + 12852;
134     double * const HRR_INT__i_f_k_s = hrrwork + 20412;
135     double * const HRR_INT__i_f_l_s = hrrwork + 30492;
136     double * const HRR_INT__i_g_k_s = hrrwork + 43092;
137     double * const HRR_INT__i_g_l_s = hrrwork + 58212;
138     double * const HRR_INT__i_h_k_s = hrrwork + 77112;
139     double * const HRR_INT__i_h_l_s = hrrwork + 98280;
140     double * const HRR_INT__k_p_k_s = hrrwork + 124740;
141     double * const HRR_INT__k_p_l_s = hrrwork + 128628;
142     double * const HRR_INT__k_d_k_s = hrrwork + 133488;
143     double * const HRR_INT__k_d_l_s = hrrwork + 141264;
144     double * const HRR_INT__k_f_k_s = hrrwork + 150984;
145     double * const HRR_INT__k_f_l_s = hrrwork + 163944;
146     double * const HRR_INT__k_g_k_s = hrrwork + 180144;
147     double * const HRR_INT__k_g_l_s = hrrwork + 199584;
148     double * const HRR_INT__l_p_k_s = hrrwork + 223884;
149     double * const HRR_INT__l_p_l_s = hrrwork + 228744;
150     double * const HRR_INT__l_d_k_s = hrrwork + 234819;
151     double * const HRR_INT__l_d_l_s = hrrwork + 244539;
152     double * const HRR_INT__l_f_k_s = hrrwork + 256689;
153     double * const HRR_INT__l_f_l_s = hrrwork + 272889;
154     double * const HRR_INT__m_p_k_s = hrrwork + 293139;
155     double * const HRR_INT__m_p_l_s = hrrwork + 299079;
156     double * const HRR_INT__m_d_k_s = hrrwork + 306504;
157     double * const HRR_INT__m_d_l_s = hrrwork + 318384;
158     double * const HRR_INT__n_p_k_s = hrrwork + 333234;
159     double * const HRR_INT__n_p_l_s = hrrwork + 340362;
160 
161 
162     // Create constants
163     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
164     const SIMINT_DBLTYPE const_10 = SIMINT_DBLSET1(10);
165     const SIMINT_DBLTYPE const_11 = SIMINT_DBLSET1(11);
166     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
167     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
168     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
169     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
170     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
171     const SIMINT_DBLTYPE const_7 = SIMINT_DBLSET1(7);
172     const SIMINT_DBLTYPE const_8 = SIMINT_DBLSET1(8);
173     const SIMINT_DBLTYPE const_9 = SIMINT_DBLSET1(9);
174     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
175 
176 
177     ////////////////////////////////////////
178     // Loop over shells and primitives
179     ////////////////////////////////////////
180 
181     real_abcd = 0;
182     istart = 0;
183     for(ab = 0; ab < P.nshell12_clip; ++ab)
184     {
185         const int iend = istart + P.nprim12[ab];
186 
187         cd = 0;
188         jstart = 0;
189 
190         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
191         {
192             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
193             int jend = jstart;
194             for(i = 0; i < nshellbatch; i++)
195                 jend += Q.nprim12[cd+i];
196 
197             // Clear the beginning of the workspace (where we are accumulating integrals)
198             memset(work, 0, SIMINT_NSHELL_SIMD * 24948 * sizeof(double));
199             abcd = 0;
200 
201 
202             for(i = istart; i < iend; ++i)
203             {
204                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
205 
206                 if(check_screen)
207                 {
208                     // Skip this whole thing if always insignificant
209                     if((P.screen[i] * Q.screen_max) < screen_tol)
210                         continue;
211                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
212                 }
213 
214                 icd = 0;
215                 iprimcd = 0;
216                 nprim_icd = Q.nprim12[cd];
217                 double * restrict PRIM_PTR_INT__i_s_k_s = INT__i_s_k_s + abcd * 1008;
218                 double * restrict PRIM_PTR_INT__i_s_l_s = INT__i_s_l_s + abcd * 1260;
219                 double * restrict PRIM_PTR_INT__k_s_k_s = INT__k_s_k_s + abcd * 1296;
220                 double * restrict PRIM_PTR_INT__k_s_l_s = INT__k_s_l_s + abcd * 1620;
221                 double * restrict PRIM_PTR_INT__l_s_k_s = INT__l_s_k_s + abcd * 1620;
222                 double * restrict PRIM_PTR_INT__l_s_l_s = INT__l_s_l_s + abcd * 2025;
223                 double * restrict PRIM_PTR_INT__m_s_k_s = INT__m_s_k_s + abcd * 1980;
224                 double * restrict PRIM_PTR_INT__m_s_l_s = INT__m_s_l_s + abcd * 2475;
225                 double * restrict PRIM_PTR_INT__n_s_k_s = INT__n_s_k_s + abcd * 2376;
226                 double * restrict PRIM_PTR_INT__n_s_l_s = INT__n_s_l_s + abcd * 2970;
227                 double * restrict PRIM_PTR_INT__o_s_k_s = INT__o_s_k_s + abcd * 2808;
228                 double * restrict PRIM_PTR_INT__o_s_l_s = INT__o_s_l_s + abcd * 3510;
229 
230 
231 
232                 // Load these one per loop over i
233                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
234                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
235                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
236 
237                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
238 
239                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
240                 {
241                     // calculate the shell offsets
242                     // these are the offset from the shell pointed to by cd
243                     // for each element
244                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
245                     int lastoffset = 0;
246                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
247 
248                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
249                     {
250                         // Handle if the first element of the vector is a new shell
251                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
252                         {
253                             nprim_icd += Q.nprim12[cd + (++icd)];
254                             PRIM_PTR_INT__i_s_k_s += 1008;
255                             PRIM_PTR_INT__i_s_l_s += 1260;
256                             PRIM_PTR_INT__k_s_k_s += 1296;
257                             PRIM_PTR_INT__k_s_l_s += 1620;
258                             PRIM_PTR_INT__l_s_k_s += 1620;
259                             PRIM_PTR_INT__l_s_l_s += 2025;
260                             PRIM_PTR_INT__m_s_k_s += 1980;
261                             PRIM_PTR_INT__m_s_l_s += 2475;
262                             PRIM_PTR_INT__n_s_k_s += 2376;
263                             PRIM_PTR_INT__n_s_l_s += 2970;
264                             PRIM_PTR_INT__o_s_k_s += 2808;
265                             PRIM_PTR_INT__o_s_l_s += 3510;
266                         }
267                         iprimcd++;
268                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
269                         {
270                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
271                             {
272                                 shelloffsets[n] = shelloffsets[n-1] + 1;
273                                 lastoffset++;
274                                 nprim_icd += Q.nprim12[cd + (++icd)];
275                             }
276                             else
277                                 shelloffsets[n] = shelloffsets[n-1];
278                             iprimcd++;
279                         }
280                     }
281                     else
282                         iprimcd += SIMINT_SIMD_LEN;
283 
284                     // Do we have to compute this vector (or has it been screened out)?
285                     // (not_screened != 0 means we have to do this vector)
286                     if(check_screen)
287                     {
288                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
289                         if(vmax < screen_tol)
290                         {
291                             PRIM_PTR_INT__i_s_k_s += lastoffset*1008;
292                             PRIM_PTR_INT__i_s_l_s += lastoffset*1260;
293                             PRIM_PTR_INT__k_s_k_s += lastoffset*1296;
294                             PRIM_PTR_INT__k_s_l_s += lastoffset*1620;
295                             PRIM_PTR_INT__l_s_k_s += lastoffset*1620;
296                             PRIM_PTR_INT__l_s_l_s += lastoffset*2025;
297                             PRIM_PTR_INT__m_s_k_s += lastoffset*1980;
298                             PRIM_PTR_INT__m_s_l_s += lastoffset*2475;
299                             PRIM_PTR_INT__n_s_k_s += lastoffset*2376;
300                             PRIM_PTR_INT__n_s_l_s += lastoffset*2970;
301                             PRIM_PTR_INT__o_s_k_s += lastoffset*2808;
302                             PRIM_PTR_INT__o_s_l_s += lastoffset*3510;
303                             continue;
304                         }
305                     }
306 
307                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
308                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
309                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
310                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
311 
312 
313                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
314                     SIMINT_DBLTYPE PQ[3];
315                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
316                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
317                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
318                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
319                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
320                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
321 
322                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
323                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
324                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
325                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
326                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
327                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
328                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
329 
330                     // NOTE: Minus sign!
331                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
332                     SIMINT_DBLTYPE aop_PQ[3];
333                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
334                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
335                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
336 
337                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
338                     SIMINT_DBLTYPE aoq_PQ[3];
339                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
340                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
341                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
342                     // Put a minus sign here so we don't have to in RR routines
343                     a_over_q = SIMINT_NEG(a_over_q);
344 
345 
346                     //////////////////////////////////////////////
347                     // Fjt function section
348                     // Maximum v value: 19
349                     //////////////////////////////////////////////
350                     // The parameter to the Fjt function
351                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
352 
353 
354                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
355 
356 
357                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 19);
358                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
359                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
360                     for(n = 0; n <= 19; n++)
361                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
362 
363                     //////////////////////////////////////////////
364                     // Primitive integrals: Vertical recurrance
365                     //////////////////////////////////////////////
366 
367                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
368                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
369                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
370                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
371                     const SIMINT_DBLTYPE vrr_const_5_over_2p = SIMINT_MUL(const_5, one_over_2p);
372                     const SIMINT_DBLTYPE vrr_const_6_over_2p = SIMINT_MUL(const_6, one_over_2p);
373                     const SIMINT_DBLTYPE vrr_const_7_over_2p = SIMINT_MUL(const_7, one_over_2p);
374                     const SIMINT_DBLTYPE vrr_const_8_over_2p = SIMINT_MUL(const_8, one_over_2p);
375                     const SIMINT_DBLTYPE vrr_const_9_over_2p = SIMINT_MUL(const_9, one_over_2p);
376                     const SIMINT_DBLTYPE vrr_const_10_over_2p = SIMINT_MUL(const_10, one_over_2p);
377                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
378                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
379                     const SIMINT_DBLTYPE vrr_const_3_over_2q = SIMINT_MUL(const_3, one_over_2q);
380                     const SIMINT_DBLTYPE vrr_const_4_over_2q = SIMINT_MUL(const_4, one_over_2q);
381                     const SIMINT_DBLTYPE vrr_const_5_over_2q = SIMINT_MUL(const_5, one_over_2q);
382                     const SIMINT_DBLTYPE vrr_const_6_over_2q = SIMINT_MUL(const_6, one_over_2q);
383                     const SIMINT_DBLTYPE vrr_const_7_over_2q = SIMINT_MUL(const_7, one_over_2q);
384                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
385                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
386                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
387                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
388                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
389                     const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
390                     const SIMINT_DBLTYPE vrr_const_7_over_2pq = SIMINT_MUL(const_7, one_over_2pq);
391                     const SIMINT_DBLTYPE vrr_const_8_over_2pq = SIMINT_MUL(const_8, one_over_2pq);
392                     const SIMINT_DBLTYPE vrr_const_9_over_2pq = SIMINT_MUL(const_9, one_over_2pq);
393                     const SIMINT_DBLTYPE vrr_const_10_over_2pq = SIMINT_MUL(const_10, one_over_2pq);
394                     const SIMINT_DBLTYPE vrr_const_11_over_2pq = SIMINT_MUL(const_11, one_over_2pq);
395 
396 
397 
398                     // Forming PRIM_INT__p_s_s_s[19 * 3];
399                     for(n = 0; n < 19; ++n)  // loop over orders of auxiliary function
400                     {
401 
402                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
403                         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]);
404 
405                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
406                         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]);
407 
408                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
409                         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]);
410 
411                     }
412 
413 
414 
415                     // Forming PRIM_INT__d_s_s_s[18 * 6];
416                     for(n = 0; n < 18; ++n)  // loop over orders of auxiliary function
417                     {
418 
419                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
420                         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]);
421                         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]);
422 
423                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
424                         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]);
425 
426                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
427                         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]);
428 
429                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
430                         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]);
431                         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]);
432 
433                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
434                         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]);
435 
436                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
437                         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]);
438                         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]);
439 
440                     }
441 
442 
443 
444                     // Forming PRIM_INT__f_s_s_s[17 * 10];
445                     for(n = 0; n < 17; ++n)  // loop over orders of auxiliary function
446                     {
447 
448                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
449                         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]);
450                         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]);
451 
452                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
453                         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]);
454 
455                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
456                         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]);
457 
458                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
459                         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]);
460 
461                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
462                         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]);
463 
464                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
465                         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]);
466 
467                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
468                         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]);
469                         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]);
470 
471                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
472                         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]);
473 
474                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
475                         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]);
476 
477                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
478                         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]);
479                         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]);
480 
481                     }
482 
483 
484                     VRR_I_g_s_s_s(
485                             PRIM_INT__g_s_s_s,
486                             PRIM_INT__f_s_s_s,
487                             PRIM_INT__d_s_s_s,
488                             P_PA,
489                             a_over_p,
490                             aop_PQ,
491                             one_over_2p,
492                             16);
493 
494 
495                     VRR_I_h_s_s_s(
496                             PRIM_INT__h_s_s_s,
497                             PRIM_INT__g_s_s_s,
498                             PRIM_INT__f_s_s_s,
499                             P_PA,
500                             a_over_p,
501                             aop_PQ,
502                             one_over_2p,
503                             15);
504 
505 
506                     ostei_general_vrr1_I(6, 14,
507                             one_over_2p, a_over_p, aop_PQ, P_PA,
508                             PRIM_INT__h_s_s_s, PRIM_INT__g_s_s_s, PRIM_INT__i_s_s_s);
509 
510 
511                     ostei_general_vrr_K(6, 0, 1, 0, 8,
512                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
513                             PRIM_INT__i_s_s_s, NULL, NULL, PRIM_INT__h_s_s_s, NULL, PRIM_INT__i_s_p_s);
514 
515 
516                     ostei_general_vrr_K(5, 0, 1, 0, 8,
517                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
518                             PRIM_INT__h_s_s_s, NULL, NULL, PRIM_INT__g_s_s_s, NULL, PRIM_INT__h_s_p_s);
519 
520 
521                     ostei_general_vrr_K(6, 0, 2, 0, 7,
522                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
523                             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);
524 
525 
526                     VRR_K_g_s_p_s(
527                             PRIM_INT__g_s_p_s,
528                             PRIM_INT__g_s_s_s,
529                             PRIM_INT__f_s_s_s,
530                             Q_PA,
531                             aoq_PQ,
532                             one_over_2pq,
533                             8);
534 
535 
536                     ostei_general_vrr_K(5, 0, 2, 0, 7,
537                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
538                             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);
539 
540 
541                     ostei_general_vrr_K(6, 0, 3, 0, 6,
542                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
543                             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);
544 
545 
546                     VRR_K_f_s_p_s(
547                             PRIM_INT__f_s_p_s,
548                             PRIM_INT__f_s_s_s,
549                             PRIM_INT__d_s_s_s,
550                             Q_PA,
551                             aoq_PQ,
552                             one_over_2pq,
553                             8);
554 
555 
556                     ostei_general_vrr_K(4, 0, 2, 0, 7,
557                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
558                             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);
559 
560 
561                     ostei_general_vrr_K(5, 0, 3, 0, 6,
562                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
563                             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);
564 
565 
566                     ostei_general_vrr_K(6, 0, 4, 0, 5,
567                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
568                             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);
569 
570 
571 
572                     // Forming PRIM_INT__d_s_p_s[8 * 18];
573                     for(n = 0; n < 8; ++n)  // loop over orders of auxiliary function
574                     {
575 
576                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
577                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__d_s_p_s[n * 18 + 0]);
578                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_p_s[n * 18 + 0]);
579 
580                         PRIM_INT__d_s_p_s[n * 18 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
581                         PRIM_INT__d_s_p_s[n * 18 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__d_s_p_s[n * 18 + 1]);
582 
583                         PRIM_INT__d_s_p_s[n * 18 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
584                         PRIM_INT__d_s_p_s[n * 18 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__d_s_p_s[n * 18 + 2]);
585 
586                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 1]);
587                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__d_s_p_s[n * 18 + 3]);
588                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_p_s[n * 18 + 3]);
589 
590                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 1]);
591                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__d_s_p_s[n * 18 + 4]);
592                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_p_s[n * 18 + 4]);
593 
594                         PRIM_INT__d_s_p_s[n * 18 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
595                         PRIM_INT__d_s_p_s[n * 18 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__d_s_p_s[n * 18 + 5]);
596 
597                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 2]);
598                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 2], PRIM_INT__d_s_p_s[n * 18 + 6]);
599                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_p_s[n * 18 + 6]);
600 
601                         PRIM_INT__d_s_p_s[n * 18 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 2]);
602                         PRIM_INT__d_s_p_s[n * 18 + 7] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 2], PRIM_INT__d_s_p_s[n * 18 + 7]);
603 
604                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 2]);
605                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 2], PRIM_INT__d_s_p_s[n * 18 + 8]);
606                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_p_s[n * 18 + 8]);
607 
608                         PRIM_INT__d_s_p_s[n * 18 + 9] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
609                         PRIM_INT__d_s_p_s[n * 18 + 9] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__d_s_p_s[n * 18 + 9]);
610 
611                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
612                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__d_s_p_s[n * 18 + 10]);
613                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_p_s[n * 18 + 10]);
614 
615                         PRIM_INT__d_s_p_s[n * 18 + 11] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
616                         PRIM_INT__d_s_p_s[n * 18 + 11] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__d_s_p_s[n * 18 + 11]);
617 
618                         PRIM_INT__d_s_p_s[n * 18 + 12] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 4]);
619                         PRIM_INT__d_s_p_s[n * 18 + 12] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 4], PRIM_INT__d_s_p_s[n * 18 + 12]);
620 
621                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 4]);
622                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 4], PRIM_INT__d_s_p_s[n * 18 + 13]);
623                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_p_s[n * 18 + 13]);
624 
625                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 4]);
626                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 4], PRIM_INT__d_s_p_s[n * 18 + 14]);
627                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_p_s[n * 18 + 14]);
628 
629                         PRIM_INT__d_s_p_s[n * 18 + 15] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
630                         PRIM_INT__d_s_p_s[n * 18 + 15] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__d_s_p_s[n * 18 + 15]);
631 
632                         PRIM_INT__d_s_p_s[n * 18 + 16] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
633                         PRIM_INT__d_s_p_s[n * 18 + 16] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__d_s_p_s[n * 18 + 16]);
634 
635                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
636                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__d_s_p_s[n * 18 + 17]);
637                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_p_s[n * 18 + 17]);
638 
639                     }
640 
641 
642                     VRR_K_f_s_d_s(
643                             PRIM_INT__f_s_d_s,
644                             PRIM_INT__f_s_p_s,
645                             PRIM_INT__f_s_s_s,
646                             PRIM_INT__d_s_p_s,
647                             Q_PA,
648                             a_over_q,
649                             aoq_PQ,
650                             one_over_2pq,
651                             one_over_2q,
652                             7);
653 
654 
655                     ostei_general_vrr_K(4, 0, 3, 0, 6,
656                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
657                             PRIM_INT__g_s_d_s, PRIM_INT__g_s_p_s, NULL, PRIM_INT__f_s_d_s, NULL, PRIM_INT__g_s_f_s);
658 
659 
660                     ostei_general_vrr_K(5, 0, 4, 0, 5,
661                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
662                             PRIM_INT__h_s_f_s, PRIM_INT__h_s_d_s, NULL, PRIM_INT__g_s_f_s, NULL, PRIM_INT__h_s_g_s);
663 
664 
665                     ostei_general_vrr_K(6, 0, 5, 0, 4,
666                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
667                             PRIM_INT__i_s_g_s, PRIM_INT__i_s_f_s, NULL, PRIM_INT__h_s_g_s, NULL, PRIM_INT__i_s_h_s);
668 
669 
670 
671                     // Forming PRIM_INT__p_s_p_s[8 * 9];
672                     for(n = 0; n < 8; ++n)  // loop over orders of auxiliary function
673                     {
674 
675                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
676                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 0]);
677                         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]);
678 
679                         PRIM_INT__p_s_p_s[n * 9 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
680                         PRIM_INT__p_s_p_s[n * 9 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 1]);
681 
682                         PRIM_INT__p_s_p_s[n * 9 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
683                         PRIM_INT__p_s_p_s[n * 9 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 2]);
684 
685                         PRIM_INT__p_s_p_s[n * 9 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 1]);
686                         PRIM_INT__p_s_p_s[n * 9 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 3]);
687 
688                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
689                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 4]);
690                         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]);
691 
692                         PRIM_INT__p_s_p_s[n * 9 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
693                         PRIM_INT__p_s_p_s[n * 9 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 5]);
694 
695                         PRIM_INT__p_s_p_s[n * 9 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 2]);
696                         PRIM_INT__p_s_p_s[n * 9 + 6] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 6]);
697 
698                         PRIM_INT__p_s_p_s[n * 9 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 2]);
699                         PRIM_INT__p_s_p_s[n * 9 + 7] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 7]);
700 
701                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
702                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 8]);
703                         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]);
704 
705                     }
706 
707 
708                     VRR_K_d_s_d_s(
709                             PRIM_INT__d_s_d_s,
710                             PRIM_INT__d_s_p_s,
711                             PRIM_INT__d_s_s_s,
712                             PRIM_INT__p_s_p_s,
713                             Q_PA,
714                             a_over_q,
715                             aoq_PQ,
716                             one_over_2pq,
717                             one_over_2q,
718                             7);
719 
720 
721                     ostei_general_vrr_K(3, 0, 3, 0, 6,
722                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
723                             PRIM_INT__f_s_d_s, PRIM_INT__f_s_p_s, NULL, PRIM_INT__d_s_d_s, NULL, PRIM_INT__f_s_f_s);
724 
725 
726                     ostei_general_vrr_K(4, 0, 4, 0, 5,
727                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
728                             PRIM_INT__g_s_f_s, PRIM_INT__g_s_d_s, NULL, PRIM_INT__f_s_f_s, NULL, PRIM_INT__g_s_g_s);
729 
730 
731                     ostei_general_vrr_K(5, 0, 5, 0, 4,
732                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
733                             PRIM_INT__h_s_g_s, PRIM_INT__h_s_f_s, NULL, PRIM_INT__g_s_g_s, NULL, PRIM_INT__h_s_h_s);
734 
735 
736                     ostei_general_vrr_K(6, 0, 6, 0, 3,
737                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
738                             PRIM_INT__i_s_h_s, PRIM_INT__i_s_g_s, NULL, PRIM_INT__h_s_h_s, NULL, PRIM_INT__i_s_i_s);
739 
740 
741 
742                     // Forming PRIM_INT__s_s_p_s[8 * 3];
743                     for(n = 0; n < 8; ++n)  // loop over orders of auxiliary function
744                     {
745 
746                         PRIM_INT__s_s_p_s[n * 3 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
747                         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]);
748 
749                         PRIM_INT__s_s_p_s[n * 3 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
750                         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]);
751 
752                         PRIM_INT__s_s_p_s[n * 3 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
753                         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]);
754 
755                     }
756 
757 
758 
759                     // Forming PRIM_INT__p_s_d_s[7 * 18];
760                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
761                     {
762 
763                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_p_s[n * 9 + 0]);
764                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_p_s[(n+1) * 9 + 0], PRIM_INT__p_s_d_s[n * 18 + 0]);
765                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]), PRIM_INT__p_s_d_s[n * 18 + 0]);
766                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_d_s[n * 18 + 0]);
767 
768                         PRIM_INT__p_s_d_s[n * 18 + 3] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_p_s[n * 9 + 1]);
769                         PRIM_INT__p_s_d_s[n * 18 + 3] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_p_s[(n+1) * 9 + 1], PRIM_INT__p_s_d_s[n * 18 + 3]);
770                         PRIM_INT__p_s_d_s[n * 18 + 3] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]), PRIM_INT__p_s_d_s[n * 18 + 3]);
771 
772                         PRIM_INT__p_s_d_s[n * 18 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_p_s[n * 9 + 2]);
773                         PRIM_INT__p_s_d_s[n * 18 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_p_s[(n+1) * 9 + 2], PRIM_INT__p_s_d_s[n * 18 + 5]);
774                         PRIM_INT__p_s_d_s[n * 18 + 5] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]), PRIM_INT__p_s_d_s[n * 18 + 5]);
775 
776                         PRIM_INT__p_s_d_s[n * 18 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_p_s[n * 9 + 3]);
777                         PRIM_INT__p_s_d_s[n * 18 + 6] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_p_s[(n+1) * 9 + 3], PRIM_INT__p_s_d_s[n * 18 + 6]);
778                         PRIM_INT__p_s_d_s[n * 18 + 6] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_s[n * 3 + 1]), PRIM_INT__p_s_d_s[n * 18 + 6]);
779 
780                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_p_s[n * 9 + 4]);
781                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_p_s[(n+1) * 9 + 4], PRIM_INT__p_s_d_s[n * 18 + 9]);
782                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_s[n * 3 + 1]), PRIM_INT__p_s_d_s[n * 18 + 9]);
783                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_d_s[n * 18 + 9]);
784 
785                         PRIM_INT__p_s_d_s[n * 18 + 11] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_p_s[n * 9 + 5]);
786                         PRIM_INT__p_s_d_s[n * 18 + 11] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_p_s[(n+1) * 9 + 5], PRIM_INT__p_s_d_s[n * 18 + 11]);
787                         PRIM_INT__p_s_d_s[n * 18 + 11] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_s[n * 3 + 1]), PRIM_INT__p_s_d_s[n * 18 + 11]);
788 
789                         PRIM_INT__p_s_d_s[n * 18 + 12] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_p_s[n * 9 + 6]);
790                         PRIM_INT__p_s_d_s[n * 18 + 12] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_p_s[(n+1) * 9 + 6], PRIM_INT__p_s_d_s[n * 18 + 12]);
791                         PRIM_INT__p_s_d_s[n * 18 + 12] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_s[n * 3 + 2]), PRIM_INT__p_s_d_s[n * 18 + 12]);
792 
793                         PRIM_INT__p_s_d_s[n * 18 + 15] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_p_s[n * 9 + 7]);
794                         PRIM_INT__p_s_d_s[n * 18 + 15] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_p_s[(n+1) * 9 + 7], PRIM_INT__p_s_d_s[n * 18 + 15]);
795                         PRIM_INT__p_s_d_s[n * 18 + 15] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_s[n * 3 + 2]), PRIM_INT__p_s_d_s[n * 18 + 15]);
796 
797                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_p_s[n * 9 + 8]);
798                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_p_s[(n+1) * 9 + 8], PRIM_INT__p_s_d_s[n * 18 + 17]);
799                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_s[n * 3 + 2]), PRIM_INT__p_s_d_s[n * 18 + 17]);
800                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_d_s[n * 18 + 17]);
801 
802                     }
803 
804 
805                     VRR_K_d_s_f_s(
806                             PRIM_INT__d_s_f_s,
807                             PRIM_INT__d_s_d_s,
808                             PRIM_INT__d_s_p_s,
809                             PRIM_INT__p_s_d_s,
810                             Q_PA,
811                             a_over_q,
812                             aoq_PQ,
813                             one_over_2pq,
814                             one_over_2q,
815                             6);
816 
817 
818                     ostei_general_vrr_K(3, 0, 4, 0, 5,
819                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
820                             PRIM_INT__f_s_f_s, PRIM_INT__f_s_d_s, NULL, PRIM_INT__d_s_f_s, NULL, PRIM_INT__f_s_g_s);
821 
822 
823                     ostei_general_vrr_K(4, 0, 5, 0, 4,
824                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
825                             PRIM_INT__g_s_g_s, PRIM_INT__g_s_f_s, NULL, PRIM_INT__f_s_g_s, NULL, PRIM_INT__g_s_h_s);
826 
827 
828                     ostei_general_vrr_K(5, 0, 6, 0, 3,
829                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
830                             PRIM_INT__h_s_h_s, PRIM_INT__h_s_g_s, NULL, PRIM_INT__g_s_h_s, NULL, PRIM_INT__h_s_i_s);
831 
832 
833                     ostei_general_vrr_K(6, 0, 7, 0, 2,
834                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
835                             PRIM_INT__i_s_i_s, PRIM_INT__i_s_h_s, NULL, PRIM_INT__h_s_i_s, NULL, PRIM_INT__i_s_k_s);
836 
837 
838 
839                     // Forming PRIM_INT__s_s_d_s[7 * 6];
840                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
841                     {
842 
843                         PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_p_s[n * 3 + 0]);
844                         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]);
845                         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]);
846 
847                         PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_p_s[n * 3 + 1]);
848                         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]);
849                         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]);
850 
851                         PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 2]);
852                         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]);
853                         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]);
854 
855                     }
856 
857 
858                     VRR_K_p_s_f_s(
859                             PRIM_INT__p_s_f_s,
860                             PRIM_INT__p_s_d_s,
861                             PRIM_INT__p_s_p_s,
862                             PRIM_INT__s_s_d_s,
863                             Q_PA,
864                             a_over_q,
865                             aoq_PQ,
866                             one_over_2pq,
867                             one_over_2q,
868                             6);
869 
870 
871                     ostei_general_vrr_K(2, 0, 4, 0, 5,
872                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
873                             PRIM_INT__d_s_f_s, PRIM_INT__d_s_d_s, NULL, PRIM_INT__p_s_f_s, NULL, PRIM_INT__d_s_g_s);
874 
875 
876                     ostei_general_vrr_K(3, 0, 5, 0, 4,
877                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
878                             PRIM_INT__f_s_g_s, PRIM_INT__f_s_f_s, NULL, PRIM_INT__d_s_g_s, NULL, PRIM_INT__f_s_h_s);
879 
880 
881                     ostei_general_vrr_K(4, 0, 6, 0, 3,
882                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
883                             PRIM_INT__g_s_h_s, PRIM_INT__g_s_g_s, NULL, PRIM_INT__f_s_h_s, NULL, PRIM_INT__g_s_i_s);
884 
885 
886                     ostei_general_vrr_K(5, 0, 7, 0, 2,
887                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
888                             PRIM_INT__h_s_i_s, PRIM_INT__h_s_h_s, NULL, PRIM_INT__g_s_i_s, NULL, PRIM_INT__h_s_k_s);
889 
890 
891                     ostei_general_vrr_K(6, 0, 8, 0, 1,
892                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
893                             PRIM_INT__i_s_k_s, PRIM_INT__i_s_i_s, NULL, PRIM_INT__h_s_k_s, NULL, PRIM_INT__i_s_l_s);
894 
895 
896                     ostei_general_vrr1_I(7, 13,
897                             one_over_2p, a_over_p, aop_PQ, P_PA,
898                             PRIM_INT__i_s_s_s, PRIM_INT__h_s_s_s, PRIM_INT__k_s_s_s);
899 
900 
901                     ostei_general_vrr_K(7, 0, 1, 0, 8,
902                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
903                             PRIM_INT__k_s_s_s, NULL, NULL, PRIM_INT__i_s_s_s, NULL, PRIM_INT__k_s_p_s);
904 
905 
906                     ostei_general_vrr_K(7, 0, 2, 0, 7,
907                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
908                             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);
909 
910 
911                     ostei_general_vrr_K(7, 0, 3, 0, 6,
912                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
913                             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);
914 
915 
916                     ostei_general_vrr_K(7, 0, 4, 0, 5,
917                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
918                             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);
919 
920 
921                     ostei_general_vrr_K(7, 0, 5, 0, 4,
922                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
923                             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);
924 
925 
926                     ostei_general_vrr_K(7, 0, 6, 0, 3,
927                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
928                             PRIM_INT__k_s_h_s, PRIM_INT__k_s_g_s, NULL, PRIM_INT__i_s_h_s, NULL, PRIM_INT__k_s_i_s);
929 
930 
931                     ostei_general_vrr_K(7, 0, 7, 0, 2,
932                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
933                             PRIM_INT__k_s_i_s, PRIM_INT__k_s_h_s, NULL, PRIM_INT__i_s_i_s, NULL, PRIM_INT__k_s_k_s);
934 
935 
936                     ostei_general_vrr_K(7, 0, 8, 0, 1,
937                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
938                             PRIM_INT__k_s_k_s, PRIM_INT__k_s_i_s, NULL, PRIM_INT__i_s_k_s, NULL, PRIM_INT__k_s_l_s);
939 
940 
941                     ostei_general_vrr1_I(8, 12,
942                             one_over_2p, a_over_p, aop_PQ, P_PA,
943                             PRIM_INT__k_s_s_s, PRIM_INT__i_s_s_s, PRIM_INT__l_s_s_s);
944 
945 
946                     ostei_general_vrr_K(8, 0, 1, 0, 8,
947                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
948                             PRIM_INT__l_s_s_s, NULL, NULL, PRIM_INT__k_s_s_s, NULL, PRIM_INT__l_s_p_s);
949 
950 
951                     ostei_general_vrr_K(8, 0, 2, 0, 7,
952                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
953                             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);
954 
955 
956                     ostei_general_vrr_K(8, 0, 3, 0, 6,
957                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
958                             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);
959 
960 
961                     ostei_general_vrr_K(8, 0, 4, 0, 5,
962                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
963                             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);
964 
965 
966                     ostei_general_vrr_K(8, 0, 5, 0, 4,
967                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
968                             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);
969 
970 
971                     ostei_general_vrr_K(8, 0, 6, 0, 3,
972                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
973                             PRIM_INT__l_s_h_s, PRIM_INT__l_s_g_s, NULL, PRIM_INT__k_s_h_s, NULL, PRIM_INT__l_s_i_s);
974 
975 
976                     ostei_general_vrr_K(8, 0, 7, 0, 2,
977                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
978                             PRIM_INT__l_s_i_s, PRIM_INT__l_s_h_s, NULL, PRIM_INT__k_s_i_s, NULL, PRIM_INT__l_s_k_s);
979 
980 
981                     ostei_general_vrr_K(8, 0, 8, 0, 1,
982                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
983                             PRIM_INT__l_s_k_s, PRIM_INT__l_s_i_s, NULL, PRIM_INT__k_s_k_s, NULL, PRIM_INT__l_s_l_s);
984 
985 
986                     ostei_general_vrr1_I(9, 11,
987                             one_over_2p, a_over_p, aop_PQ, P_PA,
988                             PRIM_INT__l_s_s_s, PRIM_INT__k_s_s_s, PRIM_INT__m_s_s_s);
989 
990 
991                     ostei_general_vrr_K(9, 0, 1, 0, 8,
992                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
993                             PRIM_INT__m_s_s_s, NULL, NULL, PRIM_INT__l_s_s_s, NULL, PRIM_INT__m_s_p_s);
994 
995 
996                     ostei_general_vrr_K(9, 0, 2, 0, 7,
997                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
998                             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);
999 
1000 
1001                     ostei_general_vrr_K(9, 0, 3, 0, 6,
1002                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1003                             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);
1004 
1005 
1006                     ostei_general_vrr_K(9, 0, 4, 0, 5,
1007                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1008                             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);
1009 
1010 
1011                     ostei_general_vrr_K(9, 0, 5, 0, 4,
1012                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1013                             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);
1014 
1015 
1016                     ostei_general_vrr_K(9, 0, 6, 0, 3,
1017                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1018                             PRIM_INT__m_s_h_s, PRIM_INT__m_s_g_s, NULL, PRIM_INT__l_s_h_s, NULL, PRIM_INT__m_s_i_s);
1019 
1020 
1021                     ostei_general_vrr_K(9, 0, 7, 0, 2,
1022                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1023                             PRIM_INT__m_s_i_s, PRIM_INT__m_s_h_s, NULL, PRIM_INT__l_s_i_s, NULL, PRIM_INT__m_s_k_s);
1024 
1025 
1026                     ostei_general_vrr_K(9, 0, 8, 0, 1,
1027                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1028                             PRIM_INT__m_s_k_s, PRIM_INT__m_s_i_s, NULL, PRIM_INT__l_s_k_s, NULL, PRIM_INT__m_s_l_s);
1029 
1030 
1031                     ostei_general_vrr1_I(10, 10,
1032                             one_over_2p, a_over_p, aop_PQ, P_PA,
1033                             PRIM_INT__m_s_s_s, PRIM_INT__l_s_s_s, PRIM_INT__n_s_s_s);
1034 
1035 
1036                     ostei_general_vrr_K(10, 0, 1, 0, 8,
1037                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1038                             PRIM_INT__n_s_s_s, NULL, NULL, PRIM_INT__m_s_s_s, NULL, PRIM_INT__n_s_p_s);
1039 
1040 
1041                     ostei_general_vrr_K(10, 0, 2, 0, 7,
1042                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1043                             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);
1044 
1045 
1046                     ostei_general_vrr_K(10, 0, 3, 0, 6,
1047                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1048                             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);
1049 
1050 
1051                     ostei_general_vrr_K(10, 0, 4, 0, 5,
1052                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1053                             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);
1054 
1055 
1056                     ostei_general_vrr_K(10, 0, 5, 0, 4,
1057                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1058                             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);
1059 
1060 
1061                     ostei_general_vrr_K(10, 0, 6, 0, 3,
1062                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1063                             PRIM_INT__n_s_h_s, PRIM_INT__n_s_g_s, NULL, PRIM_INT__m_s_h_s, NULL, PRIM_INT__n_s_i_s);
1064 
1065 
1066                     ostei_general_vrr_K(10, 0, 7, 0, 2,
1067                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1068                             PRIM_INT__n_s_i_s, PRIM_INT__n_s_h_s, NULL, PRIM_INT__m_s_i_s, NULL, PRIM_INT__n_s_k_s);
1069 
1070 
1071                     ostei_general_vrr_K(10, 0, 8, 0, 1,
1072                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1073                             PRIM_INT__n_s_k_s, PRIM_INT__n_s_i_s, NULL, PRIM_INT__m_s_k_s, NULL, PRIM_INT__n_s_l_s);
1074 
1075 
1076                     ostei_general_vrr1_I(11, 9,
1077                             one_over_2p, a_over_p, aop_PQ, P_PA,
1078                             PRIM_INT__n_s_s_s, PRIM_INT__m_s_s_s, PRIM_INT__o_s_s_s);
1079 
1080 
1081                     ostei_general_vrr_K(11, 0, 1, 0, 8,
1082                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1083                             PRIM_INT__o_s_s_s, NULL, NULL, PRIM_INT__n_s_s_s, NULL, PRIM_INT__o_s_p_s);
1084 
1085 
1086                     ostei_general_vrr_K(11, 0, 2, 0, 7,
1087                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1088                             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);
1089 
1090 
1091                     ostei_general_vrr_K(11, 0, 3, 0, 6,
1092                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1093                             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);
1094 
1095 
1096                     ostei_general_vrr_K(11, 0, 4, 0, 5,
1097                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1098                             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);
1099 
1100 
1101                     ostei_general_vrr_K(11, 0, 5, 0, 4,
1102                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1103                             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);
1104 
1105 
1106                     ostei_general_vrr_K(11, 0, 6, 0, 3,
1107                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1108                             PRIM_INT__o_s_h_s, PRIM_INT__o_s_g_s, NULL, PRIM_INT__n_s_h_s, NULL, PRIM_INT__o_s_i_s);
1109 
1110 
1111                     ostei_general_vrr_K(11, 0, 7, 0, 2,
1112                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1113                             PRIM_INT__o_s_i_s, PRIM_INT__o_s_h_s, NULL, PRIM_INT__n_s_i_s, NULL, PRIM_INT__o_s_k_s);
1114 
1115 
1116                     ostei_general_vrr_K(11, 0, 8, 0, 1,
1117                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1118                             PRIM_INT__o_s_k_s, PRIM_INT__o_s_i_s, NULL, PRIM_INT__n_s_k_s, NULL, PRIM_INT__o_s_l_s);
1119 
1120 
1121 
1122 
1123                     ////////////////////////////////////
1124                     // Accumulate contracted integrals
1125                     ////////////////////////////////////
1126                     if(lastoffset == 0)
1127                     {
1128                         contract_all(1008, PRIM_INT__i_s_k_s, PRIM_PTR_INT__i_s_k_s);
1129                         contract_all(1260, PRIM_INT__i_s_l_s, PRIM_PTR_INT__i_s_l_s);
1130                         contract_all(1296, PRIM_INT__k_s_k_s, PRIM_PTR_INT__k_s_k_s);
1131                         contract_all(1620, PRIM_INT__k_s_l_s, PRIM_PTR_INT__k_s_l_s);
1132                         contract_all(1620, PRIM_INT__l_s_k_s, PRIM_PTR_INT__l_s_k_s);
1133                         contract_all(2025, PRIM_INT__l_s_l_s, PRIM_PTR_INT__l_s_l_s);
1134                         contract_all(1980, PRIM_INT__m_s_k_s, PRIM_PTR_INT__m_s_k_s);
1135                         contract_all(2475, PRIM_INT__m_s_l_s, PRIM_PTR_INT__m_s_l_s);
1136                         contract_all(2376, PRIM_INT__n_s_k_s, PRIM_PTR_INT__n_s_k_s);
1137                         contract_all(2970, PRIM_INT__n_s_l_s, PRIM_PTR_INT__n_s_l_s);
1138                         contract_all(2808, PRIM_INT__o_s_k_s, PRIM_PTR_INT__o_s_k_s);
1139                         contract_all(3510, PRIM_INT__o_s_l_s, PRIM_PTR_INT__o_s_l_s);
1140                     }
1141                     else
1142                     {
1143                         contract(1008, shelloffsets, PRIM_INT__i_s_k_s, PRIM_PTR_INT__i_s_k_s);
1144                         contract(1260, shelloffsets, PRIM_INT__i_s_l_s, PRIM_PTR_INT__i_s_l_s);
1145                         contract(1296, shelloffsets, PRIM_INT__k_s_k_s, PRIM_PTR_INT__k_s_k_s);
1146                         contract(1620, shelloffsets, PRIM_INT__k_s_l_s, PRIM_PTR_INT__k_s_l_s);
1147                         contract(1620, shelloffsets, PRIM_INT__l_s_k_s, PRIM_PTR_INT__l_s_k_s);
1148                         contract(2025, shelloffsets, PRIM_INT__l_s_l_s, PRIM_PTR_INT__l_s_l_s);
1149                         contract(1980, shelloffsets, PRIM_INT__m_s_k_s, PRIM_PTR_INT__m_s_k_s);
1150                         contract(2475, shelloffsets, PRIM_INT__m_s_l_s, PRIM_PTR_INT__m_s_l_s);
1151                         contract(2376, shelloffsets, PRIM_INT__n_s_k_s, PRIM_PTR_INT__n_s_k_s);
1152                         contract(2970, shelloffsets, PRIM_INT__n_s_l_s, PRIM_PTR_INT__n_s_l_s);
1153                         contract(2808, shelloffsets, PRIM_INT__o_s_k_s, PRIM_PTR_INT__o_s_k_s);
1154                         contract(3510, shelloffsets, PRIM_INT__o_s_l_s, PRIM_PTR_INT__o_s_l_s);
1155                         PRIM_PTR_INT__i_s_k_s += lastoffset*1008;
1156                         PRIM_PTR_INT__i_s_l_s += lastoffset*1260;
1157                         PRIM_PTR_INT__k_s_k_s += lastoffset*1296;
1158                         PRIM_PTR_INT__k_s_l_s += lastoffset*1620;
1159                         PRIM_PTR_INT__l_s_k_s += lastoffset*1620;
1160                         PRIM_PTR_INT__l_s_l_s += lastoffset*2025;
1161                         PRIM_PTR_INT__m_s_k_s += lastoffset*1980;
1162                         PRIM_PTR_INT__m_s_l_s += lastoffset*2475;
1163                         PRIM_PTR_INT__n_s_k_s += lastoffset*2376;
1164                         PRIM_PTR_INT__n_s_l_s += lastoffset*2970;
1165                         PRIM_PTR_INT__o_s_k_s += lastoffset*2808;
1166                         PRIM_PTR_INT__o_s_l_s += lastoffset*3510;
1167                     }
1168 
1169                 }  // close loop over j
1170             }  // close loop over i
1171 
1172             //Advance to the next batch
1173             jstart = SIMINT_SIMD_ROUND(jend);
1174 
1175             //////////////////////////////////////////////
1176             // Contracted integrals: Horizontal recurrance
1177             //////////////////////////////////////////////
1178 
1179 
1180             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
1181 
1182 
1183             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
1184             {
1185                 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
1186 
1187                 // set up HRR pointers
1188                 double const * restrict HRR_INT__i_s_k_s = INT__i_s_k_s + abcd * 1008;
1189                 double const * restrict HRR_INT__i_s_l_s = INT__i_s_l_s + abcd * 1260;
1190                 double const * restrict HRR_INT__k_s_k_s = INT__k_s_k_s + abcd * 1296;
1191                 double const * restrict HRR_INT__k_s_l_s = INT__k_s_l_s + abcd * 1620;
1192                 double const * restrict HRR_INT__l_s_k_s = INT__l_s_k_s + abcd * 1620;
1193                 double const * restrict HRR_INT__l_s_l_s = INT__l_s_l_s + abcd * 2025;
1194                 double const * restrict HRR_INT__m_s_k_s = INT__m_s_k_s + abcd * 1980;
1195                 double const * restrict HRR_INT__m_s_l_s = INT__m_s_l_s + abcd * 2475;
1196                 double const * restrict HRR_INT__n_s_k_s = INT__n_s_k_s + abcd * 2376;
1197                 double const * restrict HRR_INT__n_s_l_s = INT__n_s_l_s + abcd * 2970;
1198                 double const * restrict HRR_INT__o_s_k_s = INT__o_s_k_s + abcd * 2808;
1199                 double const * restrict HRR_INT__o_s_l_s = INT__o_s_l_s + abcd * 3510;
1200                 double * restrict HRR_INT__i_h_k_p = INT__i_h_k_p + real_abcd * 63504;
1201 
1202                 // form INT__i_p_k_s
1203                 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);
1204 
1205                 // form INT__i_p_l_s
1206                 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);
1207 
1208                 // form INT__k_p_k_s
1209                 ostei_general_hrr_J(7, 1, 7, 0, hAB, HRR_INT__l_s_k_s, HRR_INT__k_s_k_s, HRR_INT__k_p_k_s);
1210 
1211                 // form INT__k_p_l_s
1212                 ostei_general_hrr_J(7, 1, 8, 0, hAB, HRR_INT__l_s_l_s, HRR_INT__k_s_l_s, HRR_INT__k_p_l_s);
1213 
1214                 // form INT__l_p_k_s
1215                 ostei_general_hrr_J(8, 1, 7, 0, hAB, HRR_INT__m_s_k_s, HRR_INT__l_s_k_s, HRR_INT__l_p_k_s);
1216 
1217                 // form INT__l_p_l_s
1218                 ostei_general_hrr_J(8, 1, 8, 0, hAB, HRR_INT__m_s_l_s, HRR_INT__l_s_l_s, HRR_INT__l_p_l_s);
1219 
1220                 // form INT__m_p_k_s
1221                 ostei_general_hrr_J(9, 1, 7, 0, hAB, HRR_INT__n_s_k_s, HRR_INT__m_s_k_s, HRR_INT__m_p_k_s);
1222 
1223                 // form INT__m_p_l_s
1224                 ostei_general_hrr_J(9, 1, 8, 0, hAB, HRR_INT__n_s_l_s, HRR_INT__m_s_l_s, HRR_INT__m_p_l_s);
1225 
1226                 // form INT__n_p_k_s
1227                 ostei_general_hrr_J(10, 1, 7, 0, hAB, HRR_INT__o_s_k_s, HRR_INT__n_s_k_s, HRR_INT__n_p_k_s);
1228 
1229                 // form INT__n_p_l_s
1230                 ostei_general_hrr_J(10, 1, 8, 0, hAB, HRR_INT__o_s_l_s, HRR_INT__n_s_l_s, HRR_INT__n_p_l_s);
1231 
1232                 // form INT__i_d_k_s
1233                 ostei_general_hrr_J(6, 2, 7, 0, hAB, HRR_INT__k_p_k_s, HRR_INT__i_p_k_s, HRR_INT__i_d_k_s);
1234 
1235                 // form INT__i_d_l_s
1236                 ostei_general_hrr_J(6, 2, 8, 0, hAB, HRR_INT__k_p_l_s, HRR_INT__i_p_l_s, HRR_INT__i_d_l_s);
1237 
1238                 // form INT__k_d_k_s
1239                 ostei_general_hrr_J(7, 2, 7, 0, hAB, HRR_INT__l_p_k_s, HRR_INT__k_p_k_s, HRR_INT__k_d_k_s);
1240 
1241                 // form INT__k_d_l_s
1242                 ostei_general_hrr_J(7, 2, 8, 0, hAB, HRR_INT__l_p_l_s, HRR_INT__k_p_l_s, HRR_INT__k_d_l_s);
1243 
1244                 // form INT__l_d_k_s
1245                 ostei_general_hrr_J(8, 2, 7, 0, hAB, HRR_INT__m_p_k_s, HRR_INT__l_p_k_s, HRR_INT__l_d_k_s);
1246 
1247                 // form INT__l_d_l_s
1248                 ostei_general_hrr_J(8, 2, 8, 0, hAB, HRR_INT__m_p_l_s, HRR_INT__l_p_l_s, HRR_INT__l_d_l_s);
1249 
1250                 // form INT__m_d_k_s
1251                 ostei_general_hrr_J(9, 2, 7, 0, hAB, HRR_INT__n_p_k_s, HRR_INT__m_p_k_s, HRR_INT__m_d_k_s);
1252 
1253                 // form INT__m_d_l_s
1254                 ostei_general_hrr_J(9, 2, 8, 0, hAB, HRR_INT__n_p_l_s, HRR_INT__m_p_l_s, HRR_INT__m_d_l_s);
1255 
1256                 // form INT__i_f_k_s
1257                 ostei_general_hrr_J(6, 3, 7, 0, hAB, HRR_INT__k_d_k_s, HRR_INT__i_d_k_s, HRR_INT__i_f_k_s);
1258 
1259                 // form INT__i_f_l_s
1260                 ostei_general_hrr_J(6, 3, 8, 0, hAB, HRR_INT__k_d_l_s, HRR_INT__i_d_l_s, HRR_INT__i_f_l_s);
1261 
1262                 // form INT__k_f_k_s
1263                 ostei_general_hrr_J(7, 3, 7, 0, hAB, HRR_INT__l_d_k_s, HRR_INT__k_d_k_s, HRR_INT__k_f_k_s);
1264 
1265                 // form INT__k_f_l_s
1266                 ostei_general_hrr_J(7, 3, 8, 0, hAB, HRR_INT__l_d_l_s, HRR_INT__k_d_l_s, HRR_INT__k_f_l_s);
1267 
1268                 // form INT__l_f_k_s
1269                 ostei_general_hrr_J(8, 3, 7, 0, hAB, HRR_INT__m_d_k_s, HRR_INT__l_d_k_s, HRR_INT__l_f_k_s);
1270 
1271                 // form INT__l_f_l_s
1272                 ostei_general_hrr_J(8, 3, 8, 0, hAB, HRR_INT__m_d_l_s, HRR_INT__l_d_l_s, HRR_INT__l_f_l_s);
1273 
1274                 // form INT__i_g_k_s
1275                 ostei_general_hrr_J(6, 4, 7, 0, hAB, HRR_INT__k_f_k_s, HRR_INT__i_f_k_s, HRR_INT__i_g_k_s);
1276 
1277                 // form INT__i_g_l_s
1278                 ostei_general_hrr_J(6, 4, 8, 0, hAB, HRR_INT__k_f_l_s, HRR_INT__i_f_l_s, HRR_INT__i_g_l_s);
1279 
1280                 // form INT__k_g_k_s
1281                 ostei_general_hrr_J(7, 4, 7, 0, hAB, HRR_INT__l_f_k_s, HRR_INT__k_f_k_s, HRR_INT__k_g_k_s);
1282 
1283                 // form INT__k_g_l_s
1284                 ostei_general_hrr_J(7, 4, 8, 0, hAB, HRR_INT__l_f_l_s, HRR_INT__k_f_l_s, HRR_INT__k_g_l_s);
1285 
1286                 // form INT__i_h_k_s
1287                 ostei_general_hrr_J(6, 5, 7, 0, hAB, HRR_INT__k_g_k_s, HRR_INT__i_g_k_s, HRR_INT__i_h_k_s);
1288 
1289                 // form INT__i_h_l_s
1290                 ostei_general_hrr_J(6, 5, 8, 0, hAB, HRR_INT__k_g_l_s, HRR_INT__i_g_l_s, HRR_INT__i_h_l_s);
1291 
1292                 // form INT__i_h_k_p
1293                 ostei_general_hrr_L(6, 5, 7, 1, hCD, HRR_INT__i_h_l_s, HRR_INT__i_h_k_s, HRR_INT__i_h_k_p);
1294 
1295 
1296             }  // close HRR loop
1297 
1298 
1299         }   // close loop cdbatch
1300 
1301         istart = iend;
1302     }  // close loop over ab
1303 
1304     return P.nshell12_clip * Q.nshell12_clip;
1305 }
1306 
ostei_h_i_k_p(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__h_i_k_p)1307 int ostei_h_i_k_p(struct simint_multi_shellpair const P,
1308                   struct simint_multi_shellpair const Q,
1309                   double screen_tol,
1310                   double * const restrict work,
1311                   double * const restrict INT__h_i_k_p)
1312 {
1313     double P_AB[3*P.nshell12];
1314     struct simint_multi_shellpair P_tmp = P;
1315     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
1316     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
1317     P_tmp.AB_x = P_AB;
1318     P_tmp.AB_y = P_AB + P.nshell12;
1319     P_tmp.AB_z = P_AB + 2*P.nshell12;
1320 
1321     for(int i = 0; i < P.nshell12; i++)
1322     {
1323         P_tmp.AB_x[i] = -P.AB_x[i];
1324         P_tmp.AB_y[i] = -P.AB_y[i];
1325         P_tmp.AB_z[i] = -P.AB_z[i];
1326     }
1327 
1328     int ret = ostei_i_h_k_p(P_tmp, Q, screen_tol, work, INT__h_i_k_p);
1329     double buffer[63504] SIMINT_ALIGN_ARRAY_DBL;
1330 
1331     for(int q = 0; q < ret; q++)
1332     {
1333         int idx = 0;
1334         for(int a = 0; a < 21; ++a)
1335         for(int b = 0; b < 28; ++b)
1336         for(int c = 0; c < 36; ++c)
1337         for(int d = 0; d < 3; ++d)
1338             buffer[idx++] = INT__h_i_k_p[q*63504+b*2268+a*108+c*3+d];
1339 
1340         memcpy(INT__h_i_k_p+q*63504, buffer, 63504*sizeof(double));
1341     }
1342 
1343     return ret;
1344 }
1345 
ostei_i_h_p_k(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__i_h_p_k)1346 int ostei_i_h_p_k(struct simint_multi_shellpair const P,
1347                   struct simint_multi_shellpair const Q,
1348                   double screen_tol,
1349                   double * const restrict work,
1350                   double * const restrict INT__i_h_p_k)
1351 {
1352     double Q_AB[3*Q.nshell12];
1353     struct simint_multi_shellpair Q_tmp = Q;
1354     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
1355     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
1356     Q_tmp.AB_x = Q_AB;
1357     Q_tmp.AB_y = Q_AB + Q.nshell12;
1358     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
1359 
1360     for(int i = 0; i < Q.nshell12; i++)
1361     {
1362         Q_tmp.AB_x[i] = -Q.AB_x[i];
1363         Q_tmp.AB_y[i] = -Q.AB_y[i];
1364         Q_tmp.AB_z[i] = -Q.AB_z[i];
1365     }
1366 
1367     int ret = ostei_i_h_k_p(P, Q_tmp, screen_tol, work, INT__i_h_p_k);
1368     double buffer[63504] SIMINT_ALIGN_ARRAY_DBL;
1369 
1370     for(int q = 0; q < ret; q++)
1371     {
1372         int idx = 0;
1373         for(int a = 0; a < 28; ++a)
1374         for(int b = 0; b < 21; ++b)
1375         for(int c = 0; c < 3; ++c)
1376         for(int d = 0; d < 36; ++d)
1377             buffer[idx++] = INT__i_h_p_k[q*63504+a*2268+b*108+d*3+c];
1378 
1379         memcpy(INT__i_h_p_k+q*63504, buffer, 63504*sizeof(double));
1380     }
1381 
1382     return ret;
1383 }
1384 
ostei_h_i_p_k(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__h_i_p_k)1385 int ostei_h_i_p_k(struct simint_multi_shellpair const P,
1386                   struct simint_multi_shellpair const Q,
1387                   double screen_tol,
1388                   double * const restrict work,
1389                   double * const restrict INT__h_i_p_k)
1390 {
1391     double P_AB[3*P.nshell12];
1392     struct simint_multi_shellpair P_tmp = P;
1393     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
1394     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
1395     P_tmp.AB_x = P_AB;
1396     P_tmp.AB_y = P_AB + P.nshell12;
1397     P_tmp.AB_z = P_AB + 2*P.nshell12;
1398 
1399     for(int i = 0; i < P.nshell12; i++)
1400     {
1401         P_tmp.AB_x[i] = -P.AB_x[i];
1402         P_tmp.AB_y[i] = -P.AB_y[i];
1403         P_tmp.AB_z[i] = -P.AB_z[i];
1404     }
1405 
1406     double Q_AB[3*Q.nshell12];
1407     struct simint_multi_shellpair Q_tmp = Q;
1408     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
1409     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
1410     Q_tmp.AB_x = Q_AB;
1411     Q_tmp.AB_y = Q_AB + Q.nshell12;
1412     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
1413 
1414     for(int i = 0; i < Q.nshell12; i++)
1415     {
1416         Q_tmp.AB_x[i] = -Q.AB_x[i];
1417         Q_tmp.AB_y[i] = -Q.AB_y[i];
1418         Q_tmp.AB_z[i] = -Q.AB_z[i];
1419     }
1420 
1421     int ret = ostei_i_h_k_p(P_tmp, Q_tmp, screen_tol, work, INT__h_i_p_k);
1422     double buffer[63504] SIMINT_ALIGN_ARRAY_DBL;
1423 
1424     for(int q = 0; q < ret; q++)
1425     {
1426         int idx = 0;
1427         for(int a = 0; a < 21; ++a)
1428         for(int b = 0; b < 28; ++b)
1429         for(int c = 0; c < 3; ++c)
1430         for(int d = 0; d < 36; ++d)
1431             buffer[idx++] = INT__h_i_p_k[q*63504+b*2268+a*108+d*3+c];
1432 
1433         memcpy(INT__h_i_p_k+q*63504, buffer, 63504*sizeof(double));
1434     }
1435 
1436     return ret;
1437 }
1438 
1439