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_p_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__i_p_h_f)8 int ostei_i_p_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__i_p_h_f)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__i_p_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__i_s_h_s = work + (SIMINT_NSHELL_SIMD * 0);
30     double * const INT__i_s_i_s = work + (SIMINT_NSHELL_SIMD * 588);
31     double * const INT__i_s_k_s = work + (SIMINT_NSHELL_SIMD * 1372);
32     double * const INT__i_s_l_s = work + (SIMINT_NSHELL_SIMD * 2380);
33     double * const INT__k_s_h_s = work + (SIMINT_NSHELL_SIMD * 3640);
34     double * const INT__k_s_i_s = work + (SIMINT_NSHELL_SIMD * 4396);
35     double * const INT__k_s_k_s = work + (SIMINT_NSHELL_SIMD * 5404);
36     double * const INT__k_s_l_s = work + (SIMINT_NSHELL_SIMD * 6700);
37     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*8320);
38     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
39     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_p_s = primwork + 16;
40     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_d_s = primwork + 61;
41     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_f_s = primwork + 145;
42     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_g_s = primwork + 275;
43     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_h_s = primwork + 455;
44     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_i_s = primwork + 686;
45     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_k_s = primwork + 966;
46     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_l_s = primwork + 1290;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_s = primwork + 1650;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_p_s = primwork + 1671;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_d_s = primwork + 1734;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_f_s = primwork + 1860;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_g_s = primwork + 2070;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_h_s = primwork + 2385;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_i_s = primwork + 2826;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_k_s = primwork + 3414;
55     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_l_s = primwork + 4170;
56     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 5115;
57     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_p_s = primwork + 5151;
58     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_d_s = primwork + 5259;
59     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_f_s = primwork + 5475;
60     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_g_s = primwork + 5835;
61     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_h_s = primwork + 6375;
62     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_i_s = primwork + 7131;
63     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_k_s = primwork + 8139;
64     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_l_s = primwork + 9435;
65     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_p_s = primwork + 11055;
66     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_d_s = primwork + 11205;
67     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_f_s = primwork + 11505;
68     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_g_s = primwork + 12005;
69     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_h_s = primwork + 12755;
70     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_i_s = primwork + 13805;
71     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_k_s = primwork + 15205;
72     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_l_s = primwork + 17005;
73     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_d_s = primwork + 19255;
74     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_f_s = primwork + 19615;
75     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_g_s = primwork + 20215;
76     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_h_s = primwork + 21115;
77     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_i_s = primwork + 22375;
78     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_k_s = primwork + 24055;
79     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_l_s = primwork + 26215;
80     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_f_s = primwork + 28915;
81     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_g_s = primwork + 29545;
82     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_h_s = primwork + 30490;
83     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_i_s = primwork + 31813;
84     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_k_s = primwork + 33577;
85     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_l_s = primwork + 35845;
86     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_g_s = primwork + 38680;
87     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_h_s = primwork + 39520;
88     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_i_s = primwork + 40696;
89     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_k_s = primwork + 42264;
90     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_l_s = primwork + 44280;
91     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_h_s = primwork + 46800;
92     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_i_s = primwork + 47556;
93     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_k_s = primwork + 48564;
94     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_l_s = primwork + 49860;
95     double * const hrrwork = (double *)(primwork + 51480);
96     double * const HRR_INT__i_p_h_s = hrrwork + 0;
97     double * const HRR_INT__i_p_h_p = hrrwork + 1764;
98     double * const HRR_INT__i_p_h_d = hrrwork + 7056;
99     double * const HRR_INT__i_p_i_s = hrrwork + 17640;
100     double * const HRR_INT__i_p_i_p = hrrwork + 19992;
101     double * const HRR_INT__i_p_i_d = hrrwork + 27048;
102     double * const HRR_INT__i_p_k_s = hrrwork + 41160;
103     double * const HRR_INT__i_p_k_p = hrrwork + 44184;
104     double * const HRR_INT__i_p_l_s = hrrwork + 53256;
105 
106 
107     // Create constants
108     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
109     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
110     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
111     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
112     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
113     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
114     const SIMINT_DBLTYPE const_7 = SIMINT_DBLSET1(7);
115     const SIMINT_DBLTYPE const_8 = SIMINT_DBLSET1(8);
116     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
117 
118 
119     ////////////////////////////////////////
120     // Loop over shells and primitives
121     ////////////////////////////////////////
122 
123     real_abcd = 0;
124     istart = 0;
125     for(ab = 0; ab < P.nshell12_clip; ++ab)
126     {
127         const int iend = istart + P.nprim12[ab];
128 
129         cd = 0;
130         jstart = 0;
131 
132         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
133         {
134             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
135             int jend = jstart;
136             for(i = 0; i < nshellbatch; i++)
137                 jend += Q.nprim12[cd+i];
138 
139             // Clear the beginning of the workspace (where we are accumulating integrals)
140             memset(work, 0, SIMINT_NSHELL_SIMD * 8320 * sizeof(double));
141             abcd = 0;
142 
143 
144             for(i = istart; i < iend; ++i)
145             {
146                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
147 
148                 if(check_screen)
149                 {
150                     // Skip this whole thing if always insignificant
151                     if((P.screen[i] * Q.screen_max) < screen_tol)
152                         continue;
153                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
154                 }
155 
156                 icd = 0;
157                 iprimcd = 0;
158                 nprim_icd = Q.nprim12[cd];
159                 double * restrict PRIM_PTR_INT__i_s_h_s = INT__i_s_h_s + abcd * 588;
160                 double * restrict PRIM_PTR_INT__i_s_i_s = INT__i_s_i_s + abcd * 784;
161                 double * restrict PRIM_PTR_INT__i_s_k_s = INT__i_s_k_s + abcd * 1008;
162                 double * restrict PRIM_PTR_INT__i_s_l_s = INT__i_s_l_s + abcd * 1260;
163                 double * restrict PRIM_PTR_INT__k_s_h_s = INT__k_s_h_s + abcd * 756;
164                 double * restrict PRIM_PTR_INT__k_s_i_s = INT__k_s_i_s + abcd * 1008;
165                 double * restrict PRIM_PTR_INT__k_s_k_s = INT__k_s_k_s + abcd * 1296;
166                 double * restrict PRIM_PTR_INT__k_s_l_s = INT__k_s_l_s + abcd * 1620;
167 
168 
169 
170                 // Load these one per loop over i
171                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
172                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
173                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
174 
175                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
176 
177                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
178                 {
179                     // calculate the shell offsets
180                     // these are the offset from the shell pointed to by cd
181                     // for each element
182                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
183                     int lastoffset = 0;
184                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
185 
186                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
187                     {
188                         // Handle if the first element of the vector is a new shell
189                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
190                         {
191                             nprim_icd += Q.nprim12[cd + (++icd)];
192                             PRIM_PTR_INT__i_s_h_s += 588;
193                             PRIM_PTR_INT__i_s_i_s += 784;
194                             PRIM_PTR_INT__i_s_k_s += 1008;
195                             PRIM_PTR_INT__i_s_l_s += 1260;
196                             PRIM_PTR_INT__k_s_h_s += 756;
197                             PRIM_PTR_INT__k_s_i_s += 1008;
198                             PRIM_PTR_INT__k_s_k_s += 1296;
199                             PRIM_PTR_INT__k_s_l_s += 1620;
200                         }
201                         iprimcd++;
202                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
203                         {
204                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
205                             {
206                                 shelloffsets[n] = shelloffsets[n-1] + 1;
207                                 lastoffset++;
208                                 nprim_icd += Q.nprim12[cd + (++icd)];
209                             }
210                             else
211                                 shelloffsets[n] = shelloffsets[n-1];
212                             iprimcd++;
213                         }
214                     }
215                     else
216                         iprimcd += SIMINT_SIMD_LEN;
217 
218                     // Do we have to compute this vector (or has it been screened out)?
219                     // (not_screened != 0 means we have to do this vector)
220                     if(check_screen)
221                     {
222                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
223                         if(vmax < screen_tol)
224                         {
225                             PRIM_PTR_INT__i_s_h_s += lastoffset*588;
226                             PRIM_PTR_INT__i_s_i_s += lastoffset*784;
227                             PRIM_PTR_INT__i_s_k_s += lastoffset*1008;
228                             PRIM_PTR_INT__i_s_l_s += lastoffset*1260;
229                             PRIM_PTR_INT__k_s_h_s += lastoffset*756;
230                             PRIM_PTR_INT__k_s_i_s += lastoffset*1008;
231                             PRIM_PTR_INT__k_s_k_s += lastoffset*1296;
232                             PRIM_PTR_INT__k_s_l_s += lastoffset*1620;
233                             continue;
234                         }
235                     }
236 
237                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
238                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
239                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
240                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
241 
242 
243                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
244                     SIMINT_DBLTYPE PQ[3];
245                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
246                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
247                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
248                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
249                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
250                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
251 
252                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
253                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
254                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
255                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
256                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
257                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
258                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
259 
260                     // NOTE: Minus sign!
261                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
262                     SIMINT_DBLTYPE aop_PQ[3];
263                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
264                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
265                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
266 
267                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
268                     SIMINT_DBLTYPE aoq_PQ[3];
269                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
270                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
271                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
272                     // Put a minus sign here so we don't have to in RR routines
273                     a_over_q = SIMINT_NEG(a_over_q);
274 
275 
276                     //////////////////////////////////////////////
277                     // Fjt function section
278                     // Maximum v value: 15
279                     //////////////////////////////////////////////
280                     // The parameter to the Fjt function
281                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
282 
283 
284                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
285 
286 
287                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 15);
288                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
289                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
290                     for(n = 0; n <= 15; n++)
291                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
292 
293                     //////////////////////////////////////////////
294                     // Primitive integrals: Vertical recurrance
295                     //////////////////////////////////////////////
296 
297                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
298                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
299                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
300                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
301                     const SIMINT_DBLTYPE vrr_const_5_over_2p = SIMINT_MUL(const_5, one_over_2p);
302                     const SIMINT_DBLTYPE vrr_const_6_over_2p = SIMINT_MUL(const_6, one_over_2p);
303                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
304                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
305                     const SIMINT_DBLTYPE vrr_const_3_over_2q = SIMINT_MUL(const_3, one_over_2q);
306                     const SIMINT_DBLTYPE vrr_const_4_over_2q = SIMINT_MUL(const_4, one_over_2q);
307                     const SIMINT_DBLTYPE vrr_const_5_over_2q = SIMINT_MUL(const_5, one_over_2q);
308                     const SIMINT_DBLTYPE vrr_const_6_over_2q = SIMINT_MUL(const_6, one_over_2q);
309                     const SIMINT_DBLTYPE vrr_const_7_over_2q = SIMINT_MUL(const_7, one_over_2q);
310                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
311                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
312                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
313                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
314                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
315                     const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
316                     const SIMINT_DBLTYPE vrr_const_7_over_2pq = SIMINT_MUL(const_7, one_over_2pq);
317                     const SIMINT_DBLTYPE vrr_const_8_over_2pq = SIMINT_MUL(const_8, one_over_2pq);
318 
319 
320 
321                     // Forming PRIM_INT__s_s_p_s[15 * 3];
322                     for(n = 0; n < 15; ++n)  // loop over orders of auxiliary function
323                     {
324 
325                         PRIM_INT__s_s_p_s[n * 3 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
326                         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]);
327 
328                         PRIM_INT__s_s_p_s[n * 3 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
329                         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]);
330 
331                         PRIM_INT__s_s_p_s[n * 3 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
332                         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]);
333 
334                     }
335 
336 
337 
338                     // Forming PRIM_INT__s_s_d_s[14 * 6];
339                     for(n = 0; n < 14; ++n)  // loop over orders of auxiliary function
340                     {
341 
342                         PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_p_s[n * 3 + 0]);
343                         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]);
344                         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]);
345 
346                         PRIM_INT__s_s_d_s[n * 6 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_p_s[n * 3 + 0]);
347                         PRIM_INT__s_s_d_s[n * 6 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_d_s[n * 6 + 1]);
348 
349                         PRIM_INT__s_s_d_s[n * 6 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 0]);
350                         PRIM_INT__s_s_d_s[n * 6 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_d_s[n * 6 + 2]);
351 
352                         PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_p_s[n * 3 + 1]);
353                         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]);
354                         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]);
355 
356                         PRIM_INT__s_s_d_s[n * 6 + 4] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 1]);
357                         PRIM_INT__s_s_d_s[n * 6 + 4] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_d_s[n * 6 + 4]);
358 
359                         PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 2]);
360                         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]);
361                         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]);
362 
363                     }
364 
365 
366 
367                     // Forming PRIM_INT__s_s_f_s[13 * 10];
368                     for(n = 0; n < 13; ++n)  // loop over orders of auxiliary function
369                     {
370 
371                         PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 0]);
372                         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]);
373                         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]);
374 
375                         PRIM_INT__s_s_f_s[n * 10 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 0]);
376                         PRIM_INT__s_s_f_s[n * 10 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__s_s_f_s[n * 10 + 1]);
377 
378                         PRIM_INT__s_s_f_s[n * 10 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 0]);
379                         PRIM_INT__s_s_f_s[n * 10 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__s_s_f_s[n * 10 + 2]);
380 
381                         PRIM_INT__s_s_f_s[n * 10 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 3]);
382                         PRIM_INT__s_s_f_s[n * 10 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__s_s_f_s[n * 10 + 3]);
383 
384                         PRIM_INT__s_s_f_s[n * 10 + 4] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 1]);
385                         PRIM_INT__s_s_f_s[n * 10 + 4] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__s_s_f_s[n * 10 + 4]);
386 
387                         PRIM_INT__s_s_f_s[n * 10 + 5] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 5]);
388                         PRIM_INT__s_s_f_s[n * 10 + 5] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__s_s_f_s[n * 10 + 5]);
389 
390                         PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 3]);
391                         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]);
392                         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]);
393 
394                         PRIM_INT__s_s_f_s[n * 10 + 7] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 3]);
395                         PRIM_INT__s_s_f_s[n * 10 + 7] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__s_s_f_s[n * 10 + 7]);
396 
397                         PRIM_INT__s_s_f_s[n * 10 + 8] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 5]);
398                         PRIM_INT__s_s_f_s[n * 10 + 8] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__s_s_f_s[n * 10 + 8]);
399 
400                         PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 5]);
401                         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]);
402                         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]);
403 
404                     }
405 
406 
407                     VRR_K_s_s_g_s(
408                             PRIM_INT__s_s_g_s,
409                             PRIM_INT__s_s_f_s,
410                             PRIM_INT__s_s_d_s,
411                             Q_PA,
412                             a_over_q,
413                             aoq_PQ,
414                             one_over_2q,
415                             12);
416 
417 
418                     VRR_K_s_s_h_s(
419                             PRIM_INT__s_s_h_s,
420                             PRIM_INT__s_s_g_s,
421                             PRIM_INT__s_s_f_s,
422                             Q_PA,
423                             a_over_q,
424                             aoq_PQ,
425                             one_over_2q,
426                             11);
427 
428 
429                     ostei_general_vrr_I(1, 0, 5, 0, 7,
430                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
431                             PRIM_INT__s_s_h_s, NULL, NULL, PRIM_INT__s_s_g_s, NULL, PRIM_INT__p_s_h_s);
432 
433 
434                     VRR_I_p_s_g_s(
435                             PRIM_INT__p_s_g_s,
436                             PRIM_INT__s_s_g_s,
437                             PRIM_INT__s_s_f_s,
438                             P_PA,
439                             aop_PQ,
440                             one_over_2pq,
441                             7);
442 
443 
444                     ostei_general_vrr_I(2, 0, 5, 0, 6,
445                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
446                             PRIM_INT__p_s_h_s, PRIM_INT__s_s_h_s, NULL, PRIM_INT__p_s_g_s, NULL, PRIM_INT__d_s_h_s);
447 
448 
449                     VRR_I_p_s_f_s(
450                             PRIM_INT__p_s_f_s,
451                             PRIM_INT__s_s_f_s,
452                             PRIM_INT__s_s_d_s,
453                             P_PA,
454                             aop_PQ,
455                             one_over_2pq,
456                             7);
457 
458 
459                     ostei_general_vrr_I(2, 0, 4, 0, 6,
460                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
461                             PRIM_INT__p_s_g_s, PRIM_INT__s_s_g_s, NULL, PRIM_INT__p_s_f_s, NULL, PRIM_INT__d_s_g_s);
462 
463 
464                     ostei_general_vrr_I(3, 0, 5, 0, 5,
465                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
466                             PRIM_INT__d_s_h_s, PRIM_INT__p_s_h_s, NULL, PRIM_INT__d_s_g_s, NULL, PRIM_INT__f_s_h_s);
467 
468 
469 
470                     // Forming PRIM_INT__p_s_d_s[7 * 18];
471                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
472                     {
473 
474                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 0]);
475                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__p_s_d_s[n * 18 + 0]);
476                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_d_s[n * 18 + 0]);
477 
478                         PRIM_INT__p_s_d_s[n * 18 + 1] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 1]);
479                         PRIM_INT__p_s_d_s[n * 18 + 1] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__p_s_d_s[n * 18 + 1]);
480                         PRIM_INT__p_s_d_s[n * 18 + 1] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_d_s[n * 18 + 1]);
481 
482                         PRIM_INT__p_s_d_s[n * 18 + 2] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 2]);
483                         PRIM_INT__p_s_d_s[n * 18 + 2] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 2], PRIM_INT__p_s_d_s[n * 18 + 2]);
484                         PRIM_INT__p_s_d_s[n * 18 + 2] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_d_s[n * 18 + 2]);
485 
486                         PRIM_INT__p_s_d_s[n * 18 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 3]);
487                         PRIM_INT__p_s_d_s[n * 18 + 3] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__p_s_d_s[n * 18 + 3]);
488 
489                         PRIM_INT__p_s_d_s[n * 18 + 4] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 4]);
490                         PRIM_INT__p_s_d_s[n * 18 + 4] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 4], PRIM_INT__p_s_d_s[n * 18 + 4]);
491 
492                         PRIM_INT__p_s_d_s[n * 18 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 5]);
493                         PRIM_INT__p_s_d_s[n * 18 + 5] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__p_s_d_s[n * 18 + 5]);
494 
495                         PRIM_INT__p_s_d_s[n * 18 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 0]);
496                         PRIM_INT__p_s_d_s[n * 18 + 6] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__p_s_d_s[n * 18 + 6]);
497 
498                         PRIM_INT__p_s_d_s[n * 18 + 7] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 1]);
499                         PRIM_INT__p_s_d_s[n * 18 + 7] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__p_s_d_s[n * 18 + 7]);
500                         PRIM_INT__p_s_d_s[n * 18 + 7] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_d_s[n * 18 + 7]);
501 
502                         PRIM_INT__p_s_d_s[n * 18 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 2]);
503                         PRIM_INT__p_s_d_s[n * 18 + 8] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 2], PRIM_INT__p_s_d_s[n * 18 + 8]);
504 
505                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 3]);
506                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__p_s_d_s[n * 18 + 9]);
507                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_d_s[n * 18 + 9]);
508 
509                         PRIM_INT__p_s_d_s[n * 18 + 10] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 4]);
510                         PRIM_INT__p_s_d_s[n * 18 + 10] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 4], PRIM_INT__p_s_d_s[n * 18 + 10]);
511                         PRIM_INT__p_s_d_s[n * 18 + 10] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_d_s[n * 18 + 10]);
512 
513                         PRIM_INT__p_s_d_s[n * 18 + 11] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 5]);
514                         PRIM_INT__p_s_d_s[n * 18 + 11] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__p_s_d_s[n * 18 + 11]);
515 
516                         PRIM_INT__p_s_d_s[n * 18 + 12] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 0]);
517                         PRIM_INT__p_s_d_s[n * 18 + 12] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__p_s_d_s[n * 18 + 12]);
518 
519                         PRIM_INT__p_s_d_s[n * 18 + 13] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 1]);
520                         PRIM_INT__p_s_d_s[n * 18 + 13] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__p_s_d_s[n * 18 + 13]);
521 
522                         PRIM_INT__p_s_d_s[n * 18 + 14] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 2]);
523                         PRIM_INT__p_s_d_s[n * 18 + 14] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 2], PRIM_INT__p_s_d_s[n * 18 + 14]);
524                         PRIM_INT__p_s_d_s[n * 18 + 14] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_d_s[n * 18 + 14]);
525 
526                         PRIM_INT__p_s_d_s[n * 18 + 15] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 3]);
527                         PRIM_INT__p_s_d_s[n * 18 + 15] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__p_s_d_s[n * 18 + 15]);
528 
529                         PRIM_INT__p_s_d_s[n * 18 + 16] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 4]);
530                         PRIM_INT__p_s_d_s[n * 18 + 16] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 4], PRIM_INT__p_s_d_s[n * 18 + 16]);
531                         PRIM_INT__p_s_d_s[n * 18 + 16] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_d_s[n * 18 + 16]);
532 
533                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 5]);
534                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__p_s_d_s[n * 18 + 17]);
535                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_d_s[n * 18 + 17]);
536 
537                     }
538 
539 
540                     VRR_I_d_s_f_s(
541                             PRIM_INT__d_s_f_s,
542                             PRIM_INT__p_s_f_s,
543                             PRIM_INT__s_s_f_s,
544                             PRIM_INT__p_s_d_s,
545                             P_PA,
546                             a_over_p,
547                             aop_PQ,
548                             one_over_2p,
549                             one_over_2pq,
550                             6);
551 
552 
553                     ostei_general_vrr_I(3, 0, 4, 0, 5,
554                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
555                             PRIM_INT__d_s_g_s, PRIM_INT__p_s_g_s, NULL, PRIM_INT__d_s_f_s, NULL, PRIM_INT__f_s_g_s);
556 
557 
558                     ostei_general_vrr_I(4, 0, 5, 0, 4,
559                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
560                             PRIM_INT__f_s_h_s, PRIM_INT__d_s_h_s, NULL, PRIM_INT__f_s_g_s, NULL, PRIM_INT__g_s_h_s);
561 
562 
563 
564                     // Forming PRIM_INT__p_s_p_s[7 * 9];
565                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
566                     {
567 
568                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_p_s[n * 3 + 0]);
569                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 0]);
570                         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]);
571 
572                         PRIM_INT__p_s_p_s[n * 9 + 1] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_p_s[n * 3 + 1]);
573                         PRIM_INT__p_s_p_s[n * 9 + 1] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 1]);
574 
575                         PRIM_INT__p_s_p_s[n * 9 + 2] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_p_s[n * 3 + 2]);
576                         PRIM_INT__p_s_p_s[n * 9 + 2] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 2]);
577 
578                         PRIM_INT__p_s_p_s[n * 9 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_p_s[n * 3 + 0]);
579                         PRIM_INT__p_s_p_s[n * 9 + 3] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 3]);
580 
581                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_p_s[n * 3 + 1]);
582                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 4]);
583                         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]);
584 
585                         PRIM_INT__p_s_p_s[n * 9 + 5] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_p_s[n * 3 + 2]);
586                         PRIM_INT__p_s_p_s[n * 9 + 5] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 5]);
587 
588                         PRIM_INT__p_s_p_s[n * 9 + 6] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_p_s[n * 3 + 0]);
589                         PRIM_INT__p_s_p_s[n * 9 + 6] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 6]);
590 
591                         PRIM_INT__p_s_p_s[n * 9 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_p_s[n * 3 + 1]);
592                         PRIM_INT__p_s_p_s[n * 9 + 7] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 7]);
593 
594                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_p_s[n * 3 + 2]);
595                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 8]);
596                         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]);
597 
598                     }
599 
600 
601                     VRR_I_d_s_d_s(
602                             PRIM_INT__d_s_d_s,
603                             PRIM_INT__p_s_d_s,
604                             PRIM_INT__s_s_d_s,
605                             PRIM_INT__p_s_p_s,
606                             P_PA,
607                             a_over_p,
608                             aop_PQ,
609                             one_over_2p,
610                             one_over_2pq,
611                             6);
612 
613 
614                     ostei_general_vrr_I(3, 0, 3, 0, 5,
615                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
616                             PRIM_INT__d_s_f_s, PRIM_INT__p_s_f_s, NULL, PRIM_INT__d_s_d_s, NULL, PRIM_INT__f_s_f_s);
617 
618 
619                     ostei_general_vrr_I(4, 0, 4, 0, 4,
620                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
621                             PRIM_INT__f_s_g_s, PRIM_INT__d_s_g_s, NULL, PRIM_INT__f_s_f_s, NULL, PRIM_INT__g_s_g_s);
622 
623 
624                     ostei_general_vrr_I(5, 0, 5, 0, 3,
625                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
626                             PRIM_INT__g_s_h_s, PRIM_INT__f_s_h_s, NULL, PRIM_INT__g_s_g_s, NULL, PRIM_INT__h_s_h_s);
627 
628 
629 
630                     // Forming PRIM_INT__p_s_s_s[7 * 3];
631                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
632                     {
633 
634                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
635                         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]);
636 
637                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
638                         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]);
639 
640                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
641                         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]);
642 
643                     }
644 
645 
646 
647                     // Forming PRIM_INT__d_s_p_s[6 * 18];
648                     for(n = 0; n < 6; ++n)  // loop over orders of auxiliary function
649                     {
650 
651                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_p_s[n * 9 + 0]);
652                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__p_s_p_s[(n+1) * 9 + 0], PRIM_INT__d_s_p_s[n * 18 + 0]);
653                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_p_s[n * 3 + 0]), PRIM_INT__d_s_p_s[n * 18 + 0]);
654                         PRIM_INT__d_s_p_s[n * 18 + 0] = 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 + 0]);
655 
656                         PRIM_INT__d_s_p_s[n * 18 + 1] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_p_s[n * 9 + 1]);
657                         PRIM_INT__d_s_p_s[n * 18 + 1] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__p_s_p_s[(n+1) * 9 + 1], PRIM_INT__d_s_p_s[n * 18 + 1]);
658                         PRIM_INT__d_s_p_s[n * 18 + 1] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_p_s[n * 3 + 1]), PRIM_INT__d_s_p_s[n * 18 + 1]);
659 
660                         PRIM_INT__d_s_p_s[n * 18 + 2] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_p_s[n * 9 + 2]);
661                         PRIM_INT__d_s_p_s[n * 18 + 2] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__p_s_p_s[(n+1) * 9 + 2], PRIM_INT__d_s_p_s[n * 18 + 2]);
662                         PRIM_INT__d_s_p_s[n * 18 + 2] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__s_s_p_s[n * 3 + 2]), PRIM_INT__d_s_p_s[n * 18 + 2]);
663 
664                         PRIM_INT__d_s_p_s[n * 18 + 9] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_p_s[n * 9 + 3]);
665                         PRIM_INT__d_s_p_s[n * 18 + 9] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__p_s_p_s[(n+1) * 9 + 3], PRIM_INT__d_s_p_s[n * 18 + 9]);
666                         PRIM_INT__d_s_p_s[n * 18 + 9] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_p_s[n * 3 + 0]), PRIM_INT__d_s_p_s[n * 18 + 9]);
667 
668                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_p_s[n * 9 + 4]);
669                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__p_s_p_s[(n+1) * 9 + 4], PRIM_INT__d_s_p_s[n * 18 + 10]);
670                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_p_s[n * 3 + 1]), PRIM_INT__d_s_p_s[n * 18 + 10]);
671                         PRIM_INT__d_s_p_s[n * 18 + 10] = 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 + 10]);
672 
673                         PRIM_INT__d_s_p_s[n * 18 + 11] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_p_s[n * 9 + 5]);
674                         PRIM_INT__d_s_p_s[n * 18 + 11] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__p_s_p_s[(n+1) * 9 + 5], PRIM_INT__d_s_p_s[n * 18 + 11]);
675                         PRIM_INT__d_s_p_s[n * 18 + 11] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__s_s_p_s[n * 3 + 2]), PRIM_INT__d_s_p_s[n * 18 + 11]);
676 
677                         PRIM_INT__d_s_p_s[n * 18 + 15] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_p_s[n * 9 + 6]);
678                         PRIM_INT__d_s_p_s[n * 18 + 15] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__p_s_p_s[(n+1) * 9 + 6], PRIM_INT__d_s_p_s[n * 18 + 15]);
679                         PRIM_INT__d_s_p_s[n * 18 + 15] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_p_s[n * 3 + 0]), PRIM_INT__d_s_p_s[n * 18 + 15]);
680 
681                         PRIM_INT__d_s_p_s[n * 18 + 16] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_p_s[n * 9 + 7]);
682                         PRIM_INT__d_s_p_s[n * 18 + 16] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__p_s_p_s[(n+1) * 9 + 7], PRIM_INT__d_s_p_s[n * 18 + 16]);
683                         PRIM_INT__d_s_p_s[n * 18 + 16] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_p_s[n * 3 + 1]), PRIM_INT__d_s_p_s[n * 18 + 16]);
684 
685                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_p_s[n * 9 + 8]);
686                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__p_s_p_s[(n+1) * 9 + 8], PRIM_INT__d_s_p_s[n * 18 + 17]);
687                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__s_s_p_s[n * 3 + 2]), PRIM_INT__d_s_p_s[n * 18 + 17]);
688                         PRIM_INT__d_s_p_s[n * 18 + 17] = 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 + 17]);
689 
690                     }
691 
692 
693                     VRR_I_f_s_d_s(
694                             PRIM_INT__f_s_d_s,
695                             PRIM_INT__d_s_d_s,
696                             PRIM_INT__p_s_d_s,
697                             PRIM_INT__d_s_p_s,
698                             P_PA,
699                             a_over_p,
700                             aop_PQ,
701                             one_over_2p,
702                             one_over_2pq,
703                             5);
704 
705 
706                     ostei_general_vrr_I(4, 0, 3, 0, 4,
707                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
708                             PRIM_INT__f_s_f_s, PRIM_INT__d_s_f_s, NULL, PRIM_INT__f_s_d_s, NULL, PRIM_INT__g_s_f_s);
709 
710 
711                     ostei_general_vrr_I(5, 0, 4, 0, 3,
712                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
713                             PRIM_INT__g_s_g_s, PRIM_INT__f_s_g_s, NULL, PRIM_INT__g_s_f_s, NULL, PRIM_INT__h_s_g_s);
714 
715 
716                     ostei_general_vrr_I(6, 0, 5, 0, 2,
717                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
718                             PRIM_INT__h_s_h_s, PRIM_INT__g_s_h_s, NULL, PRIM_INT__h_s_g_s, NULL, PRIM_INT__i_s_h_s);
719 
720 
721                     ostei_general_vrr1_K(6, 10,
722                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
723                             PRIM_INT__s_s_h_s, PRIM_INT__s_s_g_s, PRIM_INT__s_s_i_s);
724 
725 
726                     ostei_general_vrr_I(1, 0, 6, 0, 7,
727                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
728                             PRIM_INT__s_s_i_s, NULL, NULL, PRIM_INT__s_s_h_s, NULL, PRIM_INT__p_s_i_s);
729 
730 
731                     ostei_general_vrr_I(2, 0, 6, 0, 6,
732                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
733                             PRIM_INT__p_s_i_s, PRIM_INT__s_s_i_s, NULL, PRIM_INT__p_s_h_s, NULL, PRIM_INT__d_s_i_s);
734 
735 
736                     ostei_general_vrr_I(3, 0, 6, 0, 5,
737                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
738                             PRIM_INT__d_s_i_s, PRIM_INT__p_s_i_s, NULL, PRIM_INT__d_s_h_s, NULL, PRIM_INT__f_s_i_s);
739 
740 
741                     ostei_general_vrr_I(4, 0, 6, 0, 4,
742                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
743                             PRIM_INT__f_s_i_s, PRIM_INT__d_s_i_s, NULL, PRIM_INT__f_s_h_s, NULL, PRIM_INT__g_s_i_s);
744 
745 
746                     ostei_general_vrr_I(5, 0, 6, 0, 3,
747                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
748                             PRIM_INT__g_s_i_s, PRIM_INT__f_s_i_s, NULL, PRIM_INT__g_s_h_s, NULL, PRIM_INT__h_s_i_s);
749 
750 
751                     ostei_general_vrr_I(6, 0, 6, 0, 2,
752                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
753                             PRIM_INT__h_s_i_s, PRIM_INT__g_s_i_s, NULL, PRIM_INT__h_s_h_s, NULL, PRIM_INT__i_s_i_s);
754 
755 
756 
757                     // Forming PRIM_INT__d_s_s_s[6 * 6];
758                     for(n = 0; n < 6; ++n)  // loop over orders of auxiliary function
759                     {
760 
761                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
762                         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]);
763                         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]);
764 
765                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
766                         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]);
767                         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]);
768 
769                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
770                         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]);
771                         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]);
772 
773                     }
774 
775 
776                     VRR_I_f_s_p_s(
777                             PRIM_INT__f_s_p_s,
778                             PRIM_INT__d_s_p_s,
779                             PRIM_INT__p_s_p_s,
780                             PRIM_INT__d_s_s_s,
781                             P_PA,
782                             a_over_p,
783                             aop_PQ,
784                             one_over_2p,
785                             one_over_2pq,
786                             5);
787 
788 
789                     ostei_general_vrr_I(4, 0, 2, 0, 4,
790                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
791                             PRIM_INT__f_s_d_s, PRIM_INT__d_s_d_s, NULL, PRIM_INT__f_s_p_s, NULL, PRIM_INT__g_s_d_s);
792 
793 
794                     ostei_general_vrr_I(5, 0, 3, 0, 3,
795                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
796                             PRIM_INT__g_s_f_s, PRIM_INT__f_s_f_s, NULL, PRIM_INT__g_s_d_s, NULL, PRIM_INT__h_s_f_s);
797 
798 
799                     ostei_general_vrr_I(6, 0, 4, 0, 2,
800                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
801                             PRIM_INT__h_s_g_s, PRIM_INT__g_s_g_s, NULL, PRIM_INT__h_s_f_s, NULL, PRIM_INT__i_s_g_s);
802 
803 
804                     ostei_general_vrr_I(7, 0, 5, 0, 1,
805                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
806                             PRIM_INT__i_s_h_s, PRIM_INT__h_s_h_s, NULL, PRIM_INT__i_s_g_s, NULL, PRIM_INT__k_s_h_s);
807 
808 
809                     ostei_general_vrr1_K(7, 9,
810                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
811                             PRIM_INT__s_s_i_s, PRIM_INT__s_s_h_s, PRIM_INT__s_s_k_s);
812 
813 
814                     ostei_general_vrr_I(1, 0, 7, 0, 7,
815                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
816                             PRIM_INT__s_s_k_s, NULL, NULL, PRIM_INT__s_s_i_s, NULL, PRIM_INT__p_s_k_s);
817 
818 
819                     ostei_general_vrr_I(2, 0, 7, 0, 6,
820                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
821                             PRIM_INT__p_s_k_s, PRIM_INT__s_s_k_s, NULL, PRIM_INT__p_s_i_s, NULL, PRIM_INT__d_s_k_s);
822 
823 
824                     ostei_general_vrr_I(3, 0, 7, 0, 5,
825                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
826                             PRIM_INT__d_s_k_s, PRIM_INT__p_s_k_s, NULL, PRIM_INT__d_s_i_s, NULL, PRIM_INT__f_s_k_s);
827 
828 
829                     ostei_general_vrr_I(4, 0, 7, 0, 4,
830                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
831                             PRIM_INT__f_s_k_s, PRIM_INT__d_s_k_s, NULL, PRIM_INT__f_s_i_s, NULL, PRIM_INT__g_s_k_s);
832 
833 
834                     ostei_general_vrr_I(5, 0, 7, 0, 3,
835                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
836                             PRIM_INT__g_s_k_s, PRIM_INT__f_s_k_s, NULL, PRIM_INT__g_s_i_s, NULL, PRIM_INT__h_s_k_s);
837 
838 
839                     ostei_general_vrr_I(6, 0, 7, 0, 2,
840                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
841                             PRIM_INT__h_s_k_s, PRIM_INT__g_s_k_s, NULL, PRIM_INT__h_s_i_s, NULL, PRIM_INT__i_s_k_s);
842 
843 
844                     ostei_general_vrr_I(7, 0, 6, 0, 1,
845                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
846                             PRIM_INT__i_s_i_s, PRIM_INT__h_s_i_s, NULL, PRIM_INT__i_s_h_s, NULL, PRIM_INT__k_s_i_s);
847 
848 
849                     ostei_general_vrr1_K(8, 8,
850                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
851                             PRIM_INT__s_s_k_s, PRIM_INT__s_s_i_s, PRIM_INT__s_s_l_s);
852 
853 
854                     ostei_general_vrr_I(1, 0, 8, 0, 7,
855                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
856                             PRIM_INT__s_s_l_s, NULL, NULL, PRIM_INT__s_s_k_s, NULL, PRIM_INT__p_s_l_s);
857 
858 
859                     ostei_general_vrr_I(2, 0, 8, 0, 6,
860                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
861                             PRIM_INT__p_s_l_s, PRIM_INT__s_s_l_s, NULL, PRIM_INT__p_s_k_s, NULL, PRIM_INT__d_s_l_s);
862 
863 
864                     ostei_general_vrr_I(3, 0, 8, 0, 5,
865                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
866                             PRIM_INT__d_s_l_s, PRIM_INT__p_s_l_s, NULL, PRIM_INT__d_s_k_s, NULL, PRIM_INT__f_s_l_s);
867 
868 
869                     ostei_general_vrr_I(4, 0, 8, 0, 4,
870                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
871                             PRIM_INT__f_s_l_s, PRIM_INT__d_s_l_s, NULL, PRIM_INT__f_s_k_s, NULL, PRIM_INT__g_s_l_s);
872 
873 
874                     ostei_general_vrr_I(5, 0, 8, 0, 3,
875                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
876                             PRIM_INT__g_s_l_s, PRIM_INT__f_s_l_s, NULL, PRIM_INT__g_s_k_s, NULL, PRIM_INT__h_s_l_s);
877 
878 
879                     ostei_general_vrr_I(6, 0, 8, 0, 2,
880                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
881                             PRIM_INT__h_s_l_s, PRIM_INT__g_s_l_s, NULL, PRIM_INT__h_s_k_s, NULL, PRIM_INT__i_s_l_s);
882 
883 
884                     ostei_general_vrr_I(7, 0, 7, 0, 1,
885                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
886                             PRIM_INT__i_s_k_s, PRIM_INT__h_s_k_s, NULL, PRIM_INT__i_s_i_s, NULL, PRIM_INT__k_s_k_s);
887 
888 
889                     ostei_general_vrr_I(7, 0, 8, 0, 1,
890                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
891                             PRIM_INT__i_s_l_s, PRIM_INT__h_s_l_s, NULL, PRIM_INT__i_s_k_s, NULL, PRIM_INT__k_s_l_s);
892 
893 
894 
895 
896                     ////////////////////////////////////
897                     // Accumulate contracted integrals
898                     ////////////////////////////////////
899                     if(lastoffset == 0)
900                     {
901                         contract_all(588, PRIM_INT__i_s_h_s, PRIM_PTR_INT__i_s_h_s);
902                         contract_all(784, PRIM_INT__i_s_i_s, PRIM_PTR_INT__i_s_i_s);
903                         contract_all(1008, PRIM_INT__i_s_k_s, PRIM_PTR_INT__i_s_k_s);
904                         contract_all(1260, PRIM_INT__i_s_l_s, PRIM_PTR_INT__i_s_l_s);
905                         contract_all(756, PRIM_INT__k_s_h_s, PRIM_PTR_INT__k_s_h_s);
906                         contract_all(1008, PRIM_INT__k_s_i_s, PRIM_PTR_INT__k_s_i_s);
907                         contract_all(1296, PRIM_INT__k_s_k_s, PRIM_PTR_INT__k_s_k_s);
908                         contract_all(1620, PRIM_INT__k_s_l_s, PRIM_PTR_INT__k_s_l_s);
909                     }
910                     else
911                     {
912                         contract(588, shelloffsets, PRIM_INT__i_s_h_s, PRIM_PTR_INT__i_s_h_s);
913                         contract(784, shelloffsets, PRIM_INT__i_s_i_s, PRIM_PTR_INT__i_s_i_s);
914                         contract(1008, shelloffsets, PRIM_INT__i_s_k_s, PRIM_PTR_INT__i_s_k_s);
915                         contract(1260, shelloffsets, PRIM_INT__i_s_l_s, PRIM_PTR_INT__i_s_l_s);
916                         contract(756, shelloffsets, PRIM_INT__k_s_h_s, PRIM_PTR_INT__k_s_h_s);
917                         contract(1008, shelloffsets, PRIM_INT__k_s_i_s, PRIM_PTR_INT__k_s_i_s);
918                         contract(1296, shelloffsets, PRIM_INT__k_s_k_s, PRIM_PTR_INT__k_s_k_s);
919                         contract(1620, shelloffsets, PRIM_INT__k_s_l_s, PRIM_PTR_INT__k_s_l_s);
920                         PRIM_PTR_INT__i_s_h_s += lastoffset*588;
921                         PRIM_PTR_INT__i_s_i_s += lastoffset*784;
922                         PRIM_PTR_INT__i_s_k_s += lastoffset*1008;
923                         PRIM_PTR_INT__i_s_l_s += lastoffset*1260;
924                         PRIM_PTR_INT__k_s_h_s += lastoffset*756;
925                         PRIM_PTR_INT__k_s_i_s += lastoffset*1008;
926                         PRIM_PTR_INT__k_s_k_s += lastoffset*1296;
927                         PRIM_PTR_INT__k_s_l_s += lastoffset*1620;
928                     }
929 
930                 }  // close loop over j
931             }  // close loop over i
932 
933             //Advance to the next batch
934             jstart = SIMINT_SIMD_ROUND(jend);
935 
936             //////////////////////////////////////////////
937             // Contracted integrals: Horizontal recurrance
938             //////////////////////////////////////////////
939 
940 
941             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
942 
943 
944             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
945             {
946                 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
947 
948                 // set up HRR pointers
949                 double const * restrict HRR_INT__i_s_h_s = INT__i_s_h_s + abcd * 588;
950                 double const * restrict HRR_INT__i_s_i_s = INT__i_s_i_s + abcd * 784;
951                 double const * restrict HRR_INT__i_s_k_s = INT__i_s_k_s + abcd * 1008;
952                 double const * restrict HRR_INT__i_s_l_s = INT__i_s_l_s + abcd * 1260;
953                 double const * restrict HRR_INT__k_s_h_s = INT__k_s_h_s + abcd * 756;
954                 double const * restrict HRR_INT__k_s_i_s = INT__k_s_i_s + abcd * 1008;
955                 double const * restrict HRR_INT__k_s_k_s = INT__k_s_k_s + abcd * 1296;
956                 double const * restrict HRR_INT__k_s_l_s = INT__k_s_l_s + abcd * 1620;
957                 double * restrict HRR_INT__i_p_h_f = INT__i_p_h_f + real_abcd * 17640;
958 
959                 // form INT__i_p_h_s
960                 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);
961 
962                 // form INT__i_p_i_s
963                 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);
964 
965                 // form INT__i_p_k_s
966                 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);
967 
968                 // form INT__i_p_l_s
969                 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);
970 
971                 // form INT__i_p_h_p
972                 ostei_general_hrr_L(6, 1, 5, 1, hCD, HRR_INT__i_p_i_s, HRR_INT__i_p_h_s, HRR_INT__i_p_h_p);
973 
974                 // form INT__i_p_i_p
975                 ostei_general_hrr_L(6, 1, 6, 1, hCD, HRR_INT__i_p_k_s, HRR_INT__i_p_i_s, HRR_INT__i_p_i_p);
976 
977                 // form INT__i_p_k_p
978                 ostei_general_hrr_L(6, 1, 7, 1, hCD, HRR_INT__i_p_l_s, HRR_INT__i_p_k_s, HRR_INT__i_p_k_p);
979 
980                 // form INT__i_p_h_d
981                 ostei_general_hrr_L(6, 1, 5, 2, hCD, HRR_INT__i_p_i_p, HRR_INT__i_p_h_p, HRR_INT__i_p_h_d);
982 
983                 // form INT__i_p_i_d
984                 ostei_general_hrr_L(6, 1, 6, 2, hCD, HRR_INT__i_p_k_p, HRR_INT__i_p_i_p, HRR_INT__i_p_i_d);
985 
986                 // form INT__i_p_h_f
987                 ostei_general_hrr_L(6, 1, 5, 3, hCD, HRR_INT__i_p_i_d, HRR_INT__i_p_h_d, HRR_INT__i_p_h_f);
988 
989 
990             }  // close HRR loop
991 
992 
993         }   // close loop cdbatch
994 
995         istart = iend;
996     }  // close loop over ab
997 
998     return P.nshell12_clip * Q.nshell12_clip;
999 }
1000 
ostei_p_i_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__p_i_h_f)1001 int ostei_p_i_h_f(struct simint_multi_shellpair const P,
1002                   struct simint_multi_shellpair const Q,
1003                   double screen_tol,
1004                   double * const restrict work,
1005                   double * const restrict INT__p_i_h_f)
1006 {
1007     double P_AB[3*P.nshell12];
1008     struct simint_multi_shellpair P_tmp = P;
1009     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
1010     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
1011     P_tmp.AB_x = P_AB;
1012     P_tmp.AB_y = P_AB + P.nshell12;
1013     P_tmp.AB_z = P_AB + 2*P.nshell12;
1014 
1015     for(int i = 0; i < P.nshell12; i++)
1016     {
1017         P_tmp.AB_x[i] = -P.AB_x[i];
1018         P_tmp.AB_y[i] = -P.AB_y[i];
1019         P_tmp.AB_z[i] = -P.AB_z[i];
1020     }
1021 
1022     int ret = ostei_i_p_h_f(P_tmp, Q, screen_tol, work, INT__p_i_h_f);
1023     double buffer[17640] SIMINT_ALIGN_ARRAY_DBL;
1024 
1025     for(int q = 0; q < ret; q++)
1026     {
1027         int idx = 0;
1028         for(int a = 0; a < 3; ++a)
1029         for(int b = 0; b < 28; ++b)
1030         for(int c = 0; c < 21; ++c)
1031         for(int d = 0; d < 10; ++d)
1032             buffer[idx++] = INT__p_i_h_f[q*17640+b*630+a*210+c*10+d];
1033 
1034         memcpy(INT__p_i_h_f+q*17640, buffer, 17640*sizeof(double));
1035     }
1036 
1037     return ret;
1038 }
1039 
ostei_i_p_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__i_p_f_h)1040 int ostei_i_p_f_h(struct simint_multi_shellpair const P,
1041                   struct simint_multi_shellpair const Q,
1042                   double screen_tol,
1043                   double * const restrict work,
1044                   double * const restrict INT__i_p_f_h)
1045 {
1046     double Q_AB[3*Q.nshell12];
1047     struct simint_multi_shellpair Q_tmp = Q;
1048     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
1049     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
1050     Q_tmp.AB_x = Q_AB;
1051     Q_tmp.AB_y = Q_AB + Q.nshell12;
1052     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
1053 
1054     for(int i = 0; i < Q.nshell12; i++)
1055     {
1056         Q_tmp.AB_x[i] = -Q.AB_x[i];
1057         Q_tmp.AB_y[i] = -Q.AB_y[i];
1058         Q_tmp.AB_z[i] = -Q.AB_z[i];
1059     }
1060 
1061     int ret = ostei_i_p_h_f(P, Q_tmp, screen_tol, work, INT__i_p_f_h);
1062     double buffer[17640] SIMINT_ALIGN_ARRAY_DBL;
1063 
1064     for(int q = 0; q < ret; q++)
1065     {
1066         int idx = 0;
1067         for(int a = 0; a < 28; ++a)
1068         for(int b = 0; b < 3; ++b)
1069         for(int c = 0; c < 10; ++c)
1070         for(int d = 0; d < 21; ++d)
1071             buffer[idx++] = INT__i_p_f_h[q*17640+a*630+b*210+d*10+c];
1072 
1073         memcpy(INT__i_p_f_h+q*17640, buffer, 17640*sizeof(double));
1074     }
1075 
1076     return ret;
1077 }
1078 
ostei_p_i_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__p_i_f_h)1079 int ostei_p_i_f_h(struct simint_multi_shellpair const P,
1080                   struct simint_multi_shellpair const Q,
1081                   double screen_tol,
1082                   double * const restrict work,
1083                   double * const restrict INT__p_i_f_h)
1084 {
1085     double P_AB[3*P.nshell12];
1086     struct simint_multi_shellpair P_tmp = P;
1087     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
1088     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
1089     P_tmp.AB_x = P_AB;
1090     P_tmp.AB_y = P_AB + P.nshell12;
1091     P_tmp.AB_z = P_AB + 2*P.nshell12;
1092 
1093     for(int i = 0; i < P.nshell12; i++)
1094     {
1095         P_tmp.AB_x[i] = -P.AB_x[i];
1096         P_tmp.AB_y[i] = -P.AB_y[i];
1097         P_tmp.AB_z[i] = -P.AB_z[i];
1098     }
1099 
1100     double Q_AB[3*Q.nshell12];
1101     struct simint_multi_shellpair Q_tmp = Q;
1102     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
1103     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
1104     Q_tmp.AB_x = Q_AB;
1105     Q_tmp.AB_y = Q_AB + Q.nshell12;
1106     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
1107 
1108     for(int i = 0; i < Q.nshell12; i++)
1109     {
1110         Q_tmp.AB_x[i] = -Q.AB_x[i];
1111         Q_tmp.AB_y[i] = -Q.AB_y[i];
1112         Q_tmp.AB_z[i] = -Q.AB_z[i];
1113     }
1114 
1115     int ret = ostei_i_p_h_f(P_tmp, Q_tmp, screen_tol, work, INT__p_i_f_h);
1116     double buffer[17640] SIMINT_ALIGN_ARRAY_DBL;
1117 
1118     for(int q = 0; q < ret; q++)
1119     {
1120         int idx = 0;
1121         for(int a = 0; a < 3; ++a)
1122         for(int b = 0; b < 28; ++b)
1123         for(int c = 0; c < 10; ++c)
1124         for(int d = 0; d < 21; ++d)
1125             buffer[idx++] = INT__p_i_f_h[q*17640+b*630+a*210+d*10+c];
1126 
1127         memcpy(INT__p_i_f_h+q*17640, buffer, 17640*sizeof(double));
1128     }
1129 
1130     return ret;
1131 }
1132 
1133