Fix #222: Wrong result for e.g. √(2.25) - 1.5 #222 (#2424)

* add snap zero in add/sub/log

* add a comment

* Update src/CalcManager/Ratpack/exp.cpp

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

* Update src/CalcManager/Ratpack/rat.cpp

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

* Update src/CalcManager/Ratpack/rat.cpp

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

---------

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

Author: guihuaz

src/CalcManager/Ratpack/exp.cpp

@@ -78,7 +78,7 @@ void exprat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
     const int32_t intpwr = rattoi32(pint, radix, precision);
     ratpowi32(&pwr, intpwr, precision);
 
-    subrat(px, pint, precision);
+    _subrat(px, pint, precision);
 
     // It just so happens to be an integral power of e.
     if (rat_gt(*px, rat_negsmallest, precision) && rat_lt(*px, rat_smallest, precision))
@@ -97,7 +97,7 @@ void exprat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
 
 //-----------------------------------------------------------------------------
 //
-//  FUNCTION: lograt, _lograt
+//  FUNCTION: lograt, _lograt, __lograt
 //
 //  ARGUMENTS: x PRAT representation of number to logarithim
 //
@@ -119,10 +119,13 @@ void exprat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
 //   Number is scaled between one and e_to_one_half prior to taking the
 //   log. This is to keep execution time from exploding.
 //
+//   lograt tries to snap to zero. Use _lograt inside ratpak by default.
+//   __lograt is part of _lograt private implementation and should not be used.
+//
 //
 //-----------------------------------------------------------------------------
 
-void _lograt(PRAT* px, int32_t precision)
+void __lograt(PRAT* px, int32_t precision)
 
 {
     CREATETAYLOR();
@@ -149,7 +152,7 @@ void _lograt(PRAT* px, int32_t precision)
     DESTROYTAYLOR();
 }
 
-void lograt(_Inout_ PRAT* px, int32_t precision)
+void _lograt(_Inout_ PRAT* px, int32_t precision)
 
 {
     PRAT pwr = nullptr;    // pwr is the large scaling factor.
@@ -192,18 +195,18 @@ void lograt(_Inout_ PRAT* px, int32_t precision)
     while (rat_gt(*px, e_to_one_half, precision))
     {
         divrat(px, e_to_one_half, precision);
-        addrat(&offset, rat_one, precision);
+        _addrat(&offset, rat_one, precision);
     }
 
-    _lograt(px, precision);
+    __lograt(px, precision);
 
     // Add the large and small scaling factors, take into account
     // small scaling was done in e_to_one_half chunks.
     divrat(&offset, rat_two, precision);
-    addrat(&pwr, offset, precision);
+    _addrat(&pwr, offset, precision);
 
     // And add the resulting scaling factor to the answer.
-    addrat(px, pwr, precision);
+    _addrat(px, pwr, precision);
 
     trimit(px, precision);
 
@@ -217,6 +220,18 @@ void lograt(_Inout_ PRAT* px, int32_t precision)
     destroyrat(pwr);
 }
 
+void lograt(_Inout_ PRAT* px, int32_t precision)
+
+{
+    PRAT a = nullptr;
+    DUPRAT(a, *px);
+
+    _lograt(px, precision);
+
+    _snaprat(px, a, nullptr, precision);
+    destroyrat(a);
+}
+
 void log10rat(_Inout_ PRAT* px, int32_t precision)
 
 {
@@ -235,8 +250,8 @@ bool IsEven(PRAT x, uint32_t radix, int32_t precision)
     DUPRAT(tmp, x);
     divrat(&tmp, rat_two, precision);
     fracrat(&tmp, radix, precision);
-    addrat(&tmp, tmp, precision);
-    subrat(&tmp, rat_one, precision);
+    _addrat(&tmp, tmp, precision);
+    _subrat(&tmp, rat_one, precision);
     if (rat_lt(tmp, rat_zero, precision))
     {
         bRet = true;
@@ -339,11 +354,11 @@ void powratNumeratorDenominator(_Inout_ PRAT* px, _In_ PRAT y, uint32_t radix, i
         DUPRAT(roundedResult, originalResult);
         if (roundedResult->pp->sign == -1)
         {
-            subrat(&roundedResult, rat_half, precision);
+            _subrat(&roundedResult, rat_half, precision);
         }
         else
         {
-            addrat(&roundedResult, rat_half, precision);
+            _addrat(&roundedResult, rat_half, precision);
         }
         intrat(&roundedResult, radix, precision);
 
@@ -423,7 +438,7 @@ void powratcomp(_Inout_ PRAT* px, _In_ PRAT y, uint32_t radix, int32_t precision
     {
         PRAT pxint = nullptr;
         DUPRAT(pxint, *px);
-        subrat(&pxint, rat_one, precision);
+        _subrat(&pxint, rat_one, precision);
         if (rat_gt(pxint, rat_negsmallest, precision) && rat_lt(pxint, rat_smallest, precision) && (sign == 1))
         {
             // *px is one, special case a 1 return.
@@ -442,12 +457,12 @@ void powratcomp(_Inout_ PRAT* px, _In_ PRAT y, uint32_t radix, int32_t precision
                 // If power is an integer let ratpowi32 deal with it.
                 PRAT iy = nullptr;
                 DUPRAT(iy, y);
-                subrat(&iy, podd, precision);
+                _subrat(&iy, podd, precision);
                 int32_t inty = rattoi32(iy, radix, precision);
 
                 PRAT plnx = nullptr;
                 DUPRAT(plnx, *px);
-                lograt(&plnx, precision);
+                _lograt(&plnx, precision);
                 mulrat(&plnx, iy, precision);
                 if (rat_gt(plnx, rat_max_exp, precision) || rat_lt(plnx, rat_min_exp, precision))
                 {
@@ -515,7 +530,7 @@ void powratcomp(_Inout_ PRAT* px, _In_ PRAT y, uint32_t radix, int32_t precision
                     sign = 1;
                 }
 
-                lograt(px, precision);
+                _lograt(px, precision);
                 mulrat(px, y, precision);
                 exprat(px, radix, precision);
             }

src/CalcManager/Ratpack/fact.cpp

@@ -15,7 +15,6 @@
 //-----------------------------------------------------------------------------
 #include "ratpak.h"
 
-#define ABSRAT(x) (((x)->pp->sign = 1), ((x)->pq->sign = 1))
 #define NEGATE(x) ((x)->pp->sign *= -1)
 
 //-----------------------------------------------------------------------------
@@ -73,35 +72,35 @@ void _gamma(PRAT* pn, uint32_t radix, int32_t precision)
 
     // Find the best 'A' for convergence to the required precision.
     a = i32torat(radix);
-    lograt(&a, precision);
+    _lograt(&a, precision);
     mulrat(&a, ratprec, precision);
 
     // Really is -ln(n)+1, but -ln(n) will be < 1
     // if we scale n between 0.5 and 1.5
-    addrat(&a, rat_two, precision);
+    _addrat(&a, rat_two, precision);
     DUPRAT(tmp, a);
-    lograt(&tmp, precision);
+    _lograt(&tmp, precision);
     mulrat(&tmp, *pn, precision);
-    addrat(&a, tmp, precision);
-    addrat(&a, rat_one, precision);
+    _addrat(&a, tmp, precision);
+    _addrat(&a, rat_one, precision);
 
     // Calculate the necessary bump in precision and up the precision.
     // The following code is equivalent to
     // precision += ln(exp(a)*pow(a,n+1.5))-ln(radix));
     DUPRAT(tmp, *pn);
     one_pt_five = i32torat(3L);
     divrat(&one_pt_five, rat_two, precision);
-    addrat(&tmp, one_pt_five, precision);
+    _addrat(&tmp, one_pt_five, precision);
     DUPRAT(term, a);
     powratcomp(&term, tmp, radix, precision);
     DUPRAT(tmp, a);
     exprat(&tmp, radix, precision);
     mulrat(&term, tmp, precision);
-    lograt(&term, precision);
+    _lograt(&term, precision);
     const auto ratRadix = i32torat(radix);
     DUPRAT(tmp, ratRadix);
-    lograt(&tmp, precision);
-    subrat(&term, tmp, precision);
+    _lograt(&tmp, precision);
+    _subrat(&term, tmp, precision);
     precision += rattoi32(term, radix, precision);
 
     // Set up initial terms for series, refer to series in above comment block.
@@ -118,10 +117,10 @@ void _gamma(PRAT* pn, uint32_t radix, int32_t precision)
     DUPRAT(sum, rat_one);
     divrat(&sum, *pn, precision);
     DUPRAT(tmp, *pn);
-    addrat(&tmp, rat_one, precision);
+    _addrat(&tmp, rat_one, precision);
     DUPRAT(term, a);
     divrat(&term, tmp, precision);
-    subrat(&sum, term, precision);
+    _subrat(&sum, term, precision);
 
     DUPRAT(err, ratRadix);
     NEGATE(ratprec);
@@ -134,7 +133,7 @@ void _gamma(PRAT* pn, uint32_t radix, int32_t precision)
     // Loop until precision is reached, or asked to halt.
     while (!zerrat(term) && rat_gt(term, err, precision))
     {
-        addrat(pn, rat_two, precision);
+        _addrat(pn, rat_two, precision);
 
         // WARNING: mixing numbers and  rationals here.
         // for speed and efficiency.
@@ -146,22 +145,22 @@ void _gamma(PRAT* pn, uint32_t radix, int32_t precision)
         divrat(&factorial, a2, precision);
 
         DUPRAT(tmp, *pn);
-        addrat(&tmp, rat_one, precision);
+        _addrat(&tmp, rat_one, precision);
         destroyrat(term);
         createrat(term);
         DUPNUM(term->pp, count);
         DUPNUM(term->pq, num_one);
-        addrat(&term, rat_one, precision);
+        _addrat(&term, rat_one, precision);
         mulrat(&term, tmp, precision);
         DUPRAT(tmp, a);
         divrat(&tmp, term, precision);
 
         DUPRAT(term, rat_one);
         divrat(&term, *pn, precision);
-        subrat(&term, tmp, precision);
+        _subrat(&term, tmp, precision);
 
         divrat(&term, factorial, precision);
-        addrat(&sum, term, precision);
+        _addrat(&sum, term, precision);
         ABSRAT(term);
     }
 
@@ -214,7 +213,7 @@ void factrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
     while (rat_gt(*px, rat_zero, precision) && (LOGRATRADIX(*px) > -precision))
     {
         mulrat(&fact, *px, precision);
-        subrat(px, rat_one, precision);
+        _subrat(px, rat_one, precision);
     }
 
     // Added to make numbers 'close enough' to integers use integer factorial.
@@ -226,13 +225,13 @@ void factrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
 
     while (rat_lt(*px, neg_rat_one, precision))
     {
-        addrat(px, rat_one, precision);
+        _addrat(px, rat_one, precision);
         divrat(&fact, *px, precision);
     }
 
     if (rat_neq(*px, rat_zero, precision))
     {
-        addrat(px, rat_one, precision);
+        _addrat(px, rat_one, precision);
         _gamma(px, radix, precision);
         mulrat(px, fact, precision);
     }

src/CalcManager/Ratpack/itrans.cpp

@@ -95,7 +95,7 @@ void asinrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
 
     // Avoid the really bad part of the asin curve near +/-1.
     DUPRAT(phack, *px);
-    subrat(&phack, rat_one, precision);
+    _subrat(&phack, rat_one, precision);
     // Since *px might be epsilon near zero we must set it to zero.
     if (rat_le(phack, rat_smallest, precision) && rat_ge(phack, rat_negsmallest, precision))
     {
@@ -109,7 +109,7 @@ void asinrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
         {
             if (rat_gt(*px, rat_one, precision))
             {
-                subrat(px, rat_one, precision);
+                _subrat(px, rat_one, precision);
                 if (rat_gt(*px, rat_smallest, precision))
                 {
                     throw(CALC_E_DOMAIN);
@@ -122,11 +122,11 @@ void asinrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
             DUPRAT(pret, *px);
             mulrat(px, pret, precision);
             (*px)->pp->sign *= -1;
-            addrat(px, rat_one, precision);
+            _addrat(px, rat_one, precision);
             rootrat(px, rat_two, radix, precision);
             _asinrat(px, precision);
             (*px)->pp->sign *= -1;
-            addrat(px, pi_over_two, precision);
+            _addrat(px, pi_over_two, precision);
             destroyrat(pret);
         }
         else
@@ -214,7 +214,7 @@ void acosrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
         (*px)->pp->sign = sgn;
         asinrat(px, radix, precision);
         (*px)->pp->sign *= -1;
-        addrat(px, pi_over_two, precision);
+        _addrat(px, pi_over_two, precision);
     }
 }
 
@@ -297,15 +297,15 @@ void atanrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
             tmpx->pp->sign = sgn;
             tmpx->pq->sign = 1;
             DUPRAT(*px, pi_over_two);
-            subrat(px, tmpx, precision);
+            _subrat(px, tmpx, precision);
             destroyrat(tmpx);
         }
         else
         {
             (*px)->pp->sign = sgn;
             DUPRAT(tmpx, *px);
             mulrat(&tmpx, *px, precision);
-            addrat(&tmpx, rat_one, precision);
+            _addrat(&tmpx, rat_one, precision);
             rootrat(&tmpx, rat_two, radix, precision);
             divrat(px, tmpx, precision);
             destroyrat(tmpx);
@@ -322,6 +322,6 @@ void atanrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
     }
     if (rat_gt(*px, pi_over_two, precision))
     {
-        subrat(px, pi, precision);
+        _subrat(px, pi, precision);
     }
 }

src/CalcManager/Ratpack/itransh.cpp

@@ -59,10 +59,10 @@ void asinhrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
         PRAT ptmp = nullptr;
         DUPRAT(ptmp, (*px));
         mulrat(&ptmp, *px, precision);
-        addrat(&ptmp, rat_one, precision);
+        _addrat(&ptmp, rat_one, precision);
         rootrat(&ptmp, rat_two, radix, precision);
-        addrat(px, ptmp, precision);
-        lograt(px, precision);
+        _addrat(px, ptmp, precision);
+        _lograt(px, precision);
         destroyrat(ptmp);
     }
     else
@@ -113,10 +113,10 @@ void acoshrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
         PRAT ptmp = nullptr;
         DUPRAT(ptmp, (*px));
         mulrat(&ptmp, *px, precision);
-        subrat(&ptmp, rat_one, precision);
+        _subrat(&ptmp, rat_one, precision);
         rootrat(&ptmp, rat_two, radix, precision);
-        addrat(px, ptmp, precision);
-        lograt(px, precision);
+        _addrat(px, ptmp, precision);
+        _lograt(px, precision);
         destroyrat(ptmp);
     }
 }
@@ -143,11 +143,11 @@ void atanhrat(_Inout_ PRAT* px, int32_t precision)
 {
     PRAT ptmp = nullptr;
     DUPRAT(ptmp, (*px));
-    subrat(&ptmp, rat_one, precision);
-    addrat(px, rat_one, precision);
+    _subrat(&ptmp, rat_one, precision);
+    _addrat(px, rat_one, precision);
     divrat(px, ptmp, precision);
     (*px)->pp->sign *= -1;
-    lograt(px, precision);
+    _lograt(px, precision);
     divrat(px, rat_two, precision);
     destroyrat(ptmp);
 }

src/CalcManager/Ratpack/logic.cpp

@@ -246,7 +246,7 @@ void modrat(_Inout_ PRAT* pa, _In_ PRAT b)
 
     if (needAdjust && !zerrat(*pa))
     {
-        addrat(pa, b, BASEX);
+        _addrat(pa, b, BASEX);
     }
 
     // Get *pa back in the integer over integer form.