1 #include "simint/boys/boys.h"
2 #include "simint/ostei/gen/ostei_generated.h"
3 #include "simint/vectorization/vectorization.h"
4 #include <math.h>
5 #include <string.h>
6 
7 
ostei_g_g_h_f(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__g_g_h_f)8 int ostei_g_g_h_f(struct simint_multi_shellpair const P,
9                   struct simint_multi_shellpair const Q,
10                   double screen_tol,
11                   double * const restrict work,
12                   double * const restrict INT__g_g_h_f)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__g_g_h_f);
17     int ab, cd, abcd;
18     int istart, jstart;
19     int iprimcd, nprim_icd, icd;
20     const int check_screen = (screen_tol > 0.0);
21     int i, j;
22     int n;
23     int not_screened;
24     int real_abcd;
25     int iket;
26     int ibra;
27 
28     // partition workspace
29     double * const INT__g_s_h_s = work + (SIMINT_NSHELL_SIMD * 0);
30     double * const INT__g_s_i_s = work + (SIMINT_NSHELL_SIMD * 315);
31     double * const INT__g_s_k_s = work + (SIMINT_NSHELL_SIMD * 735);
32     double * const INT__g_s_l_s = work + (SIMINT_NSHELL_SIMD * 1275);
33     double * const INT__h_s_h_s = work + (SIMINT_NSHELL_SIMD * 1950);
34     double * const INT__h_s_i_s = work + (SIMINT_NSHELL_SIMD * 2391);
35     double * const INT__h_s_k_s = work + (SIMINT_NSHELL_SIMD * 2979);
36     double * const INT__h_s_l_s = work + (SIMINT_NSHELL_SIMD * 3735);
37     double * const INT__i_s_h_s = work + (SIMINT_NSHELL_SIMD * 4680);
38     double * const INT__i_s_i_s = work + (SIMINT_NSHELL_SIMD * 5268);
39     double * const INT__i_s_k_s = work + (SIMINT_NSHELL_SIMD * 6052);
40     double * const INT__i_s_l_s = work + (SIMINT_NSHELL_SIMD * 7060);
41     double * const INT__k_s_h_s = work + (SIMINT_NSHELL_SIMD * 8320);
42     double * const INT__k_s_i_s = work + (SIMINT_NSHELL_SIMD * 9076);
43     double * const INT__k_s_k_s = work + (SIMINT_NSHELL_SIMD * 10084);
44     double * const INT__k_s_l_s = work + (SIMINT_NSHELL_SIMD * 11380);
45     double * const INT__l_s_h_s = work + (SIMINT_NSHELL_SIMD * 13000);
46     double * const INT__l_s_i_s = work + (SIMINT_NSHELL_SIMD * 13945);
47     double * const INT__l_s_k_s = work + (SIMINT_NSHELL_SIMD * 15205);
48     double * const INT__l_s_l_s = work + (SIMINT_NSHELL_SIMD * 16825);
49     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*18850);
50     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_p_s = primwork + 17;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_d_s = primwork + 41;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_f_s = primwork + 83;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_g_s = primwork + 143;
55     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_s = primwork + 218;
56     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_p_s = primwork + 266;
57     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_d_s = primwork + 338;
58     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_f_s = primwork + 464;
59     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_g_s = primwork + 644;
60     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_h_s = primwork + 869;
61     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 1121;
62     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_p_s = primwork + 1211;
63     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_d_s = primwork + 1355;
64     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_f_s = primwork + 1607;
65     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_g_s = primwork + 1967;
66     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_h_s = primwork + 2417;
67     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_i_s = primwork + 2921;
68     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_s = primwork + 3425;
69     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_p_s = primwork + 3565;
70     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_d_s = primwork + 3805;
71     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_f_s = primwork + 4225;
72     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_g_s = primwork + 4825;
73     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_h_s = primwork + 5575;
74     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_i_s = primwork + 6415;
75     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_k_s = primwork + 7255;
76     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_s_s = primwork + 7975;
77     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_p_s = primwork + 8170;
78     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_d_s = primwork + 8530;
79     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_f_s = primwork + 9160;
80     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_g_s = primwork + 10060;
81     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_h_s = primwork + 11185;
82     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_i_s = primwork + 12445;
83     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_k_s = primwork + 13705;
84     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_l_s = primwork + 14785;
85     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_s_s = primwork + 15460;
86     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_p_s = primwork + 15712;
87     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_d_s = primwork + 16216;
88     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_f_s = primwork + 17098;
89     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_g_s = primwork + 18358;
90     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_h_s = primwork + 19933;
91     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_i_s = primwork + 21697;
92     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_k_s = primwork + 23461;
93     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_l_s = primwork + 24973;
94     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_s_s = primwork + 25918;
95     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_p_s = primwork + 26226;
96     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_d_s = primwork + 26898;
97     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_f_s = primwork + 28074;
98     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_g_s = primwork + 29754;
99     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_h_s = primwork + 31854;
100     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_i_s = primwork + 34206;
101     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_k_s = primwork + 36558;
102     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_l_s = primwork + 38574;
103     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_s_s = primwork + 39834;
104     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_p_s = primwork + 40194;
105     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_d_s = primwork + 41058;
106     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_f_s = primwork + 42570;
107     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_g_s = primwork + 44730;
108     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_h_s = primwork + 47430;
109     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_i_s = primwork + 50454;
110     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_k_s = primwork + 53478;
111     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_l_s = primwork + 56070;
112     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_s_s = primwork + 57690;
113     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_p_s = primwork + 58095;
114     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_d_s = primwork + 59175;
115     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_f_s = primwork + 61065;
116     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_g_s = primwork + 63765;
117     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_h_s = primwork + 67140;
118     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_i_s = primwork + 70920;
119     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_k_s = primwork + 74700;
120     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_l_s = primwork + 77940;
121     double * const hrrwork = (double *)(primwork + 79965);
122     double * const HRR_INT__g_p_h_s = hrrwork + 0;
123     double * const HRR_INT__g_p_i_s = hrrwork + 945;
124     double * const HRR_INT__g_p_k_s = hrrwork + 2205;
125     double * const HRR_INT__g_p_l_s = hrrwork + 3825;
126     double * const HRR_INT__g_d_h_s = hrrwork + 5850;
127     double * const HRR_INT__g_d_i_s = hrrwork + 7740;
128     double * const HRR_INT__g_d_k_s = hrrwork + 10260;
129     double * const HRR_INT__g_d_l_s = hrrwork + 13500;
130     double * const HRR_INT__g_f_h_s = hrrwork + 17550;
131     double * const HRR_INT__g_f_i_s = hrrwork + 20700;
132     double * const HRR_INT__g_f_k_s = hrrwork + 24900;
133     double * const HRR_INT__g_f_l_s = hrrwork + 30300;
134     double * const HRR_INT__g_g_h_s = hrrwork + 37050;
135     double * const HRR_INT__g_g_h_p = hrrwork + 41775;
136     double * const HRR_INT__g_g_h_d = hrrwork + 55950;
137     double * const HRR_INT__g_g_i_s = hrrwork + 84300;
138     double * const HRR_INT__g_g_i_p = hrrwork + 90600;
139     double * const HRR_INT__g_g_i_d = hrrwork + 109500;
140     double * const HRR_INT__g_g_k_s = hrrwork + 147300;
141     double * const HRR_INT__g_g_k_p = hrrwork + 155400;
142     double * const HRR_INT__g_g_l_s = hrrwork + 179700;
143     double * const HRR_INT__h_p_h_s = hrrwork + 189825;
144     double * const HRR_INT__h_p_i_s = hrrwork + 191148;
145     double * const HRR_INT__h_p_k_s = hrrwork + 192912;
146     double * const HRR_INT__h_p_l_s = hrrwork + 195180;
147     double * const HRR_INT__h_d_h_s = hrrwork + 198015;
148     double * const HRR_INT__h_d_i_s = hrrwork + 200661;
149     double * const HRR_INT__h_d_k_s = hrrwork + 204189;
150     double * const HRR_INT__h_d_l_s = hrrwork + 208725;
151     double * const HRR_INT__h_f_h_s = hrrwork + 214395;
152     double * const HRR_INT__h_f_i_s = hrrwork + 218805;
153     double * const HRR_INT__h_f_k_s = hrrwork + 224685;
154     double * const HRR_INT__h_f_l_s = hrrwork + 232245;
155     double * const HRR_INT__i_p_h_s = hrrwork + 241695;
156     double * const HRR_INT__i_p_i_s = hrrwork + 243459;
157     double * const HRR_INT__i_p_k_s = hrrwork + 245811;
158     double * const HRR_INT__i_p_l_s = hrrwork + 248835;
159     double * const HRR_INT__i_d_h_s = hrrwork + 252615;
160     double * const HRR_INT__i_d_i_s = hrrwork + 256143;
161     double * const HRR_INT__i_d_k_s = hrrwork + 260847;
162     double * const HRR_INT__i_d_l_s = hrrwork + 266895;
163     double * const HRR_INT__k_p_h_s = hrrwork + 274455;
164     double * const HRR_INT__k_p_i_s = hrrwork + 276723;
165     double * const HRR_INT__k_p_k_s = hrrwork + 279747;
166     double * const HRR_INT__k_p_l_s = hrrwork + 283635;
167 
168 
169     // Create constants
170     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
171     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
172     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
173     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
174     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
175     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
176     const SIMINT_DBLTYPE const_7 = SIMINT_DBLSET1(7);
177     const SIMINT_DBLTYPE const_8 = SIMINT_DBLSET1(8);
178     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
179 
180 
181     ////////////////////////////////////////
182     // Loop over shells and primitives
183     ////////////////////////////////////////
184 
185     real_abcd = 0;
186     istart = 0;
187     for(ab = 0; ab < P.nshell12_clip; ++ab)
188     {
189         const int iend = istart + P.nprim12[ab];
190 
191         cd = 0;
192         jstart = 0;
193 
194         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
195         {
196             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
197             int jend = jstart;
198             for(i = 0; i < nshellbatch; i++)
199                 jend += Q.nprim12[cd+i];
200 
201             // Clear the beginning of the workspace (where we are accumulating integrals)
202             memset(work, 0, SIMINT_NSHELL_SIMD * 18850 * sizeof(double));
203             abcd = 0;
204 
205 
206             for(i = istart; i < iend; ++i)
207             {
208                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
209 
210                 if(check_screen)
211                 {
212                     // Skip this whole thing if always insignificant
213                     if((P.screen[i] * Q.screen_max) < screen_tol)
214                         continue;
215                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
216                 }
217 
218                 icd = 0;
219                 iprimcd = 0;
220                 nprim_icd = Q.nprim12[cd];
221                 double * restrict PRIM_PTR_INT__g_s_h_s = INT__g_s_h_s + abcd * 315;
222                 double * restrict PRIM_PTR_INT__g_s_i_s = INT__g_s_i_s + abcd * 420;
223                 double * restrict PRIM_PTR_INT__g_s_k_s = INT__g_s_k_s + abcd * 540;
224                 double * restrict PRIM_PTR_INT__g_s_l_s = INT__g_s_l_s + abcd * 675;
225                 double * restrict PRIM_PTR_INT__h_s_h_s = INT__h_s_h_s + abcd * 441;
226                 double * restrict PRIM_PTR_INT__h_s_i_s = INT__h_s_i_s + abcd * 588;
227                 double * restrict PRIM_PTR_INT__h_s_k_s = INT__h_s_k_s + abcd * 756;
228                 double * restrict PRIM_PTR_INT__h_s_l_s = INT__h_s_l_s + abcd * 945;
229                 double * restrict PRIM_PTR_INT__i_s_h_s = INT__i_s_h_s + abcd * 588;
230                 double * restrict PRIM_PTR_INT__i_s_i_s = INT__i_s_i_s + abcd * 784;
231                 double * restrict PRIM_PTR_INT__i_s_k_s = INT__i_s_k_s + abcd * 1008;
232                 double * restrict PRIM_PTR_INT__i_s_l_s = INT__i_s_l_s + abcd * 1260;
233                 double * restrict PRIM_PTR_INT__k_s_h_s = INT__k_s_h_s + abcd * 756;
234                 double * restrict PRIM_PTR_INT__k_s_i_s = INT__k_s_i_s + abcd * 1008;
235                 double * restrict PRIM_PTR_INT__k_s_k_s = INT__k_s_k_s + abcd * 1296;
236                 double * restrict PRIM_PTR_INT__k_s_l_s = INT__k_s_l_s + abcd * 1620;
237                 double * restrict PRIM_PTR_INT__l_s_h_s = INT__l_s_h_s + abcd * 945;
238                 double * restrict PRIM_PTR_INT__l_s_i_s = INT__l_s_i_s + abcd * 1260;
239                 double * restrict PRIM_PTR_INT__l_s_k_s = INT__l_s_k_s + abcd * 1620;
240                 double * restrict PRIM_PTR_INT__l_s_l_s = INT__l_s_l_s + abcd * 2025;
241 
242 
243 
244                 // Load these one per loop over i
245                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
246                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
247                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
248 
249                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
250 
251                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
252                 {
253                     // calculate the shell offsets
254                     // these are the offset from the shell pointed to by cd
255                     // for each element
256                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
257                     int lastoffset = 0;
258                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
259 
260                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
261                     {
262                         // Handle if the first element of the vector is a new shell
263                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
264                         {
265                             nprim_icd += Q.nprim12[cd + (++icd)];
266                             PRIM_PTR_INT__g_s_h_s += 315;
267                             PRIM_PTR_INT__g_s_i_s += 420;
268                             PRIM_PTR_INT__g_s_k_s += 540;
269                             PRIM_PTR_INT__g_s_l_s += 675;
270                             PRIM_PTR_INT__h_s_h_s += 441;
271                             PRIM_PTR_INT__h_s_i_s += 588;
272                             PRIM_PTR_INT__h_s_k_s += 756;
273                             PRIM_PTR_INT__h_s_l_s += 945;
274                             PRIM_PTR_INT__i_s_h_s += 588;
275                             PRIM_PTR_INT__i_s_i_s += 784;
276                             PRIM_PTR_INT__i_s_k_s += 1008;
277                             PRIM_PTR_INT__i_s_l_s += 1260;
278                             PRIM_PTR_INT__k_s_h_s += 756;
279                             PRIM_PTR_INT__k_s_i_s += 1008;
280                             PRIM_PTR_INT__k_s_k_s += 1296;
281                             PRIM_PTR_INT__k_s_l_s += 1620;
282                             PRIM_PTR_INT__l_s_h_s += 945;
283                             PRIM_PTR_INT__l_s_i_s += 1260;
284                             PRIM_PTR_INT__l_s_k_s += 1620;
285                             PRIM_PTR_INT__l_s_l_s += 2025;
286                         }
287                         iprimcd++;
288                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
289                         {
290                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
291                             {
292                                 shelloffsets[n] = shelloffsets[n-1] + 1;
293                                 lastoffset++;
294                                 nprim_icd += Q.nprim12[cd + (++icd)];
295                             }
296                             else
297                                 shelloffsets[n] = shelloffsets[n-1];
298                             iprimcd++;
299                         }
300                     }
301                     else
302                         iprimcd += SIMINT_SIMD_LEN;
303 
304                     // Do we have to compute this vector (or has it been screened out)?
305                     // (not_screened != 0 means we have to do this vector)
306                     if(check_screen)
307                     {
308                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
309                         if(vmax < screen_tol)
310                         {
311                             PRIM_PTR_INT__g_s_h_s += lastoffset*315;
312                             PRIM_PTR_INT__g_s_i_s += lastoffset*420;
313                             PRIM_PTR_INT__g_s_k_s += lastoffset*540;
314                             PRIM_PTR_INT__g_s_l_s += lastoffset*675;
315                             PRIM_PTR_INT__h_s_h_s += lastoffset*441;
316                             PRIM_PTR_INT__h_s_i_s += lastoffset*588;
317                             PRIM_PTR_INT__h_s_k_s += lastoffset*756;
318                             PRIM_PTR_INT__h_s_l_s += lastoffset*945;
319                             PRIM_PTR_INT__i_s_h_s += lastoffset*588;
320                             PRIM_PTR_INT__i_s_i_s += lastoffset*784;
321                             PRIM_PTR_INT__i_s_k_s += lastoffset*1008;
322                             PRIM_PTR_INT__i_s_l_s += lastoffset*1260;
323                             PRIM_PTR_INT__k_s_h_s += lastoffset*756;
324                             PRIM_PTR_INT__k_s_i_s += lastoffset*1008;
325                             PRIM_PTR_INT__k_s_k_s += lastoffset*1296;
326                             PRIM_PTR_INT__k_s_l_s += lastoffset*1620;
327                             PRIM_PTR_INT__l_s_h_s += lastoffset*945;
328                             PRIM_PTR_INT__l_s_i_s += lastoffset*1260;
329                             PRIM_PTR_INT__l_s_k_s += lastoffset*1620;
330                             PRIM_PTR_INT__l_s_l_s += lastoffset*2025;
331                             continue;
332                         }
333                     }
334 
335                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
336                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
337                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
338                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
339 
340 
341                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
342                     SIMINT_DBLTYPE PQ[3];
343                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
344                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
345                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
346                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
347                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
348                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
349 
350                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
351                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
352                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
353                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
354                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
355                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
356                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
357 
358                     // NOTE: Minus sign!
359                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
360                     SIMINT_DBLTYPE aop_PQ[3];
361                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
362                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
363                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
364 
365                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
366                     SIMINT_DBLTYPE aoq_PQ[3];
367                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
368                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
369                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
370                     // Put a minus sign here so we don't have to in RR routines
371                     a_over_q = SIMINT_NEG(a_over_q);
372 
373 
374                     //////////////////////////////////////////////
375                     // Fjt function section
376                     // Maximum v value: 16
377                     //////////////////////////////////////////////
378                     // The parameter to the Fjt function
379                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
380 
381 
382                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
383 
384 
385                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 16);
386                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
387                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
388                     for(n = 0; n <= 16; n++)
389                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
390 
391                     //////////////////////////////////////////////
392                     // Primitive integrals: Vertical recurrance
393                     //////////////////////////////////////////////
394 
395                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
396                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
397                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
398                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
399                     const SIMINT_DBLTYPE vrr_const_5_over_2p = SIMINT_MUL(const_5, one_over_2p);
400                     const SIMINT_DBLTYPE vrr_const_6_over_2p = SIMINT_MUL(const_6, one_over_2p);
401                     const SIMINT_DBLTYPE vrr_const_7_over_2p = SIMINT_MUL(const_7, one_over_2p);
402                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
403                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
404                     const SIMINT_DBLTYPE vrr_const_3_over_2q = SIMINT_MUL(const_3, one_over_2q);
405                     const SIMINT_DBLTYPE vrr_const_4_over_2q = SIMINT_MUL(const_4, one_over_2q);
406                     const SIMINT_DBLTYPE vrr_const_5_over_2q = SIMINT_MUL(const_5, one_over_2q);
407                     const SIMINT_DBLTYPE vrr_const_6_over_2q = SIMINT_MUL(const_6, one_over_2q);
408                     const SIMINT_DBLTYPE vrr_const_7_over_2q = SIMINT_MUL(const_7, one_over_2q);
409                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
410                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
411                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
412                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
413                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
414                     const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
415                     const SIMINT_DBLTYPE vrr_const_7_over_2pq = SIMINT_MUL(const_7, one_over_2pq);
416                     const SIMINT_DBLTYPE vrr_const_8_over_2pq = SIMINT_MUL(const_8, one_over_2pq);
417 
418 
419 
420                     // Forming PRIM_INT__p_s_s_s[16 * 3];
421                     for(n = 0; n < 16; ++n)  // loop over orders of auxiliary function
422                     {
423 
424                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
425                         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]);
426 
427                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
428                         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]);
429 
430                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
431                         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]);
432 
433                     }
434 
435 
436 
437                     // Forming PRIM_INT__d_s_s_s[15 * 6];
438                     for(n = 0; n < 15; ++n)  // loop over orders of auxiliary function
439                     {
440 
441                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
442                         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]);
443                         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]);
444 
445                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
446                         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]);
447 
448                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
449                         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]);
450 
451                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
452                         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]);
453                         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]);
454 
455                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
456                         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]);
457 
458                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
459                         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]);
460                         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]);
461 
462                     }
463 
464 
465 
466                     // Forming PRIM_INT__f_s_s_s[14 * 10];
467                     for(n = 0; n < 14; ++n)  // loop over orders of auxiliary function
468                     {
469 
470                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
471                         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]);
472                         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]);
473 
474                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
475                         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]);
476 
477                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
478                         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]);
479 
480                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
481                         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]);
482 
483                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
484                         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]);
485 
486                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
487                         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]);
488 
489                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
490                         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]);
491                         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]);
492 
493                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
494                         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]);
495 
496                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
497                         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]);
498 
499                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
500                         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]);
501                         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]);
502 
503                     }
504 
505 
506                     VRR_I_g_s_s_s(
507                             PRIM_INT__g_s_s_s,
508                             PRIM_INT__f_s_s_s,
509                             PRIM_INT__d_s_s_s,
510                             P_PA,
511                             a_over_p,
512                             aop_PQ,
513                             one_over_2p,
514                             13);
515 
516 
517                     VRR_K_g_s_p_s(
518                             PRIM_INT__g_s_p_s,
519                             PRIM_INT__g_s_s_s,
520                             PRIM_INT__f_s_s_s,
521                             Q_PA,
522                             aoq_PQ,
523                             one_over_2pq,
524                             8);
525 
526 
527                     VRR_K_f_s_p_s(
528                             PRIM_INT__f_s_p_s,
529                             PRIM_INT__f_s_s_s,
530                             PRIM_INT__d_s_s_s,
531                             Q_PA,
532                             aoq_PQ,
533                             one_over_2pq,
534                             8);
535 
536 
537                     ostei_general_vrr_K(4, 0, 2, 0, 7,
538                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
539                             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);
540 
541 
542 
543                     // Forming PRIM_INT__d_s_p_s[8 * 18];
544                     for(n = 0; n < 8; ++n)  // loop over orders of auxiliary function
545                     {
546 
547                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
548                         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]);
549                         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]);
550 
551                         PRIM_INT__d_s_p_s[n * 18 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
552                         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]);
553 
554                         PRIM_INT__d_s_p_s[n * 18 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
555                         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]);
556 
557                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 1]);
558                         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]);
559                         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]);
560 
561                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 1]);
562                         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]);
563                         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]);
564 
565                         PRIM_INT__d_s_p_s[n * 18 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
566                         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]);
567 
568                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 2]);
569                         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]);
570                         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]);
571 
572                         PRIM_INT__d_s_p_s[n * 18 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 2]);
573                         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]);
574 
575                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 2]);
576                         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]);
577                         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]);
578 
579                         PRIM_INT__d_s_p_s[n * 18 + 9] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
580                         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]);
581 
582                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
583                         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]);
584                         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]);
585 
586                         PRIM_INT__d_s_p_s[n * 18 + 11] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
587                         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]);
588 
589                         PRIM_INT__d_s_p_s[n * 18 + 12] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 4]);
590                         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]);
591 
592                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 4]);
593                         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]);
594                         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]);
595 
596                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 4]);
597                         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]);
598                         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]);
599 
600                         PRIM_INT__d_s_p_s[n * 18 + 15] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
601                         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]);
602 
603                         PRIM_INT__d_s_p_s[n * 18 + 16] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
604                         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]);
605 
606                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
607                         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]);
608                         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]);
609 
610                     }
611 
612 
613                     VRR_K_f_s_d_s(
614                             PRIM_INT__f_s_d_s,
615                             PRIM_INT__f_s_p_s,
616                             PRIM_INT__f_s_s_s,
617                             PRIM_INT__d_s_p_s,
618                             Q_PA,
619                             a_over_q,
620                             aoq_PQ,
621                             one_over_2pq,
622                             one_over_2q,
623                             7);
624 
625 
626                     ostei_general_vrr_K(4, 0, 3, 0, 6,
627                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
628                             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);
629 
630 
631 
632                     // Forming PRIM_INT__p_s_p_s[8 * 9];
633                     for(n = 0; n < 8; ++n)  // loop over orders of auxiliary function
634                     {
635 
636                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
637                         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]);
638                         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]);
639 
640                         PRIM_INT__p_s_p_s[n * 9 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
641                         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]);
642 
643                         PRIM_INT__p_s_p_s[n * 9 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
644                         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]);
645 
646                         PRIM_INT__p_s_p_s[n * 9 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 1]);
647                         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]);
648 
649                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
650                         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]);
651                         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]);
652 
653                         PRIM_INT__p_s_p_s[n * 9 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
654                         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]);
655 
656                         PRIM_INT__p_s_p_s[n * 9 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 2]);
657                         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]);
658 
659                         PRIM_INT__p_s_p_s[n * 9 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 2]);
660                         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]);
661 
662                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
663                         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]);
664                         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]);
665 
666                     }
667 
668 
669                     VRR_K_d_s_d_s(
670                             PRIM_INT__d_s_d_s,
671                             PRIM_INT__d_s_p_s,
672                             PRIM_INT__d_s_s_s,
673                             PRIM_INT__p_s_p_s,
674                             Q_PA,
675                             a_over_q,
676                             aoq_PQ,
677                             one_over_2pq,
678                             one_over_2q,
679                             7);
680 
681 
682                     ostei_general_vrr_K(3, 0, 3, 0, 6,
683                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
684                             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);
685 
686 
687                     ostei_general_vrr_K(4, 0, 4, 0, 5,
688                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
689                             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);
690 
691 
692 
693                     // Forming PRIM_INT__s_s_p_s[8 * 3];
694                     for(n = 0; n < 8; ++n)  // loop over orders of auxiliary function
695                     {
696 
697                         PRIM_INT__s_s_p_s[n * 3 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
698                         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]);
699 
700                         PRIM_INT__s_s_p_s[n * 3 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
701                         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]);
702 
703                         PRIM_INT__s_s_p_s[n * 3 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
704                         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]);
705 
706                     }
707 
708 
709 
710                     // Forming PRIM_INT__p_s_d_s[7 * 18];
711                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
712                     {
713 
714                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_p_s[n * 9 + 0]);
715                         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]);
716                         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]);
717                         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]);
718 
719                         PRIM_INT__p_s_d_s[n * 18 + 3] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_p_s[n * 9 + 1]);
720                         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]);
721                         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]);
722 
723                         PRIM_INT__p_s_d_s[n * 18 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_p_s[n * 9 + 2]);
724                         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]);
725                         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]);
726 
727                         PRIM_INT__p_s_d_s[n * 18 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_p_s[n * 9 + 3]);
728                         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]);
729                         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]);
730 
731                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_p_s[n * 9 + 4]);
732                         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]);
733                         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]);
734                         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]);
735 
736                         PRIM_INT__p_s_d_s[n * 18 + 11] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_p_s[n * 9 + 5]);
737                         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]);
738                         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]);
739 
740                         PRIM_INT__p_s_d_s[n * 18 + 12] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_p_s[n * 9 + 6]);
741                         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]);
742                         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]);
743 
744                         PRIM_INT__p_s_d_s[n * 18 + 15] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_p_s[n * 9 + 7]);
745                         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]);
746                         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]);
747 
748                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_p_s[n * 9 + 8]);
749                         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]);
750                         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]);
751                         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]);
752 
753                     }
754 
755 
756                     VRR_K_d_s_f_s(
757                             PRIM_INT__d_s_f_s,
758                             PRIM_INT__d_s_d_s,
759                             PRIM_INT__d_s_p_s,
760                             PRIM_INT__p_s_d_s,
761                             Q_PA,
762                             a_over_q,
763                             aoq_PQ,
764                             one_over_2pq,
765                             one_over_2q,
766                             6);
767 
768 
769                     ostei_general_vrr_K(3, 0, 4, 0, 5,
770                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
771                             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);
772 
773 
774                     ostei_general_vrr_K(4, 0, 5, 0, 4,
775                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
776                             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);
777 
778 
779 
780                     // Forming PRIM_INT__s_s_d_s[7 * 6];
781                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
782                     {
783 
784                         PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_p_s[n * 3 + 0]);
785                         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]);
786                         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]);
787 
788                         PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_p_s[n * 3 + 1]);
789                         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]);
790                         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]);
791 
792                         PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 2]);
793                         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]);
794                         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]);
795 
796                     }
797 
798 
799                     VRR_K_p_s_f_s(
800                             PRIM_INT__p_s_f_s,
801                             PRIM_INT__p_s_d_s,
802                             PRIM_INT__p_s_p_s,
803                             PRIM_INT__s_s_d_s,
804                             Q_PA,
805                             a_over_q,
806                             aoq_PQ,
807                             one_over_2pq,
808                             one_over_2q,
809                             6);
810 
811 
812                     ostei_general_vrr_K(2, 0, 4, 0, 5,
813                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
814                             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);
815 
816 
817                     ostei_general_vrr_K(3, 0, 5, 0, 4,
818                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
819                             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);
820 
821 
822                     ostei_general_vrr_K(4, 0, 6, 0, 3,
823                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
824                             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);
825 
826 
827                     VRR_I_h_s_s_s(
828                             PRIM_INT__h_s_s_s,
829                             PRIM_INT__g_s_s_s,
830                             PRIM_INT__f_s_s_s,
831                             P_PA,
832                             a_over_p,
833                             aop_PQ,
834                             one_over_2p,
835                             12);
836 
837 
838                     ostei_general_vrr_K(5, 0, 1, 0, 8,
839                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
840                             PRIM_INT__h_s_s_s, NULL, NULL, PRIM_INT__g_s_s_s, NULL, PRIM_INT__h_s_p_s);
841 
842 
843                     ostei_general_vrr_K(5, 0, 2, 0, 7,
844                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
845                             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);
846 
847 
848                     ostei_general_vrr_K(5, 0, 3, 0, 6,
849                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
850                             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);
851 
852 
853                     ostei_general_vrr_K(5, 0, 4, 0, 5,
854                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
855                             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);
856 
857 
858                     ostei_general_vrr_K(5, 0, 5, 0, 4,
859                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
860                             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);
861 
862 
863 
864                     // Forming PRIM_INT__s_s_f_s[6 * 10];
865                     for(n = 0; n < 6; ++n)  // loop over orders of auxiliary function
866                     {
867 
868                         PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 0]);
869                         PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__s_s_f_s[n * 10 + 0]);
870                         PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_p_s[n * 3 + 0]), PRIM_INT__s_s_f_s[n * 10 + 0]);
871 
872                         PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 3]);
873                         PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__s_s_f_s[n * 10 + 6]);
874                         PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_p_s[n * 3 + 1]), PRIM_INT__s_s_f_s[n * 10 + 6]);
875 
876                         PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 5]);
877                         PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__s_s_f_s[n * 10 + 9]);
878                         PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__s_s_p_s[n * 3 + 2]), PRIM_INT__s_s_f_s[n * 10 + 9]);
879 
880                     }
881 
882 
883                     VRR_K_p_s_g_s(
884                             PRIM_INT__p_s_g_s,
885                             PRIM_INT__p_s_f_s,
886                             PRIM_INT__p_s_d_s,
887                             PRIM_INT__s_s_f_s,
888                             Q_PA,
889                             a_over_q,
890                             aoq_PQ,
891                             one_over_2pq,
892                             one_over_2q,
893                             5);
894 
895 
896                     ostei_general_vrr_K(2, 0, 5, 0, 4,
897                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
898                             PRIM_INT__d_s_g_s, PRIM_INT__d_s_f_s, NULL, PRIM_INT__p_s_g_s, NULL, PRIM_INT__d_s_h_s);
899 
900 
901                     ostei_general_vrr_K(3, 0, 6, 0, 3,
902                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
903                             PRIM_INT__f_s_h_s, PRIM_INT__f_s_g_s, NULL, PRIM_INT__d_s_h_s, NULL, PRIM_INT__f_s_i_s);
904 
905 
906                     ostei_general_vrr_K(4, 0, 7, 0, 2,
907                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
908                             PRIM_INT__g_s_i_s, PRIM_INT__g_s_h_s, NULL, PRIM_INT__f_s_i_s, NULL, PRIM_INT__g_s_k_s);
909 
910 
911                     ostei_general_vrr_K(5, 0, 6, 0, 3,
912                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
913                             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);
914 
915 
916                     ostei_general_vrr1_I(6, 11,
917                             one_over_2p, a_over_p, aop_PQ, P_PA,
918                             PRIM_INT__h_s_s_s, PRIM_INT__g_s_s_s, PRIM_INT__i_s_s_s);
919 
920 
921                     ostei_general_vrr_K(6, 0, 1, 0, 8,
922                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
923                             PRIM_INT__i_s_s_s, NULL, NULL, PRIM_INT__h_s_s_s, NULL, PRIM_INT__i_s_p_s);
924 
925 
926                     ostei_general_vrr_K(6, 0, 2, 0, 7,
927                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
928                             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);
929 
930 
931                     ostei_general_vrr_K(6, 0, 3, 0, 6,
932                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
933                             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);
934 
935 
936                     ostei_general_vrr_K(6, 0, 4, 0, 5,
937                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
938                             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);
939 
940 
941                     ostei_general_vrr_K(6, 0, 5, 0, 4,
942                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
943                             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);
944 
945 
946                     VRR_K_s_s_g_s(
947                             PRIM_INT__s_s_g_s,
948                             PRIM_INT__s_s_f_s,
949                             PRIM_INT__s_s_d_s,
950                             Q_PA,
951                             a_over_q,
952                             aoq_PQ,
953                             one_over_2q,
954                             5);
955 
956 
957                     ostei_general_vrr_K(1, 0, 5, 0, 4,
958                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
959                             PRIM_INT__p_s_g_s, PRIM_INT__p_s_f_s, NULL, PRIM_INT__s_s_g_s, NULL, PRIM_INT__p_s_h_s);
960 
961 
962                     ostei_general_vrr_K(2, 0, 6, 0, 3,
963                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
964                             PRIM_INT__d_s_h_s, PRIM_INT__d_s_g_s, NULL, PRIM_INT__p_s_h_s, NULL, PRIM_INT__d_s_i_s);
965 
966 
967                     ostei_general_vrr_K(3, 0, 7, 0, 2,
968                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
969                             PRIM_INT__f_s_i_s, PRIM_INT__f_s_h_s, NULL, PRIM_INT__d_s_i_s, NULL, PRIM_INT__f_s_k_s);
970 
971 
972                     ostei_general_vrr_K(4, 0, 8, 0, 1,
973                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
974                             PRIM_INT__g_s_k_s, PRIM_INT__g_s_i_s, NULL, PRIM_INT__f_s_k_s, NULL, PRIM_INT__g_s_l_s);
975 
976 
977                     ostei_general_vrr_K(5, 0, 7, 0, 2,
978                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
979                             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);
980 
981 
982                     ostei_general_vrr_K(6, 0, 6, 0, 3,
983                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
984                             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);
985 
986 
987                     ostei_general_vrr1_I(7, 10,
988                             one_over_2p, a_over_p, aop_PQ, P_PA,
989                             PRIM_INT__i_s_s_s, PRIM_INT__h_s_s_s, PRIM_INT__k_s_s_s);
990 
991 
992                     ostei_general_vrr_K(7, 0, 1, 0, 8,
993                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
994                             PRIM_INT__k_s_s_s, NULL, NULL, PRIM_INT__i_s_s_s, NULL, PRIM_INT__k_s_p_s);
995 
996 
997                     ostei_general_vrr_K(7, 0, 2, 0, 7,
998                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
999                             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);
1000 
1001 
1002                     ostei_general_vrr_K(7, 0, 3, 0, 6,
1003                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1004                             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);
1005 
1006 
1007                     ostei_general_vrr_K(7, 0, 4, 0, 5,
1008                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1009                             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);
1010 
1011 
1012                     ostei_general_vrr_K(7, 0, 5, 0, 4,
1013                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1014                             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);
1015 
1016 
1017                     ostei_general_vrr_K(5, 0, 8, 0, 1,
1018                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1019                             PRIM_INT__h_s_k_s, PRIM_INT__h_s_i_s, NULL, PRIM_INT__g_s_k_s, NULL, PRIM_INT__h_s_l_s);
1020 
1021 
1022                     ostei_general_vrr_K(6, 0, 7, 0, 2,
1023                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1024                             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);
1025 
1026 
1027                     ostei_general_vrr_K(7, 0, 6, 0, 3,
1028                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1029                             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);
1030 
1031 
1032                     ostei_general_vrr1_I(8, 9,
1033                             one_over_2p, a_over_p, aop_PQ, P_PA,
1034                             PRIM_INT__k_s_s_s, PRIM_INT__i_s_s_s, PRIM_INT__l_s_s_s);
1035 
1036 
1037                     ostei_general_vrr_K(8, 0, 1, 0, 8,
1038                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1039                             PRIM_INT__l_s_s_s, NULL, NULL, PRIM_INT__k_s_s_s, NULL, PRIM_INT__l_s_p_s);
1040 
1041 
1042                     ostei_general_vrr_K(8, 0, 2, 0, 7,
1043                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1044                             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);
1045 
1046 
1047                     ostei_general_vrr_K(8, 0, 3, 0, 6,
1048                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1049                             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);
1050 
1051 
1052                     ostei_general_vrr_K(8, 0, 4, 0, 5,
1053                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1054                             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);
1055 
1056 
1057                     ostei_general_vrr_K(8, 0, 5, 0, 4,
1058                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1059                             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);
1060 
1061 
1062                     ostei_general_vrr_K(6, 0, 8, 0, 1,
1063                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1064                             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);
1065 
1066 
1067                     ostei_general_vrr_K(7, 0, 7, 0, 2,
1068                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1069                             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);
1070 
1071 
1072                     ostei_general_vrr_K(8, 0, 6, 0, 3,
1073                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1074                             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);
1075 
1076 
1077                     ostei_general_vrr_K(7, 0, 8, 0, 1,
1078                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1079                             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);
1080 
1081 
1082                     ostei_general_vrr_K(8, 0, 7, 0, 2,
1083                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1084                             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);
1085 
1086 
1087                     ostei_general_vrr_K(8, 0, 8, 0, 1,
1088                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
1089                             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);
1090 
1091 
1092 
1093 
1094                     ////////////////////////////////////
1095                     // Accumulate contracted integrals
1096                     ////////////////////////////////////
1097                     if(lastoffset == 0)
1098                     {
1099                         contract_all(315, PRIM_INT__g_s_h_s, PRIM_PTR_INT__g_s_h_s);
1100                         contract_all(420, PRIM_INT__g_s_i_s, PRIM_PTR_INT__g_s_i_s);
1101                         contract_all(540, PRIM_INT__g_s_k_s, PRIM_PTR_INT__g_s_k_s);
1102                         contract_all(675, PRIM_INT__g_s_l_s, PRIM_PTR_INT__g_s_l_s);
1103                         contract_all(441, PRIM_INT__h_s_h_s, PRIM_PTR_INT__h_s_h_s);
1104                         contract_all(588, PRIM_INT__h_s_i_s, PRIM_PTR_INT__h_s_i_s);
1105                         contract_all(756, PRIM_INT__h_s_k_s, PRIM_PTR_INT__h_s_k_s);
1106                         contract_all(945, PRIM_INT__h_s_l_s, PRIM_PTR_INT__h_s_l_s);
1107                         contract_all(588, PRIM_INT__i_s_h_s, PRIM_PTR_INT__i_s_h_s);
1108                         contract_all(784, PRIM_INT__i_s_i_s, PRIM_PTR_INT__i_s_i_s);
1109                         contract_all(1008, PRIM_INT__i_s_k_s, PRIM_PTR_INT__i_s_k_s);
1110                         contract_all(1260, PRIM_INT__i_s_l_s, PRIM_PTR_INT__i_s_l_s);
1111                         contract_all(756, PRIM_INT__k_s_h_s, PRIM_PTR_INT__k_s_h_s);
1112                         contract_all(1008, PRIM_INT__k_s_i_s, PRIM_PTR_INT__k_s_i_s);
1113                         contract_all(1296, PRIM_INT__k_s_k_s, PRIM_PTR_INT__k_s_k_s);
1114                         contract_all(1620, PRIM_INT__k_s_l_s, PRIM_PTR_INT__k_s_l_s);
1115                         contract_all(945, PRIM_INT__l_s_h_s, PRIM_PTR_INT__l_s_h_s);
1116                         contract_all(1260, PRIM_INT__l_s_i_s, PRIM_PTR_INT__l_s_i_s);
1117                         contract_all(1620, PRIM_INT__l_s_k_s, PRIM_PTR_INT__l_s_k_s);
1118                         contract_all(2025, PRIM_INT__l_s_l_s, PRIM_PTR_INT__l_s_l_s);
1119                     }
1120                     else
1121                     {
1122                         contract(315, shelloffsets, PRIM_INT__g_s_h_s, PRIM_PTR_INT__g_s_h_s);
1123                         contract(420, shelloffsets, PRIM_INT__g_s_i_s, PRIM_PTR_INT__g_s_i_s);
1124                         contract(540, shelloffsets, PRIM_INT__g_s_k_s, PRIM_PTR_INT__g_s_k_s);
1125                         contract(675, shelloffsets, PRIM_INT__g_s_l_s, PRIM_PTR_INT__g_s_l_s);
1126                         contract(441, shelloffsets, PRIM_INT__h_s_h_s, PRIM_PTR_INT__h_s_h_s);
1127                         contract(588, shelloffsets, PRIM_INT__h_s_i_s, PRIM_PTR_INT__h_s_i_s);
1128                         contract(756, shelloffsets, PRIM_INT__h_s_k_s, PRIM_PTR_INT__h_s_k_s);
1129                         contract(945, shelloffsets, PRIM_INT__h_s_l_s, PRIM_PTR_INT__h_s_l_s);
1130                         contract(588, shelloffsets, PRIM_INT__i_s_h_s, PRIM_PTR_INT__i_s_h_s);
1131                         contract(784, shelloffsets, PRIM_INT__i_s_i_s, PRIM_PTR_INT__i_s_i_s);
1132                         contract(1008, shelloffsets, PRIM_INT__i_s_k_s, PRIM_PTR_INT__i_s_k_s);
1133                         contract(1260, shelloffsets, PRIM_INT__i_s_l_s, PRIM_PTR_INT__i_s_l_s);
1134                         contract(756, shelloffsets, PRIM_INT__k_s_h_s, PRIM_PTR_INT__k_s_h_s);
1135                         contract(1008, shelloffsets, PRIM_INT__k_s_i_s, PRIM_PTR_INT__k_s_i_s);
1136                         contract(1296, shelloffsets, PRIM_INT__k_s_k_s, PRIM_PTR_INT__k_s_k_s);
1137                         contract(1620, shelloffsets, PRIM_INT__k_s_l_s, PRIM_PTR_INT__k_s_l_s);
1138                         contract(945, shelloffsets, PRIM_INT__l_s_h_s, PRIM_PTR_INT__l_s_h_s);
1139                         contract(1260, shelloffsets, PRIM_INT__l_s_i_s, PRIM_PTR_INT__l_s_i_s);
1140                         contract(1620, shelloffsets, PRIM_INT__l_s_k_s, PRIM_PTR_INT__l_s_k_s);
1141                         contract(2025, shelloffsets, PRIM_INT__l_s_l_s, PRIM_PTR_INT__l_s_l_s);
1142                         PRIM_PTR_INT__g_s_h_s += lastoffset*315;
1143                         PRIM_PTR_INT__g_s_i_s += lastoffset*420;
1144                         PRIM_PTR_INT__g_s_k_s += lastoffset*540;
1145                         PRIM_PTR_INT__g_s_l_s += lastoffset*675;
1146                         PRIM_PTR_INT__h_s_h_s += lastoffset*441;
1147                         PRIM_PTR_INT__h_s_i_s += lastoffset*588;
1148                         PRIM_PTR_INT__h_s_k_s += lastoffset*756;
1149                         PRIM_PTR_INT__h_s_l_s += lastoffset*945;
1150                         PRIM_PTR_INT__i_s_h_s += lastoffset*588;
1151                         PRIM_PTR_INT__i_s_i_s += lastoffset*784;
1152                         PRIM_PTR_INT__i_s_k_s += lastoffset*1008;
1153                         PRIM_PTR_INT__i_s_l_s += lastoffset*1260;
1154                         PRIM_PTR_INT__k_s_h_s += lastoffset*756;
1155                         PRIM_PTR_INT__k_s_i_s += lastoffset*1008;
1156                         PRIM_PTR_INT__k_s_k_s += lastoffset*1296;
1157                         PRIM_PTR_INT__k_s_l_s += lastoffset*1620;
1158                         PRIM_PTR_INT__l_s_h_s += lastoffset*945;
1159                         PRIM_PTR_INT__l_s_i_s += lastoffset*1260;
1160                         PRIM_PTR_INT__l_s_k_s += lastoffset*1620;
1161                         PRIM_PTR_INT__l_s_l_s += lastoffset*2025;
1162                     }
1163 
1164                 }  // close loop over j
1165             }  // close loop over i
1166 
1167             //Advance to the next batch
1168             jstart = SIMINT_SIMD_ROUND(jend);
1169 
1170             //////////////////////////////////////////////
1171             // Contracted integrals: Horizontal recurrance
1172             //////////////////////////////////////////////
1173 
1174 
1175             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
1176 
1177 
1178             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
1179             {
1180                 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
1181 
1182                 // set up HRR pointers
1183                 double const * restrict HRR_INT__g_s_h_s = INT__g_s_h_s + abcd * 315;
1184                 double const * restrict HRR_INT__g_s_i_s = INT__g_s_i_s + abcd * 420;
1185                 double const * restrict HRR_INT__g_s_k_s = INT__g_s_k_s + abcd * 540;
1186                 double const * restrict HRR_INT__g_s_l_s = INT__g_s_l_s + abcd * 675;
1187                 double const * restrict HRR_INT__h_s_h_s = INT__h_s_h_s + abcd * 441;
1188                 double const * restrict HRR_INT__h_s_i_s = INT__h_s_i_s + abcd * 588;
1189                 double const * restrict HRR_INT__h_s_k_s = INT__h_s_k_s + abcd * 756;
1190                 double const * restrict HRR_INT__h_s_l_s = INT__h_s_l_s + abcd * 945;
1191                 double const * restrict HRR_INT__i_s_h_s = INT__i_s_h_s + abcd * 588;
1192                 double const * restrict HRR_INT__i_s_i_s = INT__i_s_i_s + abcd * 784;
1193                 double const * restrict HRR_INT__i_s_k_s = INT__i_s_k_s + abcd * 1008;
1194                 double const * restrict HRR_INT__i_s_l_s = INT__i_s_l_s + abcd * 1260;
1195                 double const * restrict HRR_INT__k_s_h_s = INT__k_s_h_s + abcd * 756;
1196                 double const * restrict HRR_INT__k_s_i_s = INT__k_s_i_s + abcd * 1008;
1197                 double const * restrict HRR_INT__k_s_k_s = INT__k_s_k_s + abcd * 1296;
1198                 double const * restrict HRR_INT__k_s_l_s = INT__k_s_l_s + abcd * 1620;
1199                 double const * restrict HRR_INT__l_s_h_s = INT__l_s_h_s + abcd * 945;
1200                 double const * restrict HRR_INT__l_s_i_s = INT__l_s_i_s + abcd * 1260;
1201                 double const * restrict HRR_INT__l_s_k_s = INT__l_s_k_s + abcd * 1620;
1202                 double const * restrict HRR_INT__l_s_l_s = INT__l_s_l_s + abcd * 2025;
1203                 double * restrict HRR_INT__g_g_h_f = INT__g_g_h_f + real_abcd * 47250;
1204 
1205                 // form INT__g_p_h_s
1206                 HRR_J_g_p(
1207                     HRR_INT__g_p_h_s,
1208                     HRR_INT__g_s_h_s,
1209                     HRR_INT__h_s_h_s,
1210                     hAB, 21);
1211 
1212                 // form INT__g_p_i_s
1213                 HRR_J_g_p(
1214                     HRR_INT__g_p_i_s,
1215                     HRR_INT__g_s_i_s,
1216                     HRR_INT__h_s_i_s,
1217                     hAB, 28);
1218 
1219                 // form INT__g_p_k_s
1220                 HRR_J_g_p(
1221                     HRR_INT__g_p_k_s,
1222                     HRR_INT__g_s_k_s,
1223                     HRR_INT__h_s_k_s,
1224                     hAB, 36);
1225 
1226                 // form INT__g_p_l_s
1227                 HRR_J_g_p(
1228                     HRR_INT__g_p_l_s,
1229                     HRR_INT__g_s_l_s,
1230                     HRR_INT__h_s_l_s,
1231                     hAB, 45);
1232 
1233                 // form INT__h_p_h_s
1234                 ostei_general_hrr_J(5, 1, 5, 0, hAB, HRR_INT__i_s_h_s, HRR_INT__h_s_h_s, HRR_INT__h_p_h_s);
1235 
1236                 // form INT__h_p_i_s
1237                 ostei_general_hrr_J(5, 1, 6, 0, hAB, HRR_INT__i_s_i_s, HRR_INT__h_s_i_s, HRR_INT__h_p_i_s);
1238 
1239                 // form INT__h_p_k_s
1240                 ostei_general_hrr_J(5, 1, 7, 0, hAB, HRR_INT__i_s_k_s, HRR_INT__h_s_k_s, HRR_INT__h_p_k_s);
1241 
1242                 // form INT__h_p_l_s
1243                 ostei_general_hrr_J(5, 1, 8, 0, hAB, HRR_INT__i_s_l_s, HRR_INT__h_s_l_s, HRR_INT__h_p_l_s);
1244 
1245                 // form INT__i_p_h_s
1246                 ostei_general_hrr_J(6, 1, 5, 0, hAB, HRR_INT__k_s_h_s, HRR_INT__i_s_h_s, HRR_INT__i_p_h_s);
1247 
1248                 // form INT__i_p_i_s
1249                 ostei_general_hrr_J(6, 1, 6, 0, hAB, HRR_INT__k_s_i_s, HRR_INT__i_s_i_s, HRR_INT__i_p_i_s);
1250 
1251                 // form INT__i_p_k_s
1252                 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);
1253 
1254                 // form INT__i_p_l_s
1255                 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);
1256 
1257                 // form INT__k_p_h_s
1258                 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);
1259 
1260                 // form INT__k_p_i_s
1261                 ostei_general_hrr_J(7, 1, 6, 0, hAB, HRR_INT__l_s_i_s, HRR_INT__k_s_i_s, HRR_INT__k_p_i_s);
1262 
1263                 // form INT__k_p_k_s
1264                 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);
1265 
1266                 // form INT__k_p_l_s
1267                 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);
1268 
1269                 // form INT__g_d_h_s
1270                 ostei_general_hrr_J(4, 2, 5, 0, hAB, HRR_INT__h_p_h_s, HRR_INT__g_p_h_s, HRR_INT__g_d_h_s);
1271 
1272                 // form INT__g_d_i_s
1273                 ostei_general_hrr_J(4, 2, 6, 0, hAB, HRR_INT__h_p_i_s, HRR_INT__g_p_i_s, HRR_INT__g_d_i_s);
1274 
1275                 // form INT__g_d_k_s
1276                 ostei_general_hrr_J(4, 2, 7, 0, hAB, HRR_INT__h_p_k_s, HRR_INT__g_p_k_s, HRR_INT__g_d_k_s);
1277 
1278                 // form INT__g_d_l_s
1279                 ostei_general_hrr_J(4, 2, 8, 0, hAB, HRR_INT__h_p_l_s, HRR_INT__g_p_l_s, HRR_INT__g_d_l_s);
1280 
1281                 // form INT__h_d_h_s
1282                 ostei_general_hrr_J(5, 2, 5, 0, hAB, HRR_INT__i_p_h_s, HRR_INT__h_p_h_s, HRR_INT__h_d_h_s);
1283 
1284                 // form INT__h_d_i_s
1285                 ostei_general_hrr_J(5, 2, 6, 0, hAB, HRR_INT__i_p_i_s, HRR_INT__h_p_i_s, HRR_INT__h_d_i_s);
1286 
1287                 // form INT__h_d_k_s
1288                 ostei_general_hrr_J(5, 2, 7, 0, hAB, HRR_INT__i_p_k_s, HRR_INT__h_p_k_s, HRR_INT__h_d_k_s);
1289 
1290                 // form INT__h_d_l_s
1291                 ostei_general_hrr_J(5, 2, 8, 0, hAB, HRR_INT__i_p_l_s, HRR_INT__h_p_l_s, HRR_INT__h_d_l_s);
1292 
1293                 // form INT__i_d_h_s
1294                 ostei_general_hrr_J(6, 2, 5, 0, hAB, HRR_INT__k_p_h_s, HRR_INT__i_p_h_s, HRR_INT__i_d_h_s);
1295 
1296                 // form INT__i_d_i_s
1297                 ostei_general_hrr_J(6, 2, 6, 0, hAB, HRR_INT__k_p_i_s, HRR_INT__i_p_i_s, HRR_INT__i_d_i_s);
1298 
1299                 // form INT__i_d_k_s
1300                 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);
1301 
1302                 // form INT__i_d_l_s
1303                 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);
1304 
1305                 // form INT__g_f_h_s
1306                 ostei_general_hrr_J(4, 3, 5, 0, hAB, HRR_INT__h_d_h_s, HRR_INT__g_d_h_s, HRR_INT__g_f_h_s);
1307 
1308                 // form INT__g_f_i_s
1309                 ostei_general_hrr_J(4, 3, 6, 0, hAB, HRR_INT__h_d_i_s, HRR_INT__g_d_i_s, HRR_INT__g_f_i_s);
1310 
1311                 // form INT__g_f_k_s
1312                 ostei_general_hrr_J(4, 3, 7, 0, hAB, HRR_INT__h_d_k_s, HRR_INT__g_d_k_s, HRR_INT__g_f_k_s);
1313 
1314                 // form INT__g_f_l_s
1315                 ostei_general_hrr_J(4, 3, 8, 0, hAB, HRR_INT__h_d_l_s, HRR_INT__g_d_l_s, HRR_INT__g_f_l_s);
1316 
1317                 // form INT__h_f_h_s
1318                 ostei_general_hrr_J(5, 3, 5, 0, hAB, HRR_INT__i_d_h_s, HRR_INT__h_d_h_s, HRR_INT__h_f_h_s);
1319 
1320                 // form INT__h_f_i_s
1321                 ostei_general_hrr_J(5, 3, 6, 0, hAB, HRR_INT__i_d_i_s, HRR_INT__h_d_i_s, HRR_INT__h_f_i_s);
1322 
1323                 // form INT__h_f_k_s
1324                 ostei_general_hrr_J(5, 3, 7, 0, hAB, HRR_INT__i_d_k_s, HRR_INT__h_d_k_s, HRR_INT__h_f_k_s);
1325 
1326                 // form INT__h_f_l_s
1327                 ostei_general_hrr_J(5, 3, 8, 0, hAB, HRR_INT__i_d_l_s, HRR_INT__h_d_l_s, HRR_INT__h_f_l_s);
1328 
1329                 // form INT__g_g_h_s
1330                 ostei_general_hrr_J(4, 4, 5, 0, hAB, HRR_INT__h_f_h_s, HRR_INT__g_f_h_s, HRR_INT__g_g_h_s);
1331 
1332                 // form INT__g_g_i_s
1333                 ostei_general_hrr_J(4, 4, 6, 0, hAB, HRR_INT__h_f_i_s, HRR_INT__g_f_i_s, HRR_INT__g_g_i_s);
1334 
1335                 // form INT__g_g_k_s
1336                 ostei_general_hrr_J(4, 4, 7, 0, hAB, HRR_INT__h_f_k_s, HRR_INT__g_f_k_s, HRR_INT__g_g_k_s);
1337 
1338                 // form INT__g_g_l_s
1339                 ostei_general_hrr_J(4, 4, 8, 0, hAB, HRR_INT__h_f_l_s, HRR_INT__g_f_l_s, HRR_INT__g_g_l_s);
1340 
1341                 // form INT__g_g_h_p
1342                 ostei_general_hrr_L(4, 4, 5, 1, hCD, HRR_INT__g_g_i_s, HRR_INT__g_g_h_s, HRR_INT__g_g_h_p);
1343 
1344                 // form INT__g_g_i_p
1345                 ostei_general_hrr_L(4, 4, 6, 1, hCD, HRR_INT__g_g_k_s, HRR_INT__g_g_i_s, HRR_INT__g_g_i_p);
1346 
1347                 // form INT__g_g_k_p
1348                 ostei_general_hrr_L(4, 4, 7, 1, hCD, HRR_INT__g_g_l_s, HRR_INT__g_g_k_s, HRR_INT__g_g_k_p);
1349 
1350                 // form INT__g_g_h_d
1351                 ostei_general_hrr_L(4, 4, 5, 2, hCD, HRR_INT__g_g_i_p, HRR_INT__g_g_h_p, HRR_INT__g_g_h_d);
1352 
1353                 // form INT__g_g_i_d
1354                 ostei_general_hrr_L(4, 4, 6, 2, hCD, HRR_INT__g_g_k_p, HRR_INT__g_g_i_p, HRR_INT__g_g_i_d);
1355 
1356                 // form INT__g_g_h_f
1357                 ostei_general_hrr_L(4, 4, 5, 3, hCD, HRR_INT__g_g_i_d, HRR_INT__g_g_h_d, HRR_INT__g_g_h_f);
1358 
1359 
1360             }  // close HRR loop
1361 
1362 
1363         }   // close loop cdbatch
1364 
1365         istart = iend;
1366     }  // close loop over ab
1367 
1368     return P.nshell12_clip * Q.nshell12_clip;
1369 }
1370 
ostei_g_g_f_h(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__g_g_f_h)1371 int ostei_g_g_f_h(struct simint_multi_shellpair const P,
1372                   struct simint_multi_shellpair const Q,
1373                   double screen_tol,
1374                   double * const restrict work,
1375                   double * const restrict INT__g_g_f_h)
1376 {
1377     double Q_AB[3*Q.nshell12];
1378     struct simint_multi_shellpair Q_tmp = Q;
1379     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
1380     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
1381     Q_tmp.AB_x = Q_AB;
1382     Q_tmp.AB_y = Q_AB + Q.nshell12;
1383     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
1384 
1385     for(int i = 0; i < Q.nshell12; i++)
1386     {
1387         Q_tmp.AB_x[i] = -Q.AB_x[i];
1388         Q_tmp.AB_y[i] = -Q.AB_y[i];
1389         Q_tmp.AB_z[i] = -Q.AB_z[i];
1390     }
1391 
1392     int ret = ostei_g_g_h_f(P, Q_tmp, screen_tol, work, INT__g_g_f_h);
1393     double buffer[47250] SIMINT_ALIGN_ARRAY_DBL;
1394 
1395     for(int q = 0; q < ret; q++)
1396     {
1397         int idx = 0;
1398         for(int a = 0; a < 15; ++a)
1399         for(int b = 0; b < 15; ++b)
1400         for(int c = 0; c < 10; ++c)
1401         for(int d = 0; d < 21; ++d)
1402             buffer[idx++] = INT__g_g_f_h[q*47250+a*3150+b*210+d*10+c];
1403 
1404         memcpy(INT__g_g_f_h+q*47250, buffer, 47250*sizeof(double));
1405     }
1406 
1407     return ret;
1408 }
1409 
1410