Line data Source code
1 : /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 : /*
3 : * This file is part of the LibreOffice project.
4 : *
5 : * This Source Code Form is subject to the terms of the Mozilla Public
6 : * License, v. 2.0. If a copy of the MPL was not distributed with this
7 : * file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 : */
9 :
10 : #include "op_addin.hxx"
11 :
12 : #include "formulagroup.hxx"
13 : #include "document.hxx"
14 : #include "formulacell.hxx"
15 : #include "tokenarray.hxx"
16 : #include "compiler.hxx"
17 : #include "interpre.hxx"
18 : #include "formula/vectortoken.hxx"
19 : #include <sstream>
20 :
21 : using namespace formula;
22 :
23 : namespace sc { namespace opencl {
24 :
25 0 : void OpBesselj::GenSlidingWindowFunction(std::stringstream &ss,
26 : const std::string &sSymName, SubArguments &vSubArguments)
27 : {
28 0 : ss << "\ndouble " << sSymName;
29 0 : ss << "_" << BinFuncName() << "(";
30 0 : for (unsigned i = 0; i < vSubArguments.size(); i++)
31 : {
32 0 : if (i)
33 0 : ss << ",";
34 0 : vSubArguments[i]->GenSlidingWindowDecl(ss);
35 : }
36 0 : ss << ") {\n";
37 0 : ss << " int gid0 = get_global_id(0);\n";
38 0 : ss << " double x = 0.0;\n";
39 0 : ss << " double N = 0.0;\n";
40 0 : if(vSubArguments.size() != 2)
41 : {
42 0 : ss << " return DBL_MAX;\n" << "}\n";
43 0 : return ;
44 : }
45 0 : FormulaToken *tmpCur0 = vSubArguments[0]->GetFormulaToken();
46 : assert(tmpCur0);
47 0 : if(ocPush == vSubArguments[0]->GetFormulaToken()->GetOpCode())
48 : {
49 0 : if(tmpCur0->GetType() == formula::svSingleVectorRef)
50 : {
51 : const formula::SingleVectorRefToken*tmpCurSVR0 =
52 0 : static_cast<const formula::SingleVectorRefToken *>(tmpCur0);
53 : #ifdef ISNAN
54 0 : ss << " if (gid0 < " << tmpCurSVR0->GetArrayLength() << ")\n";
55 0 : ss << " {\n";
56 : #endif
57 0 : ss << " x = ";
58 0 : ss << vSubArguments[0]->GenSlidingWindowDeclRef() << ";\n";
59 : #ifdef ISNAN
60 0 : ss << " if (isNan(x))\n";
61 0 : ss << " x = 0.0;\n";
62 0 : ss << " }\n";
63 : #endif
64 : }
65 0 : else if(tmpCur0->GetType() == formula::svDouble)
66 : {
67 0 : ss << " x = " << tmpCur0->GetDouble() << ";\n";
68 : }
69 : else
70 : {
71 0 : ss << " return DBL_MAX;\n" << "}\n";
72 0 : return ;
73 : }
74 : }
75 : else
76 : {
77 0 : ss << " x = ";
78 0 : ss << vSubArguments[0]->GenSlidingWindowDeclRef() << ";\n";
79 : }
80 :
81 0 : FormulaToken *tmpCur1 = vSubArguments[1]->GetFormulaToken();
82 : assert(tmpCur1);
83 0 : if(ocPush == vSubArguments[1]->GetFormulaToken()->GetOpCode())
84 : {
85 0 : if(tmpCur1->GetType() == formula::svSingleVectorRef)
86 : {
87 : const formula::SingleVectorRefToken*tmpCurSVR1 =
88 0 : static_cast<const formula::SingleVectorRefToken *>(tmpCur1);
89 : #ifdef ISNAN
90 0 : ss << " if (gid0 < " << tmpCurSVR1->GetArrayLength() << ")\n";
91 0 : ss << " {\n";
92 : #endif
93 0 : ss << " N = ";
94 0 : ss << vSubArguments[1]->GenSlidingWindowDeclRef() << ";\n";
95 : #ifdef ISNAN
96 0 : ss << " if (isNan(N))\n";
97 0 : ss << " N = 0.0;\n";
98 0 : ss << " }\n";
99 : #endif
100 : }
101 0 : else if(tmpCur1->GetType() == formula::svDouble)
102 : {
103 0 : ss << " N = " << tmpCur1->GetDouble() << ";\n";
104 : }
105 : else
106 : {
107 0 : ss << " return DBL_MAX;\n" << "}\n";
108 0 : return ;
109 : }
110 : }
111 : else
112 : {
113 0 : ss << " N = ";
114 0 : ss << vSubArguments[1]->GenSlidingWindowDeclRef() << ";\n";
115 : }
116 0 : ss << " double f_PI = 3.1415926535897932385;\n";
117 0 : ss << " double f_2_DIV_PI = 2.0 / f_PI;\n";
118 0 : ss << " double f_PI_DIV_2 = f_PI / 2.0;\n";
119 0 : ss << " double f_PI_DIV_4 = f_PI / 4.0;\n";
120 0 : ss << " if( N < 0.0 )\n";
121 0 : ss << " return DBL_MAX;\n";
122 0 : ss << " if (x == 0.0)\n";
123 0 : ss << " return (N == 0.0) ? 1.0 : 0.0;\n";
124 0 : ss << " double fSign = ((int)N % 2 == 1 && x < 0.0) ? -1.0 : 1.0;\n";
125 0 : ss << " double fX = fabs(x);\n";
126 0 : ss << " double fMaxIteration = 9000000.0;\n";
127 0 : ss << " double fEstimateIteration = fX * 1.5 + N;\n";
128 0 : ss << " bool bAsymptoticPossible = pow(fX,0.4) > N;\n";
129 0 : ss << " if (fEstimateIteration > fMaxIteration)\n";
130 0 : ss << " {\n";
131 0 : ss << " if (bAsymptoticPossible)\n";
132 0 : ss << " return fSign * sqrt(f_2_DIV_PI/fX)";
133 0 : ss << "* cos(fX-N*f_PI_DIV_2-f_PI_DIV_4);\n";
134 0 : ss << " else\n";
135 0 : ss << " return DBL_MAX;\n";
136 0 : ss << " }\n";
137 0 : ss << " double epsilon = 1.0e-15;\n";
138 0 : ss << " bool bHasfound = false;\n";
139 0 : ss << " double k= 0.0;\n";
140 0 : ss << " double u ;\n";
141 0 : ss << " double m_bar;\n";
142 0 : ss << " double g_bar;\n";
143 0 : ss << " double g_bar_delta_u;\n";
144 0 : ss << " double g = 0.0;\n";
145 0 : ss << " double delta_u = 0.0;\n";
146 0 : ss << " double f_bar = -1.0;\n";
147 0 : ss << " if (N==0)\n";
148 0 : ss << " {\n";
149 0 : ss << " u = 1.0;\n";
150 0 : ss << " g_bar_delta_u = 0.0;\n";
151 0 : ss << " g_bar = - 2.0/fX; \n";
152 0 : ss << " delta_u = g_bar_delta_u / g_bar;\n";
153 0 : ss << " u = u + delta_u ;\n";
154 0 : ss << " g = -1.0 / g_bar; \n";
155 0 : ss << " f_bar = f_bar * g;\n";
156 0 : ss << " k = 2.0;\n";
157 0 : ss << " }\n";
158 0 : ss << " if (N!=0)\n";
159 0 : ss << " {\n";
160 0 : ss << " u=0.0;\n";
161 0 : ss << " for (k =1.0; k<= N-1; k = k + 1.0)\n";
162 0 : ss << " {\n";
163 0 : ss << " m_bar=2.0 * fmod(k-1.0, 2.0) * f_bar;\n";
164 0 : ss << " g_bar_delta_u = - g * delta_u - m_bar * u;\n";
165 0 : ss << " g_bar = m_bar - 2.0*k/fX + g;\n";
166 0 : ss << " delta_u = g_bar_delta_u / g_bar;\n";
167 0 : ss << " u = u + delta_u;\n";
168 0 : ss << " g = -1.0/g_bar;\n";
169 0 : ss << " f_bar=f_bar * g;\n";
170 0 : ss << " }\n";
171 0 : ss << " m_bar=2.0 * fmod(k-1.0, 2.0) * f_bar;\n";
172 0 : ss << " g_bar_delta_u = f_bar - g * delta_u - m_bar * u;\n";
173 0 : ss << " g_bar = m_bar - 2.0*k/fX + g;\n";
174 0 : ss << " delta_u = g_bar_delta_u / g_bar;\n";
175 0 : ss << " u = u + delta_u;\n";
176 0 : ss << " g = -1.0/g_bar;\n";
177 0 : ss << " f_bar = f_bar * g;\n";
178 0 : ss << " k = k + 1.0;\n";
179 0 : ss << " }\n";
180 0 : ss << " do\n";
181 0 : ss << " {\n";
182 0 : ss << " m_bar = 2.0 * fmod(k-1.0, 2.0) * f_bar;\n";
183 0 : ss << " g_bar_delta_u = - g * delta_u - m_bar * u;\n";
184 0 : ss << " g_bar = m_bar - 2.0*k/fX + g;\n";
185 0 : ss << " delta_u = g_bar_delta_u / g_bar;\n";
186 0 : ss << " u = u + delta_u;\n";
187 0 : ss << " g = -pow(g_bar,-1.0);\n";
188 0 : ss << " f_bar = f_bar * g;\n";
189 0 : ss << " bHasfound = (fabs(delta_u)<=fabs(u)*epsilon);\n";
190 0 : ss << " k = k + 1.0;\n";
191 0 : ss << " }\n";
192 0 : ss << " while (!bHasfound && k <= fMaxIteration);\n";
193 0 : ss << " if (bHasfound)\n";
194 0 : ss << " return u * fSign;\n";
195 0 : ss << " else\n";
196 0 : ss << " return DBL_MAX;\n";
197 0 : ss << "}";
198 : }
199 :
200 : }}
201 :
202 : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|
