1 #pragma once
2
3 #include "simint/boys/boys_shortgrid.h"
4 #include "simint/vectorization/vectorization.h"
5
6 #ifdef __cplusplus
7 #include "simint/cpp_restrict.hpp"
8 extern "C" {
9 #endif
10
11 extern double boys_shortgrid[BOYS_SHORTGRID_NPOINT][BOYS_SHORTGRID_MAXN+1];
12
13 static inline
boys_F_taylor(double * restrict F,double x,int n)14 void boys_F_taylor(double * restrict F, double x, int n)
15 {
16
17 const int lookup_idx = (int)(BOYS_SHORTGRID_LOOKUPFAC*(x+BOYS_SHORTGRID_LOOKUPFAC2));
18 const double xi = ((double)lookup_idx * BOYS_SHORTGRID_SPACE);
19 const double dx = xi-x; // -delta x
20
21 double const * restrict gridpts = &(boys_shortgrid[lookup_idx][0]);
22
23 #ifndef SIMINT_SCALAR
24 #pragma omp simd
25 #endif
26 for(int i = 0; i <= n; ++i)
27 {
28 double const * restrict gridpts2 = gridpts + i;
29
30 F[i*SIMINT_SIMD_LEN] = gridpts2[0]
31 + dx * ( gridpts2[1]
32 + dx * ( (1.0/2.0 ) * gridpts2[2]
33 + dx * ( (1.0/6.0 ) * gridpts2[3]
34 + dx * ( (1.0/24.0 ) * gridpts2[4]
35 + dx * ( (1.0/120.0 ) * gridpts2[5]
36 + dx * ( (1.0/720.0 ) * gridpts2[6]
37 + dx * ( (1.0/5040.0) * gridpts2[7]
38 )))))));
39 }
40 }
41
42 static inline
boys_F_taylor_single(double x,int n)43 double boys_F_taylor_single(double x, int n)
44 {
45 const int lookup_idx = (int)(BOYS_SHORTGRID_LOOKUPFAC*(x+BOYS_SHORTGRID_LOOKUPFAC2));
46 const double xi = ((double)lookup_idx * BOYS_SHORTGRID_SPACE);
47 const double dx = xi-x; // -delta x
48
49 double const * restrict gridpts = &(boys_shortgrid[lookup_idx][n]);
50
51 return gridpts[0]
52 + dx * ( gridpts[1]
53 + dx * ( (1.0/2.0 ) * gridpts[2]
54 + dx * ( (1.0/6.0 ) * gridpts[3]
55 + dx * ( (1.0/24.0 ) * gridpts[4]
56 + dx * ( (1.0/120.0 ) * gridpts[5]
57 + dx * ( (1.0/720.0 ) * gridpts[6]
58 + dx * ( (1.0/5040.0) * gridpts[7]
59 )))))));
60 }
61
62 static inline
boys_F_taylor_vec(SIMINT_DBLTYPE * restrict F,SIMINT_DBLTYPE x,int n)63 void boys_F_taylor_vec(SIMINT_DBLTYPE * restrict F, SIMINT_DBLTYPE x, int n)
64 {
65 double * restrict Fd = (double *)F;
66 double * restrict xd = (double *)&x;
67
68 #ifndef SIMINT_SCALAR
69 // workaround for GCC missing simdlen??
70 #if defined __clang__ || defined __INTEL_COMPILER
71 #pragma omp simd simdlen(SIMINT_SIMD_LEN)
72 #else
73 #pragma omp simd
74 #endif
75 #endif
76 for(int v = 0; v < SIMINT_SIMD_LEN; v++)
77 {
78 const int lookup_idx = (int)(BOYS_SHORTGRID_LOOKUPFAC*(xd[v]+BOYS_SHORTGRID_LOOKUPFAC2));
79 const double xi = ((double)lookup_idx * BOYS_SHORTGRID_SPACE);
80 const double dx = xi-xd[v]; // -delta x
81
82 double const * restrict gridpts = &(boys_shortgrid[lookup_idx][0]);
83
84 for(int i = 0; i <= n; ++i)
85 {
86 double const * restrict gridpts2 = gridpts + i;
87
88 Fd[i*SIMINT_SIMD_LEN + v] = gridpts2[0]
89 + dx * ( gridpts2[1]
90 + dx * ( (1.0/2.0 ) * gridpts2[2]
91 + dx * ( (1.0/6.0 ) * gridpts2[3]
92 + dx * ( (1.0/24.0 ) * gridpts2[4]
93 + dx * ( (1.0/120.0 ) * gridpts2[5]
94 + dx * ( (1.0/720.0 ) * gridpts2[6]
95 + dx * ( (1.0/5040.0) * gridpts2[7]
96 )))))));
97 }
98 }
99 }
100
101 static inline
boys_F_taylor_single_vec(SIMINT_DBLTYPE x,int n)102 SIMINT_DBLTYPE boys_F_taylor_single_vec(SIMINT_DBLTYPE x, int n)
103 {
104 SIMINT_DBLTYPE ret;
105 double * retd = (double *)&ret;
106 double * restrict xd = (double *)&x;
107
108 #ifndef SIMINT_SCALAR
109 // workaround for GCC missing simdlen??
110 #if defined __clang__ || defined __INTEL_COMPILER
111 #pragma omp simd simdlen(SIMINT_SIMD_LEN)
112 #else
113 #pragma omp simd
114 #endif
115 #endif
116 for(int v = 0; v < SIMINT_SIMD_LEN; v++)
117 {
118 const int lookup_idx = (int)(BOYS_SHORTGRID_LOOKUPFAC*(xd[v]+BOYS_SHORTGRID_LOOKUPFAC2));
119 const double xi = ((double)lookup_idx * BOYS_SHORTGRID_SPACE);
120 const double dx = xi-xd[v]; // -delta x
121
122 double const * restrict gridpts = &(boys_shortgrid[lookup_idx][n]);
123
124 retd[v] = gridpts[0]
125 + dx * ( gridpts[1]
126 + dx * ( (1.0/2.0 ) * gridpts[2]
127 + dx * ( (1.0/6.0 ) * gridpts[3]
128 + dx * ( (1.0/24.0 ) * gridpts[4]
129 + dx * ( (1.0/120.0 ) * gridpts[5]
130 + dx * ( (1.0/720.0 ) * gridpts[6]
131 + dx * ( (1.0/5040.0) * gridpts[7]
132 )))))));
133 }
134
135 return ret;
136 }
137
138
139 #ifdef __cplusplus
140 }
141 #endif
142
143