Finding divisors from factors

[danaj]

Assuming I have a (sorted) list of prime factors of a number, I would like to create a sorted list of divisors. How can we do this as fast as possible with Perl code?

For example:

use ntheory qw/:all/;
my $n = 183092192580;
print "Factors:  @{[factor($n)]}\n";
print "Divisors: @{[divisors($n)]}\n";

produces

Factors:  2 2 3 5 11 277412413
Divisors: 1 2 3 4 5 6 10 11 12 15 20 22 30 33 44 55 60 66 110 132 165 220 330 660 277412413 554824826 832237239 1109649652 1387062065 1664474478 2774124130 3051536543 3328948956 4161186195 5548248260 6103073086 8322372390 9154609629 12206146172 15257682715 16644744780 18309219258 30515365430 36618438516 45773048145 61030730860 91546096290 183092192580

We can view this as the unique list of products of the powersets. For instance using the simple powerset from RosettaCode yields this rather slow method:

sub divisors {
  my($n,@factors) = @_;
  my %divisors;
  foreach my $l ( powerset(@factors) ) {
    my $d = 1;
    $d *= $_ for @$l;
    undef $divisors{$d};
  }
  my @d = sort { $a<=>$b } keys %divisors;
  @d;
}
sub powerset {@_ ? map { $_,[$_[0], @$_] } powerset(@_[1..$#_]) : [];}

Here is one possibility that runs about 4x faster:

sub divisors {
  my($n,@factors) = @_;

  my %all_factors;
  foreach my $f (@factors) {
    my @to_add = grep { $_ < $n }
                 map { $f * $_ }
                 keys %all_factors;
    undef @all_factors{ $f, @to_add };
  }
  # 3. Add 1 and n, sort.
  undef $all_factors{1};
  undef $all_factors{$n};
  my @divisors = sort {$a<=>$b} keys %all_factors;
  @divisors;
}

Note that the divisors are symmetric about sqrt(n) so if we had the bottom portion 1 .. sqrt(n) in the sorted array @d then we could quickly fill the rest using push @d, map { $_*$_ == $n ? () : int($n/$_) } reverse @d; This can be crudely applied to the solution above to yield this slightly faster solution:

sub divisors {
  my($n,@factors) = @_;
  my $sqrtn = int(sqrt($n));

  my %all_factors;
  foreach my $f ( grep { $_ <= $sqrtn } @factors) {
    my @to_add = grep { $_ <= $sqrtn }
                 map { $f * $_ }
                 keys %all_factors;
    undef @all_factors{ $f, @to_add };
  }
  undef $all_factors{1};
  my @d = sort {$a<=>$b} keys %all_factors;
  @d, map { $_*$_ == $n ? () : int($n/$_) } reverse @d;
}

Anything faster and/or simpler? We can use something crude like this to test:

#!/usr/bin/env perl
use warnings;
use strict;
use ntheory qw/factor/;
srand(1);

for (1..100000) {
  my $n = int(rand(2**36));
  my @d = divisors($n,factor($n));
  print "$n @d\n";
}

run as time perl test.pl | md5sum. The factoring and printing take some time, but less than 15% of the total for the functions above.

[LanX]

First of all powerset is wrong cause equal factors will produce duplicates, e.g.

P{2,2} = { (,),(2,),(,2),(2,2)}

Then I doubt you can multiply in a way which produces all divisors in a sorted way, so final sorting can't be avoided.

I'd try a recursive/iterative approach,

1. initialize the result set D with 1
2. shift the smallest prime p' from P
3. multiply all divisors in the temporary result set D with p' to get D'
4. push the new divisors D' into D
5. continue with step 2 till P empty
6. sort D
but note you'll need a test if p' and p'' are duplicates to avoid the aforementioned problem! :)

update

looks good! :)

use warnings;
use strict;

my @P = qw/2 2 3 5 11 277412413/;
@P= sort { $a<=>$b } @P;              # sort to be sure

my @D     = (1);
my @D_new = ();
my $p_old = 0;

for my $p (@P) {
    if ( $p == $p_old ) {
        @D_new = map { $_* $p } @D_new;
    } else {
        @D_new = map { $_* $p } @D;
    }
    push @D,@D_new;
    $p_old = $p;
}

@D= sort { $a<=>$b } @D;

print "Factors @P\n";
print "Divisors: @D\n";

out

Factors 2 2 3 5 11 277412413
Divisors: 1 2 3 4 5 6 10 11 12 15 20 22 30 33 44 55 60 66 110 132 165 
[phone] 277412413 554824826 832237239 1109649652 1387062065 16644
[phone] 3051536543 3328948956 4161186195 5548248260 61030730
[phone] 9154609629 12206146172 15257682715 16644744780 18309219
[phone] 36618438516 45773048145 61030730860 91546096290 18309
+2192580

[danaj]

Thanks Rolf.

The powerset solution was mainly to point out this simple way. It does work -- one just needs to remove duplicates using a hash.

It looks like your solution is very similar to my followup, we just do the multiply through a little different. The time is pretty close, and both faster than my earlier solutions.

They also have the advantage of not doing excess computation, which is important when we move to bigints where every operation is expensive (with Math::BigInt at least).

[LanX]

Of course you can still speed up this solution, but I doubt it would be more than micro optimizations or just related to Perl specific features/bottlenecks (like overhead for allocating temporary arrays) °

This needs only one multiplication per divisor, hence its O(#D) i.e. <= O(2^#P) for worst case, that's pretty optimal.

But I like the simplicity and readability compared to other suggestions, so personally I wouldn't start micro optimization now just to save 10 or 20 %.

update

Of course I'm not comparing to solutions using XS modules. In that case just use C directly.

°) or moving data or linearizing loops

update

E.G. you could try to replace the maps with direct

push @D_new, ($_*$p) for @D

expressions!

This could should be faster... but I'm not very motivated to benchmark myself. :)

[roboticus]

danaj:

For my version, I took the powerset version, used tail recursion to turn it into a loop, and used a hash to remove duplicates continuously. I don't expect that it's all that fast with the extra hash operations and a sort at the end, but it's simple enough:

#!/usr/bin/perl
use strict;
use warnings;

my @all_factors = (1, 2, 2, 3, 5, 11, 277412413);
my @divisors = build_factors(@all_factors);
print join(", ", @divisors),"\n";

sub build_factors {
    my @orig = @_;
    my %new = (1=>0);

    while (my $factor = shift @orig) {
        @new{map{$factor*$_} keys %new}=0;
    }
    return sort {$a<=>$b} keys %new;
}

...roboticus

When your only tool is a hammer, all problems look like your thumb.

[tye]

I make no claim that this is the fastest approach. It is likely slower because rather than an optimized call to sort, it does merge sorts as it goes (which might be faster if I were doing those merge sorts in C more like sort works). It also likely uses more memory than many solutions.

Though it appears to run in an instant on the size of problems I've thrown at it.

I only post it because I liked that I was able to avoid having to detect duplicates. I don't bother to do the multiplication to compute a divisor if it would end up being a duplicate.

#!/usr/bin/perl -w
use strict;

my %f = ( 2 => 2, 3 => 1, 5 => 1, 11 => 1, 277412413 => 1 );
my @f = sort { $a <=> $b } keys %f;

my @d = @f;
my( %d, %x );
for my $i ( 0..$#f ) {
    my $f = $f[$i];
    ( $d{$f} = [(0)x@f] )->[$i]++;
    $x{$f} = $i;
}

for my $i ( 0..$#f ) {
    my $f = $f[$i];
    my $p = $f{$f};
    my @t = @d;
    for my $e ( 1..$f{$f} ) {
        my( @t2, @m, @n );
        while( @t ) {
            my $t = shift @t;
            merge( \@n, $t, \@m, 1 < $e ? \@d : () );
            push @n, $t
                if  1 == $e;
            if(     $d{$t}[$i] < $p
                &&  $x{$t} <= $i    # Avoid duplicates
            ) {
                my $m = $t*$f;
                push @m, $m;
                push @t2, $m;
                ( $d{$m} = [@{ $d{$t} }] )->[$i]++;
                $x{$m} = $i
                    if  ! $x{$m} || $x{$m} < $i;
            }
        }
        @t = @t2;
        merge( \@n, 0, \@m, 1 < $e ? \@d : () );
        @d = @n;
    }
}

print "Factors:";
for my $f ( @f ) {
    if( 1 < $f{$f} ) {
        print " $f^$f{$f}";
    } else {
        print " $f";
    }
}
print "\nDivisors: 1 @d\n";
exit;


sub merge {
    my( $to_av, $max, $from_av, @rest ) = @_;
    if( $max ) {
        while( @$from_av && $from_av->[0] < $max ) {
            my $next = shift @$from_av;
            merge( $to_av, $next, @rest )
                if  @rest;
            push @$to_av, $next;
        }
        merge( $to_av, $max, @rest )
            if  @rest;
    } elsif( @rest ) {
        while( @$from_av ) {
            my $next = shift @$from_av;
            merge( $to_av, $next, @rest );
            push @$to_av, $next;
        }
        merge( $to_av, 0, @rest );
    } else {
        push @$to_av, @$from_av;
        @$from_av = ();
    }
}

Read more... (892 Bytes)

[danaj]

> Though it appears to run in an instant on the size of problems I've thrown at it.

For more complication, we can use a pathological case:

#!/usr/bin/env perl
use warnings;
use strict;
use bigint;

my @f = (2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59);
my @d = divisors(1922760350154212639070,@f);
print "$_\n" for @d;

Generating the 131,072 divisors is taking my machine anywhere from 5 seconds to 50 seconds for the various solutions. The performance ordering is a little different for bigints than for native ints. Roboticus's simple loop is the fastest, with my non-sqrt hash solution just behind. The least memory is my sqrt hash solution but not by a lot.

[tye]

And here is the simplest way I would solve this:

#!/usr/bin/perl -w
use strict;

use Algorithm::Loops 'NestedLoops';

my %f = ( 2 => 2, 3 => 1, 5 => 1, 11 => 1, 277412413 => 1 );
my @f = sort { $b <=> $a } keys %f;

my @d = NestedLoops(
    [ map [ 0 .. $f{$_} ], @f ],
    sub {
        my $p = 1;
        for my $i ( 0 .. $#_ ) {
            $p *= $f[$i]
                for 1..$_[$i];
        }
        return $p;
    },
);
@d = sort { $a <=> $b } @d;
print "Divisors: @d\n";

Read more... (862 Bytes)

[danaj]

Munging my C code, here is one that is faster, though not simpler:

sub divisors {
  my($n,@f) = @_;
  my @d = (1);
  while (@f) {
    my $p = shift(@f);
    if (!@f || $p != $f[0]) {
      # Factor appears once.  Multiply through.
      push @d, map { $p * $_ } @d;
    } else {
      my $e = 1;
      do { shift(@f); $e++; } while @f && $p == $f[0];
      # Factor appears $e times.  Multiply through p, p^2, p^3, ...
      my @t = @d;
      for (1 .. $e) {
        @t = map { $p * $_ } @t;
        push @d, @t;
      }
    }
  }
  @d = sort { $a <=> $b } @d;
  @d;
}

I'm basically converting to factor/exponent form, then pushing the divisors on the list. Doing it by prime factor means we avoid duplication so no need for hashes. Improvements welcome.

An aside, we can see some expected results using something like:
perl -MMath::Pari=divisors -E 'for my $n (1..1_000_000) { my $d = divisors($n); print "$n @$d\n"; }'
or
perl -Mntheory=divisors -E 'for my $n (1..1_000_000) { my @d = divisors($n); print "$n @d\n"; }'