Mercurial > jwe > octave

comparison scripts/specfun/factor.m @ 31112:bbb59cc6698c stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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.
author Arun Giridhar <arungiridhar@gmail.com>
date Sat, 25 Jun 2022 16:40:36 -0400
parents 796f54d4ddbf
children
comparison
equal deleted inserted replaced
31108:7797481038fc 31112:bbb59cc6698c
148 ## repetitions, in ascending order. 148 ## repetitions, in ascending order.
149 ## 149 ##
150 ## q itself will be divided by those prime factors to become smaller, 150 ## q itself will be divided by those prime factors to become smaller,
151 ## unless q was prime to begin with. 151 ## unless q was prime to begin with.
152 152
153 ## Now go all the way to sqrt(q), where q is smaller than the original q in 153 if (isprime (q))
154 ## most cases. 154 ## This is an optimization for numbers like 18446744073709551566
155 ## 155 ## == 2 * 9223372036854775783, where the small factors can be pulled
156 ## Note: Do not try to weed out the smallprimes inside largeprimes, whether 156 ## out easily and the remaining is prime. This optimization reduces
157 ## using length(smallprimes) or max(smallprimes) -- it slows it down! 157 ## 14.3 s to 1.8 ms (8000X faster) for such cases.
158 largeprimes = primes (sqrt (q));
159 [pf, q] = reducefactors (q, pf, largeprimes);
160
161 ## At this point, all prime factors <= the sqrt of the original q have been
162 ## pulled out in ascending order.
163 ##
164 ## If q = 1, then no further primes are left.
165 ## If q > 1, then q itself must be prime, and it must be the single prime
166 ## factor that was larger than the sqrt of the original q.
167 if (q > 1)
168 pf(end+1) = q; 158 pf(end+1) = q;
169 endif 159 else
160 ## Now go all the way to sqrt(q), where q is smaller than the original q in
161 ## most cases.
162 ##
163 ## Note: Do not try to weed out the smallprimes inside largeprimes, whether
164 ## using length(smallprimes) or max(smallprimes) -- it slows it down!
165 largeprimes = primes (sqrt (q));
166 [pf, q] = reducefactors (q, pf, largeprimes);
167
168 ## At this point, all prime factors <= the sqrt of the original q have been
169 ## pulled out in ascending order.
170 ##
171 ## If q = 1, then no further primes are left.
172 ## If q > 1, then q itself must be prime, and it must be the single prime
173 ## factor that was larger than the sqrt of the original q.
174 if (q > 1)
175 pf(end+1) = q;
176 endif
177 end
170 178
171 ## At this point, all prime factors have been pulled out of q in ascending 179 ## At this point, all prime factors have been pulled out of q in ascending
172 ## order. There is no need to sort(pf). 180 ## order. There is no need to sort(pf).
173 181
174 ## Determine multiplicity. 182 ## Determine multiplicity.