why does factor(1) return 1? shouldn't it return null?

I'm using factor() to find the prime factors of a range of numbers.

As far as I understand it (which may be very limited as I am new to this!), factor(n) should accept a value of n=1 as 1 is a "real, nonnegative integer scalar", but should only return prime factors.

I believe that 1 is not regarded as a prime.

Given this, shouldn't the function return an empty/null row vector?

Thanks.

[Andrew]

Hi All who have responded - as previously mentioned, I am new to MATLAB and also to this Answers forum, so I may not have followed the response procedures correctly - and I'm not sure if putting a general reply to all responses (so far) in this comment is the correct way to go. If not please accept my apologies in advance.

However, just to clarify:

I think all the three answers I have seen are helpful and correct, but could only accept one, so I accepted the first as it clarified the function for me.

In general, it seems like a documentation update would be very helpful as the returned value for factor(1) did puzzle me.

For context, I was using the function to find the number of unique prime factors for non negative integers and then count them, hence my confusion when I got a count of 1 for the value 1.

I agree that it wouldn't be a good idea to change the existing function due to potential impact on existing code unless there was some way to add a switch to return empty values for 0 and 1 (as there are no prime factors for these). However I think this would probably be overkill given that these values can easily be excluded by code if necessary.

Thanks to all for the rapid and helpful responses.

[John D'Errico]

A virtue of MATLAB is it is easily extensible, even trivially so. If the behavior of a function on a boundary case like this is a problem, then just write your own code, call it primeFactors. It could be as simple as this:

function Xfac = primeFactors(X)
  % primeFactors - generate the list of all prime factors of X
  
  if x >= 2
    Xfac = factor(X);
  else
    Xfac = [];
  end
  
end

And of course, primeFactors has your desired behavior, and is pedantically correct, in that 1 or 0 have no prime factorization under the fundamental theorem of arithmetic.

But as you recognize, changing the code would seem to make little sense. Even putting a switch in there, to return empty when the users desires that behavior for the edge cases seems needlessly complex. A comment in the help indicating the behavior on edge cases? Yes, of couse. That never hurts to remove any confusion.

Anyway, this is a good catch on your part, that the help needs to be updated to be made clear.

Answer 1

[dpb]

help factor

 factor - Prime factors

    Syntax
      f = factor(n)

    Input Arguments
      n - Input value
        real, nonnegative integer scalar

    Examples
      openExample('matlab/PrimeFactorsofDoubleIntegerValueExample')
      openExample('matlab/PrimeFactorsofUnsignedIntegerValueExample')

    See also isprime, primes

    Introduced in MATLAB before R2006a
    Documentation for factor
       doc factor

    Other uses of factor

       sym/factor

From factor, _"F = factor(x) returns all irreducible factors of x in vector F. If x is an integer, factor returns the prime factorization of x. "_

The two sentences are conflicting; the first description is what is actually the implementation while the second is true only for integers >=2. I think the doc should be clarified.

For your purposes, use

v=0:10;         % a range of integers
pfactors=arrayfun(@primes,v,'uni',0)

pfactors = 1×11 cell array
    {1×0 double}    {1×0 double}    {[2]}    {[2 3]}    {[2 3]}    {[2 3 5]}    {[2 3 5]}    {[2 3 5 7]}    {[2 3 5 7]}    {[2 3 5 7]}    {[2 3 5 7]}

Oooops! The above doesn't do what you want/need, either, since primes returns all primes less than the number, not the factors in the number.

You'll have to make up your own version to fix up the shortcut factor takes to maintain multiplication rule.

pfactors=arrayfun(@factors,v,'uni',0)

pfactors = 1×11 cell array
    {0×0 double}    {0×0 double}    {[2]}    {[3]}    {[2 2]}    {[5]}    {[2 3]}    {[7]}    {[2 2 2]}    {[3 3]}    {[2 5]}

function f = factors(n)
%FACTOR Prime factors.
%   FACTOR(N) returns a vector containing the prime factors of N
%   Unlike MATHWORKS version 0, 1 return the empty result, not the value

  f=factor(n);      % let the builtin do its thing first
  if n<2            % fixup the nonmath definition ones
    f=[];
    return
  end
end

[Stephen23]

I wonder whether part of the problem is the terminology. Users searching for a function named factor are possibly looking for a useful factorization of an integer, not specifically a prime factorization. One possible approach would be to define the default output as a canonical factorization whose product reproduces the input:

prod(factor(n))==n

where e.g.

factor(0) = []

factor(1) = [1]

factor(-1) = [-1]

factor(6) = [1,2,3]

factor(-6) = [-1,2,3]

This has the benefit of not requiring any special cases for positive and negative values, nor for 1. Users who want strictly prime factors could then request that explicitly via an option. This would make the behavior consistent across all nonzero integers and avoid the current situation where the documentation discusses prime factors but the edge cases (0 and 1) are inconsistent with that description.

[dpb]

I like, but doesn't it require factor(0) = 0 to achieve? The null would return [] instead of 0.

[Stephen23]

@dpb: yes, you are right:

prod([])

ans = 1

[dpb]

Clearly, I didn't check prod for prod([]) -- superficially that seems a peculiar result. It's documented, but no rationale supplied. I suppose it's the same idea as @John D'Errico noted above to not have a product return a noncomputable value requiring special-case coding.

[Steven Lord]

From https://en.wikipedia.org/wiki/Empty_set#Operations_on_the_empty_set:

"When speaking of the sum of the elements of a finite set, one is inevitably led to the convention that the sum of the elements of the empty set (the empty sum) is zero. The reason for this is that zero is the identity element for addition. Similarly, the product of the elements of the empty set (the empty product) should be considered to be one, since one is the identity element for multiplication."

If you implemented sum and prod using a loop in the natural way, 0 and 1 are kind of the "base cases".

s = loopsum([])

s = 0

p = loopprod([])

p = 1

function s = loopsum(A)
s = 0;
for k = 1:numel(A)
    s = s + A(k);
end
end
function p = loopprod(A)
p = 1;
for k = 1:numel(A)
    p = p * A(k);
end
end

[dpb]

@Steven Lord, thanks for the link to the doc; I was just coming back to add the comment to note that property for sum([]) as well...

Answer 2

[Steven Lord]

Technically, yes. If you want prod(factor(x)) to return x and the output of factor to only include prime numbers, having factor(1) return the empty matrix would be consistent.

But there would still be an inconsistency -- what should factor(0) return? 0 isn't prime. Having it error feels like it could be technically correct but practically annoying. [The fundamental theorem of arithmetic covers integer values greater than 1, so technically it doesn't apply to either of those cases.]

Doing a little archaeology and going back to the earliest version of factor.m that I have access to (over 30 years old) I believe factor(1) would actually error (with an incorrect error message.)

Almost 27 years ago, factor.m was modified to just return the input if that input is non-negative and is less than 4. That code is still present in factor.m today. This optimization works for n = 3 and n = 2. It's this that causes factor to return 1 for n = 1 and 0 for n = 0. [Negative values are checked for and cause factor to error out on an earlier line of code.]

dbtype 29:34 factor.m

29    if n < 4
30       f = floor(n); 
31       return
32    else
33       f = [];
34    end

We could change how factor.m behaves in the n = 1 case (and in the n = 0 case as well.) But given how long this behavior has been in the product I'd be hesitant to change it without a good use case and an analysis of the potential impact of the incompatibility.

[John D'Errico]

@Steven Lord - I agree, that factor as applied to 0 or 1 should remain as they are. But I do think the help should reflect that.

[Steven Lord]

That's fair. They are edge cases and as such it's reasonable to ask that they be documented.

Answer 3

[John D'Errico]

Yes, I agree, in that

f1 = factor(1)

f1 = 1

could arguably return empty. 1 is not indeed prime, and so the result is inconsistent with the documentation.

I'd also note that

f0 = factor(0)

f0 = 0

fails to return a set of prime factors. Again, 0 is not prime.

However, you could also argue that 1 is indeed a "factor" of the number 1. In fact, it is the only such number that divides 1, and is not itself reducible. It is just that 1 is technically not a prime factor of 1, which lacks any truly prime factors at all. It gets even more difficult to decide what should be returned for the factor(s) of 0. But the help for factor allows any non-negative arguement.

It is also true that the product of the vector returned by factor reproduces the input. And so we see that

prod(f0)

ans = 0

prod(f1)

ans = 1

do indeed reproduce the original number as factored. And this could be a useful property of factor, relied upon in legacy code, going back as much as 40 years now. Yes, factor has been doing this for a long time. And any change to factor would potentially cause legacy code to fail were they to arbitrarily change the code. That would mean any change to factor would need to happen with a long lead time. They would probably keep reminding us for multiple years to change any code that relies on factor(1).

Personally, I think the argument about what should be returned as the factors of either 0 or 1 as being more a pedantic thing, up there with the number of angels that fit on the head of a pin. I don't think it merits a need to change the code, because of the issues that might potentially cause.

Better I think is to just introduce a documentation change, indicating that factor(1) will return 1 as a special case (the same applies to factor(0)). That can happen quite quickly, probably seeing such a change reflected in R2027, or even sooner. A change to the documentation IS DEFINITELY a good idea here. I would recommend this be introduced to support. I'll post it as a request today.

[dpb]

I don't suppose the archeology uncovered any other rationale notes besides the shortcut, but I presumed that keeping the product of the results of factor() reproducing the inputs was the motivation. I agree also on only a documentation update as noted before. I hadn't ever thought of it specifically before as I don't know that had ever had reason to inspect factor(1) before...

[John D'Errico]

Funny thing is, I don't think it ever bothered me either. Returning a factorization that contains 1 (or 0) just seems to make more sense. Do I really want factor to generate an error message for a boundary case like this? Not really. And even returning an empty vector for factor(1) or factor(0), that now would fail the requirement that prod(factor(n))==N. So I sort of like the return of a non-prime for those boundary cases.

And that means, in my eyes, really the only right thing to do is to update the documentation. Just have it say that well, yeah, factor does not return primes for these special cases, but what can we do there?

[Sam Chak]

Thanks to everyone for their input on this matter. I agree that the documentation should be updated for the exceptions: 1, 0, -1. Prime factorization also extends to negative integers. However, it extracts the negative sign as an irreducible unit (-1), rather than treating it as a prime.

F = factor(sym(-6))

F = 

It is also intriguing to investigate the computational limits of the factor command. I test the special class of numbers called Mersenne primes, which take the form , where p is also a prime.

F = factor(sym(2^53 - 1))

F = 

It correctly returns 3 distinct prime factors.

I also test this case further. Consider the value 2^53 + 2. This is not a Mersenne prime, but it is known to produce 4 distinct prime factors: 2 × 17 × 858001 × 308761441. Yet, the factor command returns this result:

F = factor(sym(2^53 + 2))

F = 

The issue is likely related to flintmax, which equals 2^53, the largest consecutive integer in IEEE double precision. However, I feel it should not produce a run of 2s. Therefore, the documentation should specifically state that for every nonnegative integer x with , factor returns the prime factorization of x.

flintmax

ans = 9.0072e+15

[Steven Lord]

For numbers that large, don't operate using double. Use a 64 bit integer instead, but be careful how you compute it in the first place.

two = uint64(2);
x = two^53+two;
factor(x)

ans = 1×4

2 17 858001 308761441

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Or if you want to do this symbolically:

two = sym(2);
x = two^53+two;
factor(x)

ans =