# HG changeset patch # User Arun Giridhar # Date 1656189636 14400 # Node ID bbb59cc6698c6c10db40dcd49bb3bf746ca75cf1 # Parent 7797481038fc0b99ab5f51a11d83c0e56d41ea3c factor.m: Performance tweak to avoid division in certain cases factor.m: For large numbers where only one factor lies outside the fast first round of divisions, check if it is prime before calling primes () trying to factorize it. This is up to 8000X faster for such numbers, and for an "average" input it gives a 22% to 40% speedup over a large number of trials. diff -r 7797481038fc -r bbb59cc6698c scripts/specfun/factor.m --- a/scripts/specfun/factor.m Tue Jun 21 15:43:27 2022 +0200 +++ b/scripts/specfun/factor.m Sat Jun 25 16:40:36 2022 -0400 @@ -150,23 +150,31 @@ ## q itself will be divided by those prime factors to become smaller, ## unless q was prime to begin with. - ## Now go all the way to sqrt(q), where q is smaller than the original q in - ## most cases. - ## - ## Note: Do not try to weed out the smallprimes inside largeprimes, whether - ## using length(smallprimes) or max(smallprimes) -- it slows it down! - largeprimes = primes (sqrt (q)); - [pf, q] = reducefactors (q, pf, largeprimes); + if (isprime (q)) + ## This is an optimization for numbers like 18446744073709551566 + ## == 2 * 9223372036854775783, where the small factors can be pulled + ## out easily and the remaining is prime. This optimization reduces + ## 14.3 s to 1.8 ms (8000X faster) for such cases. + pf(end+1) = q; + else + ## Now go all the way to sqrt(q), where q is smaller than the original q in + ## most cases. + ## + ## Note: Do not try to weed out the smallprimes inside largeprimes, whether + ## using length(smallprimes) or max(smallprimes) -- it slows it down! + largeprimes = primes (sqrt (q)); + [pf, q] = reducefactors (q, pf, largeprimes); - ## At this point, all prime factors <= the sqrt of the original q have been - ## pulled out in ascending order. - ## - ## If q = 1, then no further primes are left. - ## If q > 1, then q itself must be prime, and it must be the single prime - ## factor that was larger than the sqrt of the original q. - if (q > 1) - pf(end+1) = q; - endif + ## At this point, all prime factors <= the sqrt of the original q have been + ## pulled out in ascending order. + ## + ## If q = 1, then no further primes are left. + ## If q > 1, then q itself must be prime, and it must be the single prime + ## factor that was larger than the sqrt of the original q. + if (q > 1) + pf(end+1) = q; + endif + end ## At this point, all prime factors have been pulled out of q in ascending ## order. There is no need to sort(pf).