Changeset 42846 for trunk/psLib

Timestamp:

May 9, 2025, 12:38:05 PM (14 months ago)

Author:

eugene

Message:

fixing an error with psBinomialCoeff due to integerizing a fraction; this function is only used by the 1D vector fitting function (regular and clipped) and only manifests for polynomials of order > 2; I have scanned the code and the only place this happens is in the overscan correction for which the order is set in the recipe; some cameras other than gpc1 and gpc2 have higher order overscan models and would have been affected

File:

• trunk/psLib/src/math/psBinomialCoeff.c (modified) (1 diff)

trunk/psLib/src/math/psBinomialCoeff.c

@@ -22 +22 @@
   if (k == 0) { return (b); }
 
+  psS64 num = 1;
+  psS64 den = 1;
+
   for (i = 1; i <= k; i++) {
-    b *= (N - (k - i))/i;
+
+      // old version: b *= (N - (k - i))/i;
+      // wikipedia recursion version: b *= (N + 1 - i)/i;
+
+      // the old version just counts the elements in the opposite order from the wikipedia version:
+
+      // numerators would be, e.g., for C(6,4):
+      // old  version (6 - (4 - 1)), (6 - (4 - 2)), (6 - (4 - 3)), (6 - (4 - 4)) = 3, 4, 5, 6
+      // wiki version (6 + 1 - 1), (6 + 1 - 2), (6 + 1 - 3), (6 + 1 - 4) = 6, 5, 4, 3
+
+      // denominators would be the same: 1, 2, 3, 4.
+
+      // since these are multiplied, the order of the numerators is independent of the order of the denominators
+
+      // but both of these formulae are bad because the factors can be fractional, but that means we cannot use an
+      // integer to carry the intermediate result.  
+
+      // instead, we need to calculate the integer product of the denominators and the integer
+      // product of the numerators and divide at the end
+      
+      // this causes an error for C(3,2), C(5,2), C(4,3), etc.
+
+      num *= (N + 1 - i);
+      den *= i;
   }
+
+  b = num / den;
 
   return(b);