http://qs1969.pair.com?node_id=482105


in reply to Algorithm for cancelling common factors between two lists of multiplicands

There are better ways to compute factorial quotients without overflow or loss of precision. But I won't go into that here. Instead I'll address the problem you stated originally.

After examining all of the replies to date, it seems you want to reduce the quotient. Since you're dealing with large numbers (products of factorials), reducing it won't be so easy. Here's a school-boy method that just relies on breaking the process down into steps, nothing fancy:

First in pseudo-code:

1) Factor each term in the numerator and denominator into its prime factors. (Since these are factorials, this is pretty quick).

2) Collect all of the numerator factors into one array, and all of the denominator factors in another array.

3) Sort both arrays.

4) Do a modified merge sort to cancel like terms in both arrays. Keep any terms that don't have a match in the other array.

5) You can compute the quotient to a high degree of accuracy without overflowing or underflowing by judicious choices in the next term to grab.

Now let's see if I can do this in code (untested):
my @num = ( 3, 7, 11 ); my @den = ( 5, 12, 19 ); # get prime factors (you have to write this yourself :) # it returns the list of prime factors for each argument my @num_fac = get_prime_factors_list( 2..$num[-1] ); my @den_fac = get_prime_factors_list( 2..$den[-1] ); @num_fac = sort { $a <=> $b } @num_fac; @den_fac = sort { $a <=> $b } @den_fac; # modified merge sort my $n, $d; while (( $n < $#num_fac ) and ( $d < $#den_fac )) { if ( $num_fac[$n] == $den_fac[$d] ) { # drop both $n++; $d++; } elsif ( $num_fac[$n] < $den_fac[$d] ) { push @num_fin, $num_fac[$n]; $n++; } else { push @den_fin, $den_fac[$d]; $d++; } } # any remaining factors are not shared push @num_fin, @num_fac[$n..$#num_fac]; push @den_fin, @den_fac[$d..$#den_fac]; # edge cases @num_fin = (1) unless @num_fin; @den_fin = (1) unless @den_fin; # carefully multiply and divide my $q = shift @num_fin; while ( @num_fin and @den_fin ) { if ( $q > 1 ) { $q = $q / (shift @den_fin); } if ( $q < 1 ) { $q = $q * (shift @num_fin); } } # only one of these will be true: $q = $q * (shift @num_fin) while @num_fin; $q = $q / (shift @den_fin) while @den_fin;
You can probably improve this considerably. For instance, instead of keeping each prime factor (including duplicates), only keep a count.

-QM
--
Quantum Mechanics: The dreams stuff is made of

  • Comment on Re: Algorithm for cancelling common factors between two lists of multiplicands
  • Download Code

Replies are listed 'Best First'.
Re^2: Algorithm for cancelling common factors between two lists of multiplicands
by BrowserUk (Patriarch) on Aug 09, 2005 at 09:06 UTC
    There are better ways to compute factorial quotients without overflow or loss of precision.

    So, now I've got my "schoolboy approach" working, care to point me at the methods clever people would use?


    Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
    Lingua non convalesco, consenesco et abolesco. -- Rule 1 has a caveat! -- Who broke the cabal?
    "Science is about questioning the status quo. Questioning authority".
    The "good enough" maybe good enough for the now, and perfection maybe unobtainable, but that should not preclude us from striving for perfection, when time, circumstance or desire allow.
      So, now I've got my "schoolboy approach" working, care to point me at the methods clever people would use?
      I was too lazy to look before, but now that I have, I was a bit off. On the other hand, I snuck in something similar in step 5 anyway. Here's what I remembered when I said that:

      In QOTW #2 MJD considers the n_choose_k function that computes:

      sub n_choose_k { my ($n, $k) = @_; # f($n) = $n! return f($n)/f($k)/f($n-$k); }
      He notes that the intermediate values are very large, even though the final result is not. He submits this replacement:
      sub n_choose_k { my ($n, $k) = @_; my $t = 1; for my $i (1 .. $k) { $t *= $n - $k + $i; $t /= $i; } return $t; }
      and he notes:
      $t here is always an integer, and never gets bigger than necessary.
      For your formula, this only applies to a small portion of the terms. You may be better off just ignoring MJD's improvement and proceeding as I outlined before.

      -QM
      --
      Quantum Mechanics: The dreams stuff is made of