diff options
Diffstat (limited to 'sc')
-rw-r--r-- | sc/source/core/opencl/op_statistical.cxx | 13 | ||||
-rw-r--r-- | sc/source/core/opencl/opinlinefun_statistical.cxx | 202 |
2 files changed, 122 insertions, 93 deletions
diff --git a/sc/source/core/opencl/op_statistical.cxx b/sc/source/core/opencl/op_statistical.cxx index 85cd762eda5c..d2e3c391a1e6 100644 --- a/sc/source/core/opencl/op_statistical.cxx +++ b/sc/source/core/opencl/op_statistical.cxx @@ -6743,7 +6743,8 @@ void OpGammaInv::GenSlidingWindowFunction(std::stringstream &ss, ss << "))\n"; ss << " arg"<<i<<"= 0;\n"; ss << " else\n"; - ss << " arg"<<i<<"="<<vSubArguments[i]->GenSlidingWindowDeclRef(); + ss << " arg"<<i<<"="; + ss <<vSubArguments[i]->GenSlidingWindowDeclRef(); ss << ";\n"; ss << " }\n"; ss << " else\n"; @@ -6758,7 +6759,8 @@ void OpGammaInv::GenSlidingWindowFunction(std::stringstream &ss, ss << "))\n"; ss << " arg"<<i<<"= 0;\n"; ss << " else\n"; - ss << " arg"<<i<<"="<<vSubArguments[i]->GenSlidingWindowDeclRef(); + ss << " arg"<<i<<"="; + ss <<vSubArguments[i]->GenSlidingWindowDeclRef(); ss << ";\n"; #endif } @@ -7043,9 +7045,10 @@ void OpFInv::GenSlidingWindowFunction(std::stringstream &ss, " {\n" " if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n" " {\n" - " fSx = fPx * fRy * fQy / (fRy-fPy) / (fQy-fPy)+" - " fRx * fQy * fPy" - "/ (fQy-fRy) / (fPy-fRy)+ fQx * fPy * fRy / (fPy-fQy) / (fRy-fQy);\n" + " fSx = fPx * fRy * fQy *pow( (fRy-fPy),-1)" + " *pow( (fQy-fPy),-1)+fRx * fQy * fPy*pow( (fQy-fRy),-1) *" + "pow( (fPy-fRy),-1)+ fQx * fPy * fRy *pow( (fPy-fQy),-1)" + " *pow((fRy-fQy),-1);\n" " bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n" " }\n" " else\n" diff --git a/sc/source/core/opencl/opinlinefun_statistical.cxx b/sc/source/core/opencl/opinlinefun_statistical.cxx index dce8f6692e07..b75be381469c 100644 --- a/sc/source/core/opencl/opinlinefun_statistical.cxx +++ b/sc/source/core/opencl/opinlinefun_statistical.cxx @@ -252,9 +252,9 @@ std::string GetFInvValueDecl = "double GetFInvValue(double fF1,double fF2" std::string GetFInvValue = "double GetFInvValue(double fF1,double fF2,double fX1 )\n" "{\n" -" double arg = fF2/(fF2+fF1*fX1);\n" -" double alpha = fF2/2.0;\n" -" double beta = fF1/2.0;\n" +" double arg = fF2*pow((fF2+fF1*fX1),-1.0);\n" +" double alpha = fF2*pow(2.0,-1.0);\n" +" double beta = fF1*pow(2.0,-1.0);\n" " double fXin,fAlpha,fBeta;\n" " fXin=arg;fAlpha=alpha;fBeta=beta;\n" " if (fXin <= 0.0)\n" @@ -272,7 +272,7 @@ std::string GetFInvValue = " double flnX = log(fXin);\n" " double fA = fAlpha;\n" " double fB = fBeta;\n" -" bool bReflect = fXin > fAlpha/(fAlpha+fBeta);\n" +" bool bReflect = fXin > fAlpha*pow((fAlpha+fBeta),-1.0);\n" " if (bReflect)\n" " {\n" " fA = fBeta;\n" @@ -283,9 +283,9 @@ std::string GetFInvValue = " flnY = log(fXin);\n" " }\n" " fResult = lcl_GetBetaHelperContFrac(fX,fA,fB);\n" -" fResult = fResult/fA;\n" -" double fP = fA/(fA+fB);\n" -" double fQ = fB/(fA+fB);\n" +" fResult = fResult*pow(fA,-1.0);\n" +" double fP = fA*pow((fA+fB),-1.0);\n" +" double fQ = fB*pow((fA+fB),-1.0);\n" " double fTemp;\n" " if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97)\n" " fTemp = GetBetaDistPDF(fX,fA,fB)*fX*fY;\n" @@ -951,30 +951,88 @@ std::string lcl_getLanczosSum = " int nI;\n" " if (fZ<=1.0)\n" " {\n" -" fSumNum = fNum[12];\n" -" fSumDenom = fDenom[12];\n" -" for (nI = 11; nI >= 0; --nI)\n" -" {\n" -" fSumNum *= fZ;\n" -" fSumNum += fNum[nI];\n" -" fSumDenom *= fZ;\n" -" fSumDenom += fDenom[nI];\n" -" }\n" +" fSumNum = fNum[12];\n" +" fSumDenom = fDenom[12];\n" +" nI = 11;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 10;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 9;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 8;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 7;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 6;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 5;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 4;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 3;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 2;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 1;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 0;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" " }\n" -" else\n" +" if (fZ>1.0)\n" " {\n" -" double fZInv = 1/fZ;\n" -" fSumNum = fNum[0];\n" -" fSumDenom = fDenom[0];\n" -" for (nI = 1; nI <=12; ++nI)\n" -" {\n" -" fSumNum *= fZInv;\n" -" fSumNum += fNum[nI];\n" -" fSumDenom *= fZInv;\n" -" fSumDenom += fDenom[nI];\n" -" }\n" +" double fZInv = pow(fZ,-1.0);\n" +" fSumNum = fNum[0];\n" +" fSumDenom = fDenom[0];\n" +" nI = 1;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 2;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 3;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 4;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 5;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 6;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 7;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 8;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 9;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 10;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 11;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 12;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" " }\n" -" return fSumNum/fSumDenom;\n" +" return fSumNum*pow(fSumDenom,-1.0);\n" "}\n"; std::string GetUpRegIGammaDecl= @@ -1011,28 +1069,17 @@ std::string GetBeta = "{\n" " double fA;\n" " double fB;\n" -" if (fAlpha > fBeta)\n" -" {\n" -" fA = fAlpha; fB = fBeta;\n" -" }\n" -" else\n" -" {\n" -" fA = fBeta; fB = fAlpha;\n" -" }\n" -" if (fA+fB < fMaxGammaArgument)\n" -" return tgamma(fA)/tgamma(fA + fB)*tgamma(fB);\n" -" double fg = 6.024680040776729583740234375;\n" -" double fgm = fg - 0.5;\n" -" double fLanczos = lcl_getLanczosSum(fA);\n" -" fLanczos /= lcl_getLanczosSum(fA + fB);\n" -" fLanczos *= lcl_getLanczosSum(fB);\n" -" double fABgm = fA + fB + fgm;\n" -" fLanczos *= sqrt((fABgm/(fA + fgm))/(fB + fgm));\n" -" double fTempA = fB/(fA + fgm);\n" -" double fTempB = fA/(fB + fgm);\n" -" double fResult = exp(-fA*log1p(fTempA) - fB*log1p(fTempB) - fgm);\n" -" fResult *= fLanczos;\n" -" return fResult;\n" +" fAlpha>fBeta?(fA = fAlpha,fB = fBeta):(fA = fBeta,fB = fAlpha);\n" +" double fAB = fA + fB;\n" + +" if (fAB < fMaxGammaArgument)\n" +" return tgamma(fA)*pow(tgamma(fAB),-1.0)*tgamma(fB);\n" +" double fgm = 5.524680040776729583740234375;\n" +" double fLanczos = lcl_getLanczosSum(fA)*lcl_getLanczosSum(fB)\n" +" *pow(lcl_getLanczosSum(fAB),-1.0);\n" +" fLanczos *= sqrt(((fAB + fgm)*pow(fA + fgm,-1.0))*pow(fB + fgm,-1.0));\n" +" return fLanczos * exp(-fA*log1p(fB*pow(fA + fgm,-1.0))\n" +" - fB*log1p(fA*pow(fB + fgm,-1.0)) - fgm);\n" "}\n"; std::string GetLogBetaDecl= @@ -1042,27 +1089,15 @@ std::string GetLogBeta = "{\n" " double fA;\n" " double fB;\n" -" if (fAlpha > fBeta)\n" -" {\n" -" fA = fAlpha; fB = fBeta;\n" -" }\n" -" else\n" -" {\n" -" fA = fBeta; fB = fAlpha;\n" -" }\n" -" double fg = 6.024680040776729583740234375;\n" -" double fgm = fg - 0.5;\n" -" double fLanczos = lcl_getLanczosSum(fA);\n" -" fLanczos /= lcl_getLanczosSum(fA + fB);\n" -" fLanczos *= lcl_getLanczosSum(fB);\n" -" double fLogLanczos = log(fLanczos);\n" -" double fABgm = fA + fB + fgm;\n" -" fLogLanczos += 0.5*(log(fABgm) - log(fA + fgm) - log(fB + fgm));\n" -" double fTempA = fB/(fA + fgm);\n" -" double fTempB = fA/(fB + fgm);\n" -" double fResult = -fA * log1p(fTempA)\n" -" -fB * log1p(fTempB)-fgm;\n" -" fResult += fLogLanczos;\n" +" fAlpha>fBeta?(fA = fAlpha,fB = fBeta):(fA = fBeta,fB = fAlpha);\n" +" double fgm = 5.524680040776729583740234375;\n" + +" double fLanczos = lcl_getLanczosSum(fA)*lcl_getLanczosSum(fB)*\n" +" pow(lcl_getLanczosSum(fA + fB),-1.0);\n" +" double fResult= -fA *log1p(fB*pow(fA + fgm,-1.0))" +"-fB *log1p(fA*pow(fB + fgm,-1.0))-fgm;\n" +" fResult += log(fLanczos)+0.5*(log(fA + fB + fgm) - log(fA + fgm)\n" +" - log(fB + fgm));\n" " return fResult;\n" "}\n"; @@ -1126,7 +1161,8 @@ std::string GetBetaDistPDF = " && fLogBeta < fLogDblMax && fLogBeta > fLogDblMin\n" " && fAm1LogX + fBm1LogY < fLogDblMax && fAm1LogX + fBm1LogY > \n" " fLogDblMin)\n" -" return pow(fX,fA-1.0) * pow(0.5-fX+0.5,fB-1.0) / GetBeta(fA,fB);\n" +" return pow(fX,fA-1.0)*pow(0.5-fX+0.5,fB-1.0)" +"*pow(GetBeta(fA,fB),-1.0);\n" " else \n" " return exp( fAm1LogX + fBm1LogY - fLogBeta);\n" "}\n"; @@ -1139,19 +1175,9 @@ std::string lcl_GetBetaHelperContFrac = " double a1, b1, a2, b2, fnorm, apl2m, d2m, d2m1, cfnew, cf;\n" " a1 = 1.0; b1 = 1.0;\n" -" b2 = 1.0 - (fA+fB)/(fA+1.0)*fX;\n" -" if (b2 == 0.0)\n" -" {\n" -" a2 = 0.0;\n" -" fnorm = 1.0;\n" -" cf = 1.0;\n" -" }\n" -" else\n" -" {\n" -" a2 = 1.0;\n" -" fnorm = 1.0/b2;\n" -" cf = a2*fnorm;\n" -" }\n" +" b2 = 1.0 - (fA+fB)*pow(fA+1.0,-1.0)*fX;\n" +" b2==0.0?(a2 = 0.0,fnorm = 1.0,cf = 1.0):\n" +" (a2 = 1.0,fnorm = pow(b2,-1.0),cf = a2*fnorm);\n" " cfnew = 1.0;\n" " double rm = 1.0;\n" " double fMaxIter = 50000.0;\n" @@ -1159,15 +1185,15 @@ std::string lcl_GetBetaHelperContFrac = " do\n" " {\n" " apl2m = fA + 2.0*rm;\n" -" d2m = rm*(fB-rm)*fX/((apl2m-1.0)*apl2m);\n" -" d2m1 = -(fA+rm)*(fA+fB+rm)*fX/(apl2m*(apl2m+1.0));\n" +" d2m = (rm*(fB-rm))*fX*pow(apl2m*(apl2m-1.0),-1.0);\n" +" d2m1 = -((fA+rm)*(fA+rm+fB))*fX*pow(apl2m*(apl2m+1.0),-1.0);\n" " a1 = (a2+d2m*a1)*fnorm;\n" " b1 = (b2+d2m*b1)*fnorm;\n" " a2 = a1 + d2m1*a2*fnorm;\n" " b2 = b1 + d2m1*b2*fnorm;\n" " if (b2 != 0.0) \n" " {\n" -" fnorm = 1.0/b2;\n" +" fnorm = pow(b2,-1.0);\n" " cfnew = a2*fnorm;\n" " bfinished = (fabs(cf-cfnew) < fabs(cf)*fMachEps);\n" " }\n" |