Prime decomposition

From Rosetta Code

The prime decomposition of a number is defined as a list of prime numbers which when all multiplied together, are equal to that number. Example: 12 = 2 * 2 * 3, so its prime decomposition is {2, 2, 3}

Write a function which returns an array or collection which contains the prime decomposition of a given number, n, greater than 1. If your language does not have an isPrime-like function available, you may assume that you have a function which determines whether a number is prime (note its name before your code).

If you would like to test code from this task, you may use code from trial division or the Sieve of Eratosthenes.

Note: The program must not be limited by the word size of your computer or some other artificial limit; it should work for any number regardless of size (ignoring the physical limits of RAM etc).

Contents 1 Ada 2 C 3 C++ 4 Common Lisp 5 Forth 6 Haskell 7 J 8 Java 9 JavaScript 10 Logo 11 OCaml 12 Python 13 Smalltalk 14 Scheme 15 V

[edit] Ada

 
with Ada.Text_IO;  use Ada.Text_IO;
 
procedure Test_Prime is
   generic
      type Number is private;
      Zero : Number;
      One  : Number;
      Two  : Number;
      with function Image (X : Number) return String is <>;
      with function "+"   (X, Y : Number) return Number is <>;
      with function "/"   (X, Y : Number) return Number is <>;
      with function "mod" (X, Y : Number) return Number is <>;
      with function ">="  (X, Y : Number) return Boolean is <>;
   package Prime_Numbers is
      type Number_List is array (Positive range <>) of Number;
      function Decompose (N : Number) return Number_List;
      procedure Put (List : Number_List);
   end Prime_Numbers;
 
   package body Prime_Numbers is
      function Decompose (N : Number) return Number_List is
         Size : Natural := 0;
         M    : Number  := N;
         K    : Number  := Two;
      begin
         -- Estimation of the result length from above
         while M >= Two loop
            M := (M + One) / Two;
            Size := Size + 1;
         end loop;
         M := N;
         -- Filling the result with prime numbers
         declare
            Result : Number_List (1..Size);
            Index  : Positive := 1;
         begin
            while N >= K loop -- Divisors loop
               while Zero = (M mod K) loop -- While divides
                  Result (Index) := K;
                  Index := Index + 1;
                  M := M / K;
               end loop;
               K := K + One;
            end loop;
            return Result (1..Index - 1);
         end;
      end Decompose;
 
      procedure Put (List : Number_List) is
      begin
         for Index in List'Range loop
            Put (Image (List (Index)));
         end loop;
      end Put;
   end Prime_Numbers;
   package Integer_Numbers is new Prime_Numbers (Natural, 0, 1, 2, Positive'Image);
   use Integer_Numbers;
begin
   Put (Decompose (12));
end Test_Prime;
 

The solution is generic. The package is instantiated by a type that supports necessary operations +, /, mod, >=. The constants 0, 1, 2 are parameters too, because the type might have no literals. The package also provides a procedure to output an array of prime numbers and a function to convert a number to string (as a parameter). The function Decompose first estimates the maximal result length as log₂ of the argument. Then it allocates the result and starts to enumerate divisors. It does not care to check if the divisors are prime, because non-prime divisors will be automatically excluded. In the example provided, the package is instantiated with plain integer type. Sample output:

 2 2 3

[edit] C

Library: GMP

#include <stdio.h>
#include <stdlib.h>
#include <GMP/gmp.h>

void decompose(mpz_t n) {
  int i;
  mpz_t tmp, d;
  
  i = 0;
  mpz_init(tmp);
  mpz_init(d);
  
  while(mpz_cmp_si(n, 1)) {
    mpz_set_ui(d, 1);
    do {
      mpz_add_ui(tmp, d, 1);
      mpz_swap(tmp, d);
    } while(!mpz_divisible_p(n, d));
    mpz_divexact(tmp, n, d);
    mpz_swap(tmp, n);
    gmp_printf("%s%Zd", i++?" * ":"", d);
  }
  gmp_printf("\n");
}

int main(int argc, char ** argv) {
  mpz_t n;
  
  if(argc != 2) {
    puts("Pass a parameter");
    return EXIT_SUCCESS;
  }
  
  mpz_init_set_str(n, argv[1], 10);
  decompose(n);
  mpz_clear(n);
  
  return EXIT_SUCCESS;
}

[edit] C++

Works with: g++ version 4.1.2 20061115 (prerelease) (Debian 4.1.1-21)

Library: GMP

 
#include <iostream>
#include <gmpxx.h>
 
// This function template works for any type representing integers or
// nonnegative integers, and has the standard operator overloads for
// arithmetic and comparison operators, as well as explicit conversion
// from int.
//
// OutputIterator must be an output iterator with value_type Integer. 
// It receives the prime factors.
template<typename Integer, typename OutputIterator>
 void decompose(Integer n, OutputIterator out)
{
  Integer i(2);
 
  while (n != 1)
  {
    while (n % i == Integer(0))
    {
      *out++ = i;
      n /= i;
    }
    ++i;
  }
}
 
// this is an output iterator similar to std::ostream_iterator, except
// that it outputs the separation string *before* the value, but not
// before the first value (i.e. it produces an infix notation).
template<typename T> class infix_ostream_iterator:
  public std::iterator<T, std::output_iterator_tag>
{
  class Proxy;
  friend class Proxy;
  class Proxy
  {
  public:
    Proxy(infix_ostream_iterator& iter): iterator(iter) {}
    Proxy& operator=(T const& value)
    {
      if (!iterator.first)
      {
        iterator.stream << iterator.infix;
      }
      iterator.stream << value;
    }
  private:
    infix_ostream_iterator& iterator;
  };
public:
  infix_ostream_iterator(std::ostream& os, char const* inf):
    stream(os),
    first(true),
    infix(inf)
  {
  }
  infix_ostream_iterator& operator++() { first = false; return *this; }
  infix_ostream_iterator operator++(int)
  {
    infix_ostream_iterator prev(*this);
    ++*this;
    return prev;
  }
  Proxy operator*() { return Proxy(*this); }
private:
  std::ostream& stream;
  bool first;
  char const* infix;
};
 
int main()
{
  std::cout << "please enter a positive number: ";
  mpz_class number;
  std::cin >> number;
 
  if (number <= 0)
    std::cout << "this number is not positive!\n;";
  else
  {
    std::cout << "decomposition: ";
    decompose(number, infix_ostream_iterator<mpz_class>(std::cout, " * "));
    std::cout << "\n";
  }
}
 

[edit] Common Lisp

(defun prime-decomposition (n)
  (when (> n 1)
    (loop for d from 2 upto n
          until (zerop (mod n d))
          finally (return (cons d (prime-decomposition (truncate n d)))))))

[edit] Forth

: decomp ( n -- )
  2
  begin  2dup dup * >=
  while  2dup /mod swap
         if   drop  1+ 1 or    \ next odd number
         else -rot nip  dup .
         then
  repeat
  drop . ;

[edit] Haskell

primes = sieve [2..]
  where
    sieve (p:xs) = p : sieve [x|x <- xs, x `mod` p > 0]

factorize n pps@(p:ps) = case n `divMod` p of
  (0,1)         -> []
  (remainder,0) -> p : factorize remainder pps
  _             -> factorize n ps

[edit] J

q:

Example use:

   q: 3684
2 2 3 307

   _1+2^128x
340282366920938463463374607431768211455
   q: _1+2^128x
3 5 17 257 641 65537 274177 6700417 67280421310721
   */ q: _1+2^128x
340282366920938463463374607431768211455

[edit] Java

Works with: Java version 1.5+

This is a version for arbitrary-precision integers which assumes the existence of a function with the signature:

public boolean prime(BigInteger i);

You will need to import java.util.List, java.util.LinkedList, and java.math.BigInteger.

public static List<BigInteger> primeFactorBig(BigInteger a){
    List<BigInteger> ans = new LinkedList<BigInteger>();
    //loop until we test the number itself or the number is 1
    for (BigInteger i = BigInteger.valueOf(2); i.compareTo(a) <= 0 && !a.equals(BigInteger.ONE);
         i = i.add(BigInteger.ONE)){
        while (a.remainder(i).equals(BigInteger.ZERO) && prime(i)) { //if we have a prime factor
            ans.add(i); //put it in the list
            a = a.divide(i); //factor it out of the number
        }
    }
    return ans;
}

[edit] JavaScript

This code uses the BigInteger Library jsbn and jsbn2
http://xenon.stanford.edu/~tjw/jsbn/jsbn.js
http://xenon.stanford.edu/~tjw/jsbn/jsbn2.js

  function run_factorize(input, output) {
     var n = new BigInteger(input.value, 10);
     var prod = false;
     var two = new BigInteger("2", 10);
     var divisor = new BigInteger("3", 10);
     
     if(n.compareTo(two) < 0) { return; }
     
     output.value = "";
     
     while(true) {
        var qr = n.divideAndRemainder(two);
        if(qr[1].compareTo(BigInteger.ZERO) == 0) {
           if(prod) { output.value += "*"; } else { prod = true; }
           output.value += "2";
           n = qr[0];
        }
        else { break; }
     }
     
     while(n.compareTo(BigInteger.ONE) != 0) {
        var qr = n.divideAndRemainder(divisor);
        if(qr[1].compareTo(BigInteger.ZERO) == 0) {
           if(prod) { output.value += "*"; } else { prod = true; }
           output.value += divisor;
           n = qr[0];
        }
        else { divisor = divisor.add(two); }
     }
  }

[edit] Logo

to decomp :p :n
  if :p*:p > :n [output (list :n)]
  if less? 0 modulo :n :p [output decomp :p+1 :n]
  output fput :p decomp :p :n/:p
end
to decompose :n
  output decomp 2 :n
end

[edit] OCaml

open Big_int;;
 
let prime_decomposition x =
  let rec inner c p =
    if lt_big_int p (square_big_int c) then
      [p]
    else if eq_big_int (mod_big_int p c) zero_big_int then
      c :: inner c (div_big_int p c)
    else
      inner (succ_big_int c) p
  in
    inner (succ_big_int (succ_big_int zero_big_int)) x;;
 

[edit] Python

import sys, time
 
def primes():
    yield 2
    n = 3
    while n < sys.maxint - 2:
        yield n
        n += 2
        while n < sys.maxint - 2 and len(zip(range(2),decompose(n))) > 1:
            n += 2
 
def decompose(n):
    for p in primes():
        if p*p > n: break
        while n % p == 0:
            yield p
            n /=p
    if n > 1:
        yield n
 
# Example: calculate factors of Mersenne numbers to M59 #
for m in xrange(1,60):
    if len(list(decompose(m))) != 1: continue
    p = 2 ** m - 1
    print "2**%d-1 = %d, with factors:"%(m, p),
    start = time.time()
    for factor in decompose(p):
        print factor,
        sys.stdout.flush()
    print "=> %.2fs"%(time.time()-start)
 

Output:

2**2-1 = 3, with factors: 3 => 0.00s
2**3-1 = 7, with factors: 7 => 0.00s
2**5-1 = 31, with factors: 31 => 0.00s
2**7-1 = 127, with factors: 127 => 0.00s
2**11-1 = 2047, with factors: 23 89 => 0.00s
2**13-1 = 8191, with factors: 8191 => 0.00s
2**17-1 = 131071, with factors: 131071 => 0.02s
2**19-1 = 524287, with factors: 524287 => 0.05s
2**23-1 = 8388607, with factors: 47 178481 => 0.02s
2**29-1 = 536870911, with factors: 233 1103 2089 => 0.09s
2**31-1 = 2147483647, with factors: 2147483647 => 33.57s
2**37-1 = 137438953471, with factors: 223 616318177 => 12.12s
2**41-1 = 2199023255551, with factors: 13367 164511353 => 4.50s
2**43-1 = 8796093022207, with factors: 431 9719 2099863 => 2.71s
2**47-1 = 140737488355327, with factors: 2351 4513 13264529 => 0.79s
2**53-1 = 9007199254740991, with factors: 6361 69431 20394401 => 65.28s
2**59-1 = 576460752303423487, with factors: 179951 3203431780337 => 14324.83s

[edit] Smalltalk

Integer extend [
    primesDo: aBlock [
        | div next rest |
        div := 2. next := 3.
        rest := self.
        [ [ rest \\ div == 0 ]
              whileTrue: [
                  aBlock value: div.
                  rest := rest // div ].
          rest = 1] whileFalse: [
              div := next. next := next + 2 ]
    ]
]
123456 primesDo: [ :each | each printNl ]

[edit] Scheme

(define (prime-decomposition x)
  (let inner ((c 2) (p x))
    (if (> (* c c) p)
        (list p)
        (if (zero? (remainder p c))
          (cons c (inner c (quotient p c)))
          (inner (+ c 1) p)))))

[edit] V

like in scheme (using variables)

[prime-decomposition
   [inner [c p] let
       [c c * p >]
           [p unit]
           [ [p c % zero?]
                   [c c p c / inner cons]
                   [c 1 + p inner]
             ifte]
       ifte].
   2 swap inner].

(mostly) the same thing using stack (with out variables)

[prime-decomposition
   [inner
       [dup * <]
           [pop unit]
           [ [% zero?]
                   [ [p c : [c p c / c]] view i inner cons]
                   [succ inner]
             ifte]
       ifte].
   2 inner].

Using it

|1221 prime-decomposition puts
=[3 11 37]