Mercurial > octave

changeset 31402:0783418ac443

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
factor.m: Guard against non-prime factors returned from C++. factor.m: Extensive testing reveals a small chance that __pollardrho__ can return a composite number as a factor for some inputs. Call isprime to trap that case. Add BISTs.
author Arun Giridhar <arungiridhar@gmail.com>
date Sun, 06 Nov 2022 11:40:07 -0500
parents e7fc6251b698
children 821cdacfdca7
files scripts/specfun/factor.m
diffstat 1 files changed, 29 insertions(+), 4 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/specfun/factor.m	Sat Nov 05 19:15:13 2022 +0100
+++ b/scripts/specfun/factor.m	Sun Nov 06 11:40:07 2022 -0500
@@ -95,18 +95,27 @@
   sortflag = false;
   if (isprime (q))
     pf(end+1) = q;
-  else
+  elseif (q > 1)
     ## Use Pollard Rho technique to pull factors one at a time.
     while (q > 1e10 && ! isprime (q))
-      pr = feval (cls, __pollardrho__ (q));  # pr is a prime factor.
+      pr = feval (cls, __pollardrho__ (q));  # pr is a factor of q.
+
+      ## There is a small chance (13 in 1e5) that pr is not actually prime.
+      ## To guard against that, factorize pr, which will force smaller factors
+      ## to be found.  The use of isprime above guards against infinite
+      ## recursion.
+      if (! isprime (pr))
+        pr = factor (pr);
+      end
+
       [pf, q] = reducefactors (q, pf, pr);
-      ## q is now divided by all occurrences of factor pr.
+      ## q is now divided by all occurrences of factor(s) pr.
       sortflag = true;
     endwhile
 
     if (isprime (q))
       pf(end+1) = q;
-    else
+    elseif (q > 1)
       ## If we are here, then q is composite but less than 1e10,
       ## and that is fast enough to test by division.
       largeprimes = primes (feval (cls, sqrt (q)));
@@ -167,6 +176,22 @@
 %!   assert (all ([0,pf] != [pf,0]));
 %! endfor
 
+## Make sure that all factors returned are indeed prime, even when
+## __pollardrho__ returns a composite factor.
+%!assert (all (isprime (factor (uint64 (18446744073707633197)))))
+%!assert (all (isprime (factor (uint64 (18446744073707551733)))))
+%!assert (all (isprime (factor (uint64 (18446744073709427857)))))
+%!assert (all (isprime (factor (uint64 (18446744073709396891)))))
+%!assert (all (isprime (factor (uint64 (18446744073708666563)))))
+%!assert (all (isprime (factor (uint64 (18446744073708532009)))))
+%!assert (all (isprime (factor (uint64 (18446744073708054211)))))
+%!assert (all (isprime (factor (uint64 (18446744073707834741)))))
+%!assert (all (isprime (factor (uint64 (18446744073707298053)))))
+%!assert (all (isprime (factor (uint64 (18446744073709407383)))))
+%!assert (all (isprime (factor (uint64 (18446744073708730121)))))
+%!assert (all (isprime (factor (uint64 (18446744073708104447)))))
+%!assert (all (isprime (factor (uint64 (18446744073709011493)))))
+
 %!assert (factor (uint8 (8)), uint8 ([2 2 2]))
 %!assert (factor (single (8)), single ([2 2 2]))
 %!test