changelog shortlog tags changeset files revisions annotate raw

scripts/specfun/factor.m

changeset 10289: 4b124317dc38
parent:c1fff751b5a8
author: John W. Eaton <jwe@octave.org>
date: Tue Feb 09 20:58:55 2010 -0500 (45 minutes ago)
permissions: -rw-r--r--
description: base_properties::set_children: account for hidden children 
     1## Copyright (C) 2000, 2006, 2007, 2009 Paul Kienzle

     2##

     3## This file is part of Octave.

     4##

     5## Octave is free software; you can redistribute it and/or modify it

     6## under the terms of the GNU General Public License as published by

     7## the Free Software Foundation; either version 3 of the License, or (at

     8## your option) any later version.

     9##

    10## Octave is distributed in the hope that it will be useful, but

    11## WITHOUT ANY WARRANTY; without even the implied warranty of

    12## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

    13## General Public License for more details.

    14##

    15## You should have received a copy of the GNU General Public License

    16## along with Octave; see the file COPYING.  If not, see

    17## <http://www.gnu.org/licenses/>.

    18

    19## -*- texinfo -*-

    20## @deftypefn  {Function File} {@var{p} =} factor (@var{q})

    21## @deftypefnx {Function File} {[@var{p}, @var{n}] =} factor (@var{q})

    22##

    23## Return prime factorization of @var{q}.  That is, @code{prod (@var{p})

    24## == @var{q}} and every element of @var{p} is a prime number.  If

    25## @code{@var{q} == 1}, returns 1. 

    26##

    27## With two output arguments, return the unique primes @var{p} and

    28## their multiplicities.  That is, @code{prod (@var{p} .^ @var{n}) ==

    29## @var{q}}.

    30## @seealso{gcd, lcm}

    31## @end deftypefn

    32

    33## Author: Paul Kienzle

    34

    35## 2002-01-28 Paul Kienzle

    36## * remove recursion; only check existing primes for multiplicity > 1

    37## * return multiplicity as suggested by Dirk Laurie

    38## * add error handling

    39

    40function [x, m] = factor (n)

    41

    42  if (nargin < 1)

    43    print_usage ();

    44  endif

    45

    46  if (! isscalar (n) || n != fix (n))

    47    error ("factor: n must be a scalar integer");

    48  endif

    49

    50  ## Special case of no primes less than sqrt(n).

    51  if (n < 4)

    52    x = n;

    53    m = 1;

    54    return;

    55  endif 

    56

    57  x = [];

    58  ## There is at most one prime greater than sqrt(n), and if it exists,

    59  ## it has multiplicity 1, so no need to consider any factors greater

    60  ## than sqrt(n) directly. [If there were two factors p1, p2 > sqrt(n),

    61  ## then n >= p1*p2 > sqrt(n)*sqrt(n) == n. Contradiction.]

    62  p = primes (sqrt (n));

    63  while (n > 1)

    64    ## Find prime factors in remaining n.

    65    q = n ./ p;

    66    p = p (q == fix (q));

    67    if (isempty (p))

    68      ## Can't be reduced further, so n must itself be a prime.

    69      p = n;

    70    endif

    71    x = [x, p];

    72    ## Reduce n.

    73    n = n / prod (p);

    74  endwhile

    75  x = sort (x);

    76

    77  ## Determine muliplicity.

    78  if (nargout > 1)

    79    idx = find ([0, x] != [x, 0]);

    80    x = x(idx(1:length(idx)-1));

    81    m = diff (idx);

    82  endif

    83

    84endfunction

    85

    86## test:

    87##   assert(factor(1),1);

    88##   for i=2:20

    89##      p = factor(i);

    90##      assert(prod(p),i);

    91##      assert(all(isprime(p)));

    92##      [p,n] = factor(i);

    93##      assert(prod(p.^n),i);

    94##      assert(all([0,p]!=[p,0]));

    95##   end