octave: scripts/specfun/factor.m annotate

annotate scripts/specfun/factor.m @ 31107:bbb59cc6698c stable

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	068342cc93b8

rev	line source
27923 bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	1 ########################################################################
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	2 ##
30564 796f54d4ddbf update Octave Project Developers copyright for the new year John W. Eaton <jwe@octave.org> parents: 30379 diff changeset	3 ## Copyright (C) 2000-2022 The Octave Project Developers
27918 b442ec6dda5c use centralized file for copyright info for individual contributors John W. Eaton <jwe@octave.org> parents: 26376 diff changeset	4 ##
27923 bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	5 ## See the file COPYRIGHT.md in the top-level directory of this
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	6 ## distribution or <https://octave.org/copyright/>.
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	7 ##
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	8 ## This file is part of Octave.
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	9 ##
24534 194eb4bd202b maint: Update punctuation for GPL v3 license text. Rik <rik@octave.org> parents: 23220 diff changeset	10 ## Octave is free software: you can redistribute it and/or modify it
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	11 ## under the terms of the GNU General Public License as published by
24534 194eb4bd202b maint: Update punctuation for GPL v3 license text. Rik <rik@octave.org> parents: 23220 diff changeset	12 ## the Free Software Foundation, either version 3 of the License, or
22755 3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22323 diff changeset	13 ## (at your option) any later version.
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	14 ##
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	15 ## Octave is distributed in the hope that it will be useful, but
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	16 ## WITHOUT ANY WARRANTY; without even the implied warranty of
22755 3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22323 diff changeset	17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22323 diff changeset	18 ## GNU General Public License for more details.
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	19 ##
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	20 ## You should have received a copy of the GNU General Public License
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 7001 diff changeset	21 ## along with Octave; see the file COPYING. If not, see
24534 194eb4bd202b maint: Update punctuation for GPL v3 license text. Rik <rik@octave.org> parents: 23220 diff changeset	22 ## <https://www.gnu.org/licenses/>.
27923 bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	23 ##
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	24 ########################################################################
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	25
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	26 ## -- texinfo --
20852 516bb87ea72e 2015 Code Sprint: remove class of function from docstring for all m-files. Rik <rik@octave.org> parents: 20486 diff changeset	27 ## @deftypefn {} {@var{pf} =} factor (@var{q})
516bb87ea72e 2015 Code Sprint: remove class of function from docstring for all m-files. Rik <rik@octave.org> parents: 20486 diff changeset	28 ## @deftypefnx {} {[@var{pf}, @var{n}] =} factor (@var{q})
19097 e7bffcea619f factor.m: Overhaul function. Rik <rik@octave.org> parents: 17744 diff changeset	29 ## Return the prime factorization of @var{q}.
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	30 ##
19097 e7bffcea619f factor.m: Overhaul function. Rik <rik@octave.org> parents: 17744 diff changeset	31 ## The prime factorization is defined as @code{prod (@var{pf}) == @var{q}}
e7bffcea619f factor.m: Overhaul function. Rik <rik@octave.org> parents: 17744 diff changeset	32 ## where every element of @var{pf} is a prime number. If @code{@var{q} == 1},
25270 617fe022e965 factor.m: Return output factors of same class as input for Matlab compatibility. Rik <rik@octave.org> parents: 25269 diff changeset	33 ## return 1. The output @var{pf} is of the same numeric class as the input.
19593 446c46af4b42 strip trailing whitespace from most source files John W. Eaton <jwe@octave.org> parents: 17744 diff changeset	34 ##
19097 e7bffcea619f factor.m: Overhaul function. Rik <rik@octave.org> parents: 17744 diff changeset	35 ## With two output arguments, return the unique prime factors @var{pf} and
21546 f7f97d7e9294 doc: Wrap m-file docstrings to 79 characters + newline (80 total). Rik <rik@octave.org> parents: 20852 diff changeset	36 ## their multiplicities. That is,
f7f97d7e9294 doc: Wrap m-file docstrings to 79 characters + newline (80 total). Rik <rik@octave.org> parents: 20852 diff changeset	37 ## @code{prod (@var{pf} .^ @var{n}) == @var{q}}.
19596 0e1f5a750d00 maint: Periodic merge of gui-release to default. John W. Eaton <jwe@octave.org> parents: 19152 19593 diff changeset	38 ##
25270 617fe022e965 factor.m: Return output factors of same class as input for Matlab compatibility. Rik <rik@octave.org> parents: 25269 diff changeset	39 ## Implementation Note: The input @var{q} must be less than @code{flintmax}
30228 39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	40 ## when the input is a floating-point class (double or single).
19097 e7bffcea619f factor.m: Overhaul function. Rik <rik@octave.org> parents: 17744 diff changeset	41 ## @seealso{gcd, lcm, isprime, primes}
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	42 ## @end deftypefn
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	43
19097 e7bffcea619f factor.m: Overhaul function. Rik <rik@octave.org> parents: 17744 diff changeset	44 function [pf, n] = factor (q)
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	45
28891 de5f2f9a64ff maint: Use same coding style when checking for a minimum of 1 input. Rik <rik@octave.org> parents: 27978 diff changeset	46 if (nargin < 1)
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	47 print_usage ();
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	48 endif
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	49
25269 cac96fd5310d factor.m: Emit an error if input is negative (bug #53425). Dildar Sk <dildarsk101010@gmail.com> parents: 25054 diff changeset	50 if (! isscalar (q) \|\| ! isreal (q) \|\| q < 0 \|\| q != fix (q))
cac96fd5310d factor.m: Emit an error if input is negative (bug #53425). Dildar Sk <dildarsk101010@gmail.com> parents: 25054 diff changeset	51 error ("factor: Q must be a real non-negative integer");
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	52 endif
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	53
30330 01de0045b2e3 maint: Shorten some long lines to <= 80 characters (bug #57599) Rik <rik@octave.org> parents: 30261 diff changeset	54 ## Special case if q is prime, because isprime() is now much faster than
01de0045b2e3 maint: Shorten some long lines to <= 80 characters (bug #57599) Rik <rik@octave.org> parents: 30261 diff changeset	55 ## factor(). This also absorbs the case of q < 4, where there are no primes
01de0045b2e3 maint: Shorten some long lines to <= 80 characters (bug #57599) Rik <rik@octave.org> parents: 30261 diff changeset	56 ## less than sqrt(q).
30261 a49c635b179d isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312). Arun Giridhar <arungiridhar@gmail.com> parents: 30228 diff changeset	57 if (q < 4 \|\| isprime (q))
19097 e7bffcea619f factor.m: Overhaul function. Rik <rik@octave.org> parents: 17744 diff changeset	58 pf = q;
11469 c776f063fefe Overhaul m-script files to use common variable name between code and documentation. Rik <octave@nomad.inbox5.com> parents: 10476 diff changeset	59 n = 1;
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	60 return;
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	61 endif
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	62
30261 a49c635b179d isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312). Arun Giridhar <arungiridhar@gmail.com> parents: 30228 diff changeset	63 ## If we are here, then q is composite.
a49c635b179d isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312). Arun Giridhar <arungiridhar@gmail.com> parents: 30228 diff changeset	64
29983 ecbcc4647dbe factor.m: Overhaul function to support inputs > flintmax. Rik <rik@octave.org> parents: 29359 diff changeset	65 cls = class (q); # store class
ecbcc4647dbe factor.m: Overhaul function to support inputs > flintmax. Rik <rik@octave.org> parents: 29359 diff changeset	66 if (isfloat (q) && q > flintmax (q))
29984 c633d34960b4 factor.m: Fix typo in error() message. Rik <rik@octave.org> parents: 29983 diff changeset	67 error ("factor: Q too large to factor (> flintmax)");
29983 ecbcc4647dbe factor.m: Overhaul function to support inputs > flintmax. Rik <rik@octave.org> parents: 29359 diff changeset	68 endif
ecbcc4647dbe factor.m: Overhaul function to support inputs > flintmax. Rik <rik@octave.org> parents: 29359 diff changeset	69
30228 39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	70 ## The basic idea is to divide by the prime numbers from 1 to sqrt(q).
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	71 ## But primes(sqrt(q)) can be very time-consuming to compute for q > 1e16,
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	72 ## so we divide by smaller primes first.
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	73 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	74 ## This won't make a difference for prime q, but it makes a big (100x)
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	75 ## difference for large composite q. Since there are many more composites
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	76 ## than primes, this leads overall to a speedup.
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	77 ##
11469 c776f063fefe Overhaul m-script files to use common variable name between code and documentation. Rik <octave@nomad.inbox5.com> parents: 10476 diff changeset	78 ## There is at most one prime greater than sqrt(q), and if it exists,
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	79 ## it has multiplicity 1, so no need to consider any factors greater
30228 39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	80 ## than sqrt(q) directly. If there were two factors p1, p2 > sqrt(q), then
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	81 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	82 ## q >= p1p2 > sqrt(q)sqrt(q) == q,
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	83 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	84 ## which is a contradiction.
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	85 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	86 ## The following calculation of transition and number of divisors to use
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	87 ## was determined empirically. As of now (October 2021) it gives the best
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	88 ## overall performance over the range of 1 <= q <= intmax ("uint64").
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	89 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	90 ## For future programmers: check periodically for performance improvements
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	91 ## and tune this transition as required. Trials that didn't yield success
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	92 ## in (October 2021):
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	93 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	94 ## 1.) persistent smallprimes = primes (FOO)
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	95 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	96 ## For various fixed FOO in the range 10 <= FOO <= 10e6.
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	97 ## (FOO is independent of q.) The thought had been that making it
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	98 ## persistent would cache it so it didn't need to be recomputed for
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	99 ## subsequent calls, but it slowed it down overall. It seems calling
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	100 ## primes twice with smaller q is still faster than one persistent
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	101 ## call for a large q.
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	102 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	103 ## 2.) smallprimes = primes (q ^ FOO)
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	104 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	105 ## For various values of FOO. For FOO >= 0.25 or FOO <= 0.16, the
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	106 ## performance is very poor. FOO needs to be in the 0.17 to 0.24 range,
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	107 ## somewhat. Benchmark experiments indicate it should increase gently
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	108 ## from 0.18 to 0.21 as q goes from 10^11 to 10^18.
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	109 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	110 ## But putting in such an expression would require calculating the log
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	111 ## of q, which defeats any performance improvement. Or a step-wise
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	112 ## approximation like:
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	113 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	114 ## foo = 0.18 + 0.01 * (q > 1e12) + 0.01 * (q > 1e14) ...
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	115 ## + 0.01 * (q > 1e16);
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	116 ## smallprimes = primes (feval (cls, q^foo));
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	117 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	118 ## where the RHS of foo would go from 0.18 to 0.21 over several orders
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	119 ## of magnitude without calling the log. Obviously that is overly
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	120 ## empirical, so putting in q^0.2 seems to be the most robust overall
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	121 ## for 64-bit q.
29983 ecbcc4647dbe factor.m: Overhaul function to support inputs > flintmax. Rik <rik@octave.org> parents: 29359 diff changeset	122
30228 39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	123 ## Lookup table for sufficiently small values for q.
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	124 if (q < 10e9)
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	125 ## Lookup, rather calling up to primes(100) is about 3% faster, than the
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	126 ## previous value of primes(30). Same for very small q < 1e6.
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	127 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	128 ## For 1e9 < q < 10e9 the lookup approach is about 7% faster.
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	129
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	130 smallprimes = feval (cls, ...
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	131 [2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97]);
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	132
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	133 ## Only for really small values of q, statements like
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	134 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	135 ## smallprimes(smallprimes > q) = [];
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	136 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	137 ## are relevant and slow down significantly for large values of q.
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	138 else
30379 363fb10055df maint: Style check m-files ahead of 7.1 release. Rik <rik@octave.org> parents: 30330 diff changeset	139 ## For sufficiently large q, go up to the 5th root of q for now.
30228 39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	140 smallprimes = primes (feval (cls, q^0.2));
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	141 endif
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	142
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	143 ## pf is the list of prime factors returned with type of input class.
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	144 pf = feval (cls, []);
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	145 [pf, q] = reducefactors (q, pf, smallprimes);
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	146
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	147 ## pf now contains all prime factors of q within smallprimes, including
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	148 ## repetitions, in ascending order.
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	149 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	150 ## q itself will be divided by those prime factors to become smaller,
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	151 ## unless q was prime to begin with.
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	152
31107 bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	153 if (isprime (q))
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	154 ## This is an optimization for numbers like 18446744073709551566
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	155 ## == 2 * 9223372036854775783, where the small factors can be pulled
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	156 ## out easily and the remaining is prime. This optimization reduces
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	157 ## 14.3 s to 1.8 ms (8000X faster) for such cases.
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	158 pf(end+1) = q;
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	159 else
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	160 ## Now go all the way to sqrt(q), where q is smaller than the original q in
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	161 ## most cases.
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	162 ##
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	163 ## Note: Do not try to weed out the smallprimes inside largeprimes, whether
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	164 ## using length(smallprimes) or max(smallprimes) -- it slows it down!
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	165 largeprimes = primes (sqrt (q));
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	166 [pf, q] = reducefactors (q, pf, largeprimes);
30228 39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	167
31107 bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	168 ## At this point, all prime factors <= the sqrt of the original q have been
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	169 ## pulled out in ascending order.
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	170 ##
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	171 ## If q = 1, then no further primes are left.
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	172 ## If q > 1, then q itself must be prime, and it must be the single prime
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	173 ## factor that was larger than the sqrt of the original q.
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	174 if (q > 1)
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	175 pf(end+1) = q;
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	176 endif
bbb59cc6698c factor.m: Performance tweak to avoid division in certain cases Arun Giridhar <arungiridhar@gmail.com> parents: 30564 diff changeset	177 end
30228 39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	178
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	179 ## At this point, all prime factors have been pulled out of q in ascending
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	180 ## order. There is no need to sort(pf).
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	181
25269 cac96fd5310d factor.m: Emit an error if input is negative (bug #53425). Dildar Sk <dildarsk101010@gmail.com> parents: 25054 diff changeset	182 ## Determine multiplicity.
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	183 if (nargout > 1)
19097 e7bffcea619f factor.m: Overhaul function. Rik <rik@octave.org> parents: 17744 diff changeset	184 idx = find ([0, pf] != [pf, 0]);
e7bffcea619f factor.m: Overhaul function. Rik <rik@octave.org> parents: 17744 diff changeset	185 pf = pf(idx(1:length (idx)-1));
11469 c776f063fefe Overhaul m-script files to use common variable name between code and documentation. Rik <octave@nomad.inbox5.com> parents: 10476 diff changeset	186 n = diff (idx);
5827 1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	187 endif
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	188
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	189 endfunction
1fe78adb91bc [project @ 2006-05-22 06:25:14 by jwe] jwe parents: diff changeset	190
30228 39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	191 function [pf, q] = reducefactors (qin, pfin, divisors)
30379 363fb10055df maint: Style check m-files ahead of 7.1 release. Rik <rik@octave.org> parents: 30330 diff changeset	192
30228 39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	193 pf = pfin;
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	194 q = qin;
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	195 ## The following line is a few milliseconds faster than
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	196 ## divisors (mod (q, divisors) ~= 0) = [];
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	197 divisors = divisors (mod (q, divisors) == 0);
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	198
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	199 for pp = divisors # for each factor in turn
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	200 ## Keep extracting all occurrences of that factor before going to larger
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	201 ## factors.
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	202 ##
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	203 ## Note: mod() was marginally faster than rem(), when assessed over 10e6
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	204 ## trials of the whole factor() function.
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	205 while (mod (q, pp) == 0)
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	206 pf(end+1) = pp;
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	207 q /= pp;
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	208 endwhile
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	209 endfor
30379 363fb10055df maint: Style check m-files ahead of 7.1 release. Rik <rik@octave.org> parents: 30330 diff changeset	210
30228 39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	211 endfunction
39a4ab124fd0 factor.m: significant speedup for large input quantities (> 1e14) (bug #61129). Arun Giridhar <arungiridhar@gmail.com> parents: 29984 diff changeset	212
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	213
29983 ecbcc4647dbe factor.m: Overhaul function to support inputs > flintmax. Rik <rik@octave.org> parents: 29359 diff changeset	214 ## Test special case input
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	215 %!assert (factor (1), 1)
29983 ecbcc4647dbe factor.m: Overhaul function to support inputs > flintmax. Rik <rik@octave.org> parents: 29359 diff changeset	216 %!assert (factor (2), 2)
ecbcc4647dbe factor.m: Overhaul function to support inputs > flintmax. Rik <rik@octave.org> parents: 29359 diff changeset	217 %!assert (factor (3), 3)
ecbcc4647dbe factor.m: Overhaul function to support inputs > flintmax. Rik <rik@octave.org> parents: 29359 diff changeset	218
11469 c776f063fefe Overhaul m-script files to use common variable name between code and documentation. Rik <octave@nomad.inbox5.com> parents: 10476 diff changeset	219 %!test
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	220 %! for i = 2:20
19097 e7bffcea619f factor.m: Overhaul function. Rik <rik@octave.org> parents: 17744 diff changeset	221 %! pf = factor (i);
e7bffcea619f factor.m: Overhaul function. Rik <rik@octave.org> parents: 17744 diff changeset	222 %! assert (prod (pf), i);
e7bffcea619f factor.m: Overhaul function. Rik <rik@octave.org> parents: 17744 diff changeset	223 %! assert (all (isprime (pf)));
e7bffcea619f factor.m: Overhaul function. Rik <rik@octave.org> parents: 17744 diff changeset	224 %! [pf, n] = factor (i);
e7bffcea619f factor.m: Overhaul function. Rik <rik@octave.org> parents: 17744 diff changeset	225 %! assert (prod (pf.^n), i);
e7bffcea619f factor.m: Overhaul function. Rik <rik@octave.org> parents: 17744 diff changeset	226 %! assert (all ([0,pf] != [pf,0]));
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	227 %! endfor
11469 c776f063fefe Overhaul m-script files to use common variable name between code and documentation. Rik <octave@nomad.inbox5.com> parents: 10476 diff changeset	228
25270 617fe022e965 factor.m: Return output factors of same class as input for Matlab compatibility. Rik <rik@octave.org> parents: 25269 diff changeset	229 %!assert (factor (uint8 (8)), uint8 ([2 2 2]))
617fe022e965 factor.m: Return output factors of same class as input for Matlab compatibility. Rik <rik@octave.org> parents: 25269 diff changeset	230 %!assert (factor (single (8)), single ([2 2 2]))
617fe022e965 factor.m: Return output factors of same class as input for Matlab compatibility. Rik <rik@octave.org> parents: 25269 diff changeset	231 %!test
617fe022e965 factor.m: Return output factors of same class as input for Matlab compatibility. Rik <rik@octave.org> parents: 25269 diff changeset	232 %! [pf, n] = factor (int16 (8));
617fe022e965 factor.m: Return output factors of same class as input for Matlab compatibility. Rik <rik@octave.org> parents: 25269 diff changeset	233 %! assert (pf, int16 (2));
617fe022e965 factor.m: Return output factors of same class as input for Matlab compatibility. Rik <rik@octave.org> parents: 25269 diff changeset	234 %! assert (n, double (3));
617fe022e965 factor.m: Return output factors of same class as input for Matlab compatibility. Rik <rik@octave.org> parents: 25269 diff changeset	235
19833 9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions. Rik <rik@octave.org> parents: 19697 diff changeset	236 ## Test input validation
28896 90fea9cc9caa test: Add expected error message <Invalid call> to BIST tests for nargin. Rik <rik@octave.org> parents: 28891 diff changeset	237 %!error <Invalid call> factor ()
29983 ecbcc4647dbe factor.m: Overhaul function to support inputs > flintmax. Rik <rik@octave.org> parents: 29359 diff changeset	238 %!error <Q must be a real non-negative integer> factor ([1,2])
25269 cac96fd5310d factor.m: Emit an error if input is negative (bug #53425). Dildar Sk <dildarsk101010@gmail.com> parents: 25054 diff changeset	239 %!error <Q must be a real non-negative integer> factor (6i)
29983 ecbcc4647dbe factor.m: Overhaul function to support inputs > flintmax. Rik <rik@octave.org> parents: 29359 diff changeset	240 %!error <Q must be a real non-negative integer> factor (-20)
25269 cac96fd5310d factor.m: Emit an error if input is negative (bug #53425). Dildar Sk <dildarsk101010@gmail.com> parents: 25054 diff changeset	241 %!error <Q must be a real non-negative integer> factor (1.5)
29983 ecbcc4647dbe factor.m: Overhaul function to support inputs > flintmax. Rik <rik@octave.org> parents: 29359 diff changeset	242 %!error <Q too large to factor> factor (flintmax ("single") + 2)
ecbcc4647dbe factor.m: Overhaul function to support inputs > flintmax. Rik <rik@octave.org> parents: 29359 diff changeset	243 %!error <Q too large to factor> factor (flintmax ("double") + 2)

Mercurial > octave

annotate scripts/specfun/factor.m @ 31107:bbb59cc6698c stable