sc/source/core/opencl/op_addin.cxx

   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 void OpBesselj::GenSlidingWindowFunction(std::stringstream &ss,
  26     const std::string sSymName, SubArguments &vSubArguments)
  27 {
  28     ss << "\ndouble " << sSymName;
  29     ss << "_" << BinFuncName() << "(";
  30     for (unsigned i = 0; i < vSubArguments.size(); i++)
  31     {
  32         if (i)
  33             ss << ",";
  34         vSubArguments[i]->GenSlidingWindowDecl(ss);
  35     }
  36     ss << ") {\n";
  37     ss << "    int gid0 = get_global_id(0);\n";
  38     ss << "    double x = 0.0;\n";
  39     ss << "    double N = 0.0;\n";
  40     if(vSubArguments.size() != 2)
  41     {
  42         ss << "    return DBL_MAX;\n" << "}\n";
  43         return ;
  44     }
  45     FormulaToken *tmpCur0 = vSubArguments[0]->GetFormulaToken();
  46     assert(tmpCur0);
  47     if(ocPush == vSubArguments[0]->GetFormulaToken()->GetOpCode())
  48     {
  49         if(tmpCur0->GetType() == formula::svSingleVectorRef)
  50         {
  51             const formula::SingleVectorRefToken*tmpCurSVR0 =
  52                 dynamic_cast<const formula::SingleVectorRefToken *>(tmpCur0);
  53 #ifdef  ISNAN
  54             ss << "    if (gid0 < " << tmpCurSVR0->GetArrayLength() << ")\n";
  55             ss << "    {\n";
  56 #endif
  57             ss << "        x = ";
  58             ss << vSubArguments[0]->GenSlidingWindowDeclRef() << ";\n";
  59 #ifdef  ISNAN
  60             ss << "        if (isNan(x))\n";
  61             ss << "            x = 0.0;\n";
  62             ss << "    }\n";
  63 #endif
  64         }
  65         else if(tmpCur0->GetType() == formula::svDouble)
  66         {
  67             ss << "    x = " << tmpCur0->GetDouble() << ";\n";
  68         }
  69         else
  70         {
  71             ss << "    return DBL_MAX;\n" << "}\n";
  72             return ;
  73         }
  74     }
  75     else
  76     {
  77         ss << "    x = ";
  78         ss << vSubArguments[0]->GenSlidingWindowDeclRef() << ";\n";
  79     }
  80
  81     FormulaToken *tmpCur1 = vSubArguments[1]->GetFormulaToken();
  82     assert(tmpCur1);
  83     if(ocPush == vSubArguments[1]->GetFormulaToken()->GetOpCode())
  84     {
  85         if(tmpCur1->GetType() == formula::svSingleVectorRef)
  86         {
  87             const formula::SingleVectorRefToken*tmpCurSVR1 =
  88                 dynamic_cast<const formula::SingleVectorRefToken *>(tmpCur1);
  89 #ifdef  ISNAN
  90             ss << "    if (gid0 < " << tmpCurSVR1->GetArrayLength() << ")\n";
  91             ss << "    {\n";
  92 #endif
  93             ss << "        N = ";
  94             ss << vSubArguments[1]->GenSlidingWindowDeclRef() << ";\n";
  95 #ifdef  ISNAN
  96             ss << "        if (isNan(N))\n";
  97             ss << "            N = 0.0;\n";
  98             ss << "    }\n";
  99 #endif
 100         }
 101         else if(tmpCur1->GetType() == formula::svDouble)
 102         {
 103             ss << "    N = " << tmpCur1->GetDouble() << ";\n";
 104         }
 105         else
 106         {
 107             ss << "    return DBL_MAX;\n" << "}\n";
 108             return ;
 109         }
 110     }
 111     else
 112     {
 113         ss << "    N = ";
 114         ss << vSubArguments[1]->GenSlidingWindowDeclRef() << ";\n";
 115     }
 116     ss << "    double f_PI       = 3.1415926535897932385;\n";
 117     ss << "    double f_2_DIV_PI = 2.0 / f_PI;\n";
 118     ss << "    double f_PI_DIV_2 = f_PI / 2.0;\n";
 119     ss << "    double f_PI_DIV_4 = f_PI / 4.0;\n";
 120     ss << "    if( N < 0.0 )\n";
 121     ss << "        return DBL_MAX;\n";
 122     ss << "    if (x == 0.0)\n";
 123     ss << "        return (N == 0.0) ? 1.0 : 0.0;\n";
 124     ss << "    double fSign = ((int)N % 2 == 1 && x < 0.0) ? -1.0 : 1.0;\n";
 125     ss << "    double fX = fabs(x);\n";
 126     ss << "    double fMaxIteration = 9000000.0;\n";
 127     ss << "    double fEstimateIteration = fX * 1.5 + N;\n";
 128     ss << "    bool bAsymptoticPossible = pow(fX,0.4) > N;\n";
 129     ss << "    if (fEstimateIteration > fMaxIteration)\n";
 130     ss << "    {\n";
 131     ss << "        if (bAsymptoticPossible)\n";
 132     ss << "            return fSign * sqrt(f_2_DIV_PI/fX)";
 133     ss << "* cos(fX-N*f_PI_DIV_2-f_PI_DIV_4);\n";
 134     ss << "        else\n";
 135     ss << "            return DBL_MAX;\n";
 136     ss << "    }\n";
 137     ss << "    double epsilon = 1.0e-15;\n";
 138     ss << "    bool bHasfound = false;\n";
 139     ss << "    double k= 0.0;\n";
 140     ss << "    double  u ;\n";
 141     ss << "    double m_bar;\n";
 142     ss << "    double g_bar;\n";
 143     ss << "    double g_bar_delta_u;\n";
 144     ss << "    double g = 0.0;\n";
 145     ss << "    double delta_u = 0.0;\n";
 146     ss << "    double f_bar = -1.0;\n";
 147     ss << "    if (N==0)\n";
 148     ss << "    {\n";
 149     ss << "        u = 1.0;\n";
 150     ss << "        g_bar_delta_u = 0.0;\n";
 151     ss << "        g_bar = - 2.0/fX; \n";
 152     ss << "        delta_u = g_bar_delta_u / g_bar;\n";
 153     ss << "        u = u + delta_u ;\n";
 154     ss << "        g = -1.0 / g_bar; \n";
 155     ss << "        f_bar = f_bar * g;\n";
 156     ss << "        k = 2.0;\n";
 157     ss << "    }\n";
 158     ss << "    if (N!=0)\n";
 159     ss << "    {\n";
 160     ss << "        u=0.0;\n";
 161     ss << "        for (k =1.0; k<= N-1; k = k + 1.0)\n";
 162     ss << "        {\n";
 163     ss << "            m_bar=2.0 * fmod(k-1.0, 2.0) * f_bar;\n";
 164     ss << "            g_bar_delta_u = - g * delta_u - m_bar * u;\n";
 165     ss << "            g_bar = m_bar - 2.0*k/fX + g;\n";
 166     ss << "            delta_u = g_bar_delta_u / g_bar;\n";
 167     ss << "            u = u + delta_u;\n";
 168     ss << "            g = -1.0/g_bar;\n";
 169     ss << "            f_bar=f_bar * g;\n";
 170     ss << "        }\n";
 171     ss << "        m_bar=2.0 * fmod(k-1.0, 2.0) * f_bar;\n";
 172     ss << "        g_bar_delta_u = f_bar - g * delta_u - m_bar * u;\n";
 173     ss << "        g_bar = m_bar - 2.0*k/fX + g;\n";
 174     ss << "        delta_u = g_bar_delta_u / g_bar;\n";
 175     ss << "        u = u + delta_u;\n";
 176     ss << "        g = -1.0/g_bar;\n";
 177     ss << "        f_bar = f_bar * g;\n";
 178     ss << "        k = k + 1.0;\n";
 179     ss << "    }\n";
 180     ss << "    do\n";
 181     ss << "    {\n";
 182     ss << "        m_bar = 2.0 * fmod(k-1.0, 2.0) * f_bar;\n";
 183     ss << "        g_bar_delta_u = - g * delta_u - m_bar * u;\n";
 184     ss << "        g_bar = m_bar - 2.0*k/fX + g;\n";
 185     ss << "        delta_u = g_bar_delta_u / g_bar;\n";
 186     ss << "        u = u + delta_u;\n";
 187     ss << "        g = -pow(g_bar,-1.0);\n";
 188     ss << "        f_bar = f_bar * g;\n";
 189     ss << "        bHasfound = (fabs(delta_u)<=fabs(u)*epsilon);\n";
 190     ss << "        k = k + 1.0;\n";
 191     ss << "    }\n";
 192     ss << "    while (!bHasfound && k <= fMaxIteration);\n";
 193     ss << "    if (bHasfound)\n";
 194     ss << "        return u * fSign;\n";
 195     ss << "    else\n";
 196     ss << "        return DBL_MAX;\n";
 197     ss << "}";
 198 }
 199
 200 }}
 201
 202 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */