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