? i32833_normsinv.diff
Index: inc/interpre.hxx
===================================================================
RCS file: /cvs/sc/sc/source/core/inc/interpre.hxx,v
retrieving revision 1.20.4.2
diff -u -r1.20.4.2 interpre.hxx
--- inc/interpre.hxx	18 May 2005 11:15:39 -0000	1.20.4.2
+++ inc/interpre.hxx	22 Jul 2005 16:41:44 -0000
@@ -612,6 +612,8 @@
 double taylor(double* pPolynom, USHORT nMax, double x);
 double gauss(double x);
 double gaussinv(double x);
+double normsdist(double x);
+double normsinv(double x);
 double GetGammaDist(double x, double alpha, double beta);
 double GetBetaDist(double x, double alpha, double beta);
 double GetChiDist(double fChi, double fDF);
Index: tool/interpr3.cxx
===================================================================
RCS file: /cvs/sc/sc/source/core/tool/interpr3.cxx,v
retrieving revision 1.12.136.1
diff -u -r1.12.136.1 interpr3.cxx
--- tool/interpr3.cxx	20 Jul 2005 15:20:22 -0000	1.12.136.1
+++ tool/interpr3.cxx	22 Jul 2005 16:41:45 -0000
@@ -59,12 +59,6 @@
  *
  ************************************************************************/
 
-#ifdef PCH
-#include "core_pch.hxx"
-#endif
-
-#pragma hdrstop
-
 // INCLUDE ---------------------------------------------------------------
 
 #include <tools/solar.h>
@@ -283,6 +277,79 @@
 #pragma optimize("",on)
 #endif
 
+
+double ScInterpreter::normsdist(double x)
+{
+	return 0.5 + gauss(x);
+}
+
+
+double ScInterpreter::normsinv(double x)
+{
+    double y;
+	if (x < 0.0 || x > 1.0)
+    {
+	    SetError(errIllegalArgument);
+        y = 0.0;
+    }
+	else if (x == 0.0 || x == 1.0)
+    {
+        SetError(errNoValue);
+        y = 0.0;
+    }
+	else
+    {
+        // Algorithm of #i32833#
+        //
+        // We have
+        // normsdist(-5.99) =~= 0,000000001049
+        // normsdist(5,99)  =~= 0,999999998951
+
+        /* Outside the external borders of our searching segment, we use the
+         * old method. We know the result will be with an accuracy better than
+         * 1.10E-8. Inside, we use the new approach.
+         */
+        if (x <= 0.000000001049 || x >= 0.999999998951)
+		    y = gaussinv(x);
+        else
+        {
+            const double fEps = 1.10E-12;
+            // set the searching segment
+            double fMin = -5.99;
+            double fMax =  5.99;
+            double fMid;
+            do
+            {
+                fMid = (fMin + fMax) / 2;
+                double my_y = normsdist( fMid);
+                if (::rtl::math::approxEqual( x, my_y))
+                {
+                    // jackpot, we're lucky, stop the algorithm
+                    // y = fMid;   is done below after the loop
+                    break;  // while
+                }
+                else
+                {
+                    if (x < my_y )
+                        fMax = fMid;
+                    else
+                        fMin = fMid;
+                }
+            } while (fabs( fMin - fMax) > fEps);
+
+            /* End of loop --> we know that  y  is in the segment
+             * [ y - fMin-fEps, y + fMin-fEps] (with our definition of
+             * the normal standard law normsdist, not the real law). It means
+             * that we're sure that normsdist(norminv(x)) and x have the same 7
+             * decimals, and probably more.
+             */
+            y = fMid;
+        }
+    }
+	return y;
+}
+
+
 double ScInterpreter::Fakultaet(double x)
 {
 	x = ::rtl::math::approxFloor(x);
@@ -975,7 +1042,7 @@
 
 void ScInterpreter::ScStdNormDist()
 {
-	PushDouble(0.5 + gauss(GetDouble()));
+	PushDouble( normsdist( GetDouble()));
 }
 
 void ScInterpreter::ScExpDist()
@@ -1180,37 +1247,41 @@
 		double sigma = GetDouble();
 		double mue   = GetDouble();
 		double x     = GetDouble();
-		if (sigma <= 0.0 || x < 0.0 || x > 1.0)
+		if (sigma <= 0.0)
 			SetIllegalArgument();
-		else if (x == 0.0 || x == 1.0)
-			SetNoValue();
 		else
-			PushDouble(gaussinv(x)*sigma + mue);
+        {
+            double y = normsinv(x);
+            if (nGlobalError)
+                PushInt(0);
+            else
+                PushDouble( y * sigma + mue);
+        }
 	}
 }
 
 void ScInterpreter::ScSNormInv()
 {
-	double x = GetDouble();
-	if (x < 0.0 || x > 1.0)
-		SetIllegalArgument();
-	else if (x == 0.0 || x == 1.0)
-		SetNoValue();
-	else
-		PushDouble(gaussinv(x));
+    PushDouble( normsinv( GetDouble()));
 }
 
 void ScInterpreter::ScLogNormInv()
 {
 	if ( MustHaveParamCount( GetByte(), 3 ) )
 	{
-		double sigma = GetDouble();					// Stdabw
-		double mue = GetDouble();					// Mittelwert
-		double y = GetDouble();						// y
-		if (sigma <= 0.0 || y <= 0.0 || y >= 1.0)
+		double sigma = GetDouble();
+		double mue = GetDouble();
+		double x = GetDouble();
+		if (sigma <= 0.0)
 			SetIllegalArgument();
 		else
-			PushDouble(exp(mue+sigma*gaussinv(y)));
+        {
+            double y = normsinv(x);
+            if (nGlobalError)
+                PushInt(0);
+            else
+                PushDouble( exp( y * sigma + mue));
+        }
 	}
 }
 
@@ -1417,7 +1488,7 @@
 		if (sigma <= 0.0 || alpha <= 0.0 || alpha >= 1.0 || n < 1.0)
 			SetIllegalArgument();
 		else
-			PushDouble( gaussinv(1.0-alpha/2.0) * sigma/sqrt(n) );
+			PushDouble( normsinv(1.0-alpha/2.0) * sigma/sqrt(n) );
 	}
 }