Prime decomposition: Difference between revisions

From Rosetta Code
Content added Content deleted
No edit summary
No edit summary
Line 539: Line 539:
}
}
}</lang>
}</lang>
=={{header|C sharp|C#}}==
<lang csharp>using System;
using System.Collections.Generic;


namespace PrimeDecomposition
{
class Program
{
static void Main(string[] args)
{
getPrimes(12);
}

static List<int> getPrimes(decimal n)
{
List<int> storage = new List<int>();
while (n > 1)
{
int i = 1;
while (true)
{
if (isPrime(i))
{
if (((decimal)n / i) == Math.Round((decimal) n / i))
{
n /= i;
storage.Add(i);
break;
}
}
i++;
}
}
return storage;
}

static bool isPrime(int n)
{
if (n <= 1) return false;
for (int i = 2; i <= Math.Sqrt(n); i++)
if (n % i == 0) return false;
return true;
}
}
}</lang>
=={{header|Clojure}}==
=={{header|Clojure}}==
<lang lisp>(use 'clojure.contrib.lazy-seqs) ; defines lazy sequence 'primes
<lang lisp>(use 'clojure.contrib.lazy-seqs) ; defines lazy sequence 'primes

Revision as of 16:08, 7 June 2010

Task
Prime decomposition
You are encouraged to solve this task according to the task description, using any language you may know.

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).

Ada

<lang 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;</lang> 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 log2 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

ALGOL 68

Translation of: Python

- note: This specimen retains the original Python coding style.

Works with: ALGOL 68 version Revision 1 - no extensions to language used
Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny
Works with: ELLA ALGOL 68 version Any (with appropriate job cards) - tested with release 1.8-8d

<lang algol68>#IF long int possible THEN #

MODE LINT = LONG INT; LINT lmax int = long max int; OP LLENG = (INT i)LINT: LENG i, 

  LSHORTEN = (LINT i)INT: SHORTEN i;
  1. ELSE

MODE LINT = INT; LINT lmax int = max int; OP LLENG = (INT i)LINT: i, 

  LSHORTEN = (LINT i)INT: i;

FI #

OP LLONG = (INT i)LINT: LLENG i;

MODE YIELDLINT = PROC(LINT)VOID;

PROC (LINT, YIELDLINT)VOID decompose;

INT upb cache = bits width;

BITS cache := 2r0; BITS cached := 2r0;

PROC is prime = (LINT n)BOOL: ( 

   BOOL 
       has factor := FALSE,
       out := TRUE;
 # FOR LINT factor IN # decompose(n, # ) DO ( #
 ##   (LINT factor)VOID:(
     IF has factor THEN out := FALSE; GO TO done FI;
     has factor := TRUE
 # OD # ));
   done: out

);

PROC is prime cached := (LINT n)BOOL: ( 

   LINT l half n = n OVER LLONG 2 - LLONG 1;
   IF l half n <= LLENG upb cache THEN
       INT half n = LSHORTEN l half n;
       IF half n ELEM cached THEN
           BOOL(half n ELEM cache)
       ELSE
           BOOL out = is prime(n);
           BITS mask = 2r1 SHL (upb cache - half n);
           cached := cached OR mask;
           IF out THEN cache := cache OR mask FI;
           out
       FI
   ELSE
       is prime(n) # above useful cache limit #
   FI

);


PROC primes := (YIELDLINT yield)VOID:( 

   yield(LLONG 2);
   LINT n := LLONG 3;
   WHILE n < l maxint - LLONG 2 DO
       yield(n);
       n +:= LLONG 2;
       WHILE n < l maxint - LLONG 2 AND NOT is prime cached(n) DO
           n +:= LLONG 2
       OD
   OD

);

  1. PROC # decompose := (LINT in n, YIELDLINT yield)VOID: (
   LINT n := in n;
 # FOR LINT p IN # primes( # ) DO ( #
 ##   (LINT p)VOID:
       IF p*p > n THEN
           GO TO done
       ELSE
           WHILE n MOD p = LLONG 0 DO
               yield(p);
               n := n OVER p
           OD
       FI
 # OD #  );
   done:
   IF n > LLONG 1 THEN
       yield(n)
   FI

);

main:(

  1. FOR LINT m IN # primes( # ) DO ( #
    1. (LINT m)VOID:(
     LINT p = LLONG 2 ** LSHORTEN m - LLONG 1;
     print(("2**",whole(m,0),"-1 = ",whole(p,0),", with factors:"));
   # FOR LINT factor IN # decompose(p, # ) DO ( #
   ##   (LINT factor)VOID:
         print((" ",whole(factor,0)))
   # OD # );
     print(new line);
     IF m >= LLONG 59 THEN GO TO done FI
  1. OD # ));
 done: EMPTY

)</lang> Output: 

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

Note: ALGOL 68G took 49,109,599 BogoMI and ELLA ALGOL 68RS took 1,127,634 BogoMI to complete the example.

AutoHotkey

<lang AutoHotkey>MsgBox % factor(8388607)  ; 47 * 178481

factor(n) { 

 If (n = 1)
   Return
 f = 2
 While (f <= n)
 {
   If (Mod(n, f) = 0)
   {
     next := factor(n / f)
     factors = %f%`n%next%
     Return factors
   }
   f++
 }

}</lang>

AWK

As the examples show, pretty large numbers can be factored in tolerable time: <lang awk>$ awk 'func pfac(n){r="";f=2;while(f<=n){while(!(n%f)){n=n/f;r=r" "f};f=f+2-(f==2)};return r}{print pfac($1)}' 36 

2 2 3 3

77 

7 11

536870911 

233 1103 2089

8796093022207 

431 9719 2099863</lang>

C

Library: GMP

primedecompose.h

<lang c>#ifndef _PRIMEDECOMPOSE_H_

  1. define _PRIMEDECOMPOSE_H_
  2. include <gmp.h>

int decompose(mpz_t n, mpz_t *o);

  1. endif</lang>

primedecompose.c

<lang c>#include "primedecompose.h"

int decompose(mpz_t n, mpz_t *o) { 

 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);
   mpz_init(o[i]);
   mpz_set(o[i], d);
   i++;
 }
 return i;

}</lang>

Testing

<lang c>#include <stdio.h>

  1. include <stdlib.h>
  2. include <gmp.h>
  3. include "primedecompose.h"

mpz_t dest[100]; /* must be big enough to hold all the factors! */

int main(int argc, char **argv) { 

 mpz_t n;
 int i, l;

 if(argc != 2) {
   puts("Pass a parameter");
   return EXIT_SUCCESS;
 }

 mpz_init_set_str(n, argv[1], 10);
 l = decompose(n, dest);
 
 for(i=0; i < l; i++) {
   gmp_printf("%s%Zd", i?" * ":"", dest[i]);
   mpz_clear(dest[i]);
 }
 printf("\n");

 return EXIT_SUCCESS;

}</lang>

Using GNU Compiler Collection gcc extensions

Translation of: ALGOL 68
Works with: gcc version 4.3.0 20080428 (Red Hat 4.3.0-8)

Note: The following code sample is experimental as it implements python style iterators for (potentially) infinite sequences. C is not normally written this way, and in the case of this sample it requires the GCC "nested procedure" extension to the C language. <lang c>#include <limits.h>

  1. include <stdio.h>
  2. include <math.h>

typedef enum{false=0, true=1}bool; const int max_lint = LONG_MAX;

typedef long long int lint;

  1. assert sizeof_long_long_int (LONG_MAX>=8) /* XXX */
  1. ifdef NEED_GOTO
  2. include <setjmp.h>

/* declare label otherwise it is not visible in sub-scope */

  1. define LABEL(label) jmp_buf label; if(setjmp(label))goto label;
  2. define GOTO(label) longjmp(label, true)
  3. endif

/* the following line is the only time I have ever required "auto" */

  1. define FOR(i,iterator) auto bool lambda(i); yield_init = (void *)λ iterator; bool lambda(i)
  2. define DO {
  3. define YIELD(x) if(!yield(x))return
  4. define BREAK return false
  5. define CONTINUE return true
  6. define OD CONTINUE; }

/* Warning: _Most_ FOR(,){ } loops _must_ have a CONTINUE as the last statement. 

*   Otherwise the lambda will return random value from stack, and may terminate early */

typedef void iterator, lint_iterator; /* hint at procedure purpose */ static volatile void *yield_init; /* not thread safe */

  1. define YIELDS(type) bool (*yield)(type) = yield_init

typedef unsigned int bits;

  1. define ELEM(shift, bits) ( (bits >> shift) & 0b1 )

bits cache = 0b0, cached = 0b0; const lint upb_cache = 8 * sizeof(cache);

lint_iterator decompose(lint); /* forward declaration */

bool is_prime(lint n){ 

  bool has_factor = false, out = true;

/* for factor in decompose(n) do */ 

  FOR(lint factor, decompose(n)){
      if( has_factor ){ out = false; BREAK; }
      has_factor = true;
      CONTINUE;
  }
  return out;

}

bool is_prime_cached (lint n){ 

   lint half_n = n / 2 - 2;
   if( half_n <= upb_cache){
       /* dont cache the initial four, nor the even numbers */
       if (ELEM(half_n,cached)){
           return ELEM(half_n,cache);
       } else {
           bool out = is_prime(n);
           cache = cache | out << half_n;
           cached = cached | 0b1 << half_n;
           return out;
       }
   } else {
       return is_prime(n);
   }

}

lint_iterator primes (){ 

   YIELDS(lint);
   YIELD(2);
   lint n = 3;
   while( n < max_lint - 2 ){
       YIELD(n);
       n += 2;
       while( n < max_lint - 2 && ! is_prime_cached(n) ){
           n += 2;
       }
   }

}

lint_iterator decompose (lint in_n){ 

   YIELDS(lint);
   lint n = in_n;
/* for p in primes do */
   FOR(lint p, primes()){
       if( p*p > n ){
           BREAK;
       } else {
           while( n % p == 0 ){
               YIELD(p);
               n = n / p;
           }
       }
       CONTINUE;
   }
   if( n > 1 ){
       YIELD(n);
   }

}

main(){ 

   FOR(lint m, primes()){
       lint p = powl(2, m) - 1;
       printf("2**%lld-1 = %lld, with factors:",m,p);
       FOR(lint factor, decompose(p)){
           printf(" %lld",factor);
           fflush(stdout);
           CONTINUE;
       }
       printf("\n",m);
       if( m >= 59 )BREAK;
       CONTINUE;
   }

}</lang> Output: 

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

Note: gcc took 487,719 BogoMI to complete the example.

C++

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

<lang cpp>#include <iostream>

  1. 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";
 }

}</lang>

C#

<lang csharp>using System; using System.Collections.Generic;

namespace PrimeDecomposition { 

   class Program
   {
       static void Main(string[] args)
       {
           getPrimes(12);
       }

       static List<int> getPrimes(decimal n)
       {
           List<int> storage = new List<int>();
           while (n > 1)
           {
               int i = 1;
               while (true)
               {
                   if (isPrime(i))
                   {                        
                       if (((decimal)n / i) == Math.Round((decimal) n / i))
                       {
                           n /= i;
                           storage.Add(i);                            
                           break;
                       }
                   }
                   i++;
               }
           }
           return storage;
       }

       static bool isPrime(int n)
       {
           if (n <= 1) return false;
           for (int i = 2; i <= Math.Sqrt(n); i++)
               if (n % i == 0) return false;
           return true;
       }
   }

}</lang>

Clojure

<lang lisp>(use 'clojure.contrib.lazy-seqs) ; defines lazy sequence 'primes

(defn prime-factors [arg] 

 (assert (and (integer? arg) (>= arg 2)))
 (loop [pfs [], n arg]  ; pfs is the vector of prime factors already determined
   (if (= n 1)
     pfs
     (let [dps (for [p primes :while (<= (* p p) n) :when (zero? (rem n p))] p)
           ps  (for [p dps, q (rest (iterate #(/ % p) n)) :while (integer? q)] p)]
       (if (empty? dps)
         (recur (conj pfs n), 1)
         (recur (into pfs ps), (apply / n ps)))))))</lang>

Common Lisp

<lang Lisp>;;; Recursive algorithm (defun factor (n) 

 "Return a list of factors of N."
 (when (> n 1)
   (loop with max-d = (isqrt n)

for d = 2 then (if (evenp d) (+ d 1) (+ d 2)) do (cond ((> d max-d) (return (list n))) ; n is prime ((zerop (rem n d)) (return (cons d (factor (truncate n d)))))))))</lang>

E

This example assumes a function isPrime and was tested with this one. It could use a self-referential implementation such as the Python task, but the original author of this example did not like the ordering dependency involved.

<lang e>def primes := { 

   var primesCache := [2]
   /** A collection of all prime numbers. */
   def primes {
       to iterate(f) {
           primesCache.iterate(f)
           for x in (int > primesCache.last()) {
               if (isPrime(x)) {
                   f(primesCache.size(), x)
                   primesCache with= x
               }
           }
       }
   }

}

def primeDecomposition(var x :(int > 0)) { 

   var factors := []
   for p in primes {
       while (x % p <=> 0) {
           factors with= p
           x //= p
       }
       if (x <=> 1) {
           break
       }
   }
   return factors

}</lang>

Erlang

<lang erlang>factors(N) -> 

    factors(N,2).

factors(1,_) -> []; factors(N,K) when N rem K == 0 -> 

   [K|factors(N div K,K)];

factors(N,K) -> 

   factors(N,K+1).

</lang>

Factor

Factor already offers this functionality in its standard library. Example use: <lang factor>USING: math.primes.factors ; 12 factors . { 2 2 3 } 576460752303423487 factors . { 179951 3203431780337 }</lang>

FALSE

<lang false>[2[\$@$$*@>~][\$@$@$@$@\/*=$[%$." "$@\/\0~]?~[1+1|]?]#%.]d: 27720d;! {2 2 2 3 3 5 7 11}</lang>

Forth

<lang 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 . ;</lang>

Fortran

Works with: Fortran version 90 and later

<lang fortran>module PrimeDecompose 

 implicit none

 integer, parameter :: huge = selected_int_kind(18)
 ! => integer(8) ... more fails on my 32 bit machine with gfortran(gcc) 4.3.2

contains 

 subroutine find_factors(n, d)
   integer(huge), intent(in) :: n
   integer, dimension(:), intent(out) :: d

   integer(huge) :: div, next, rest
   integer :: i

   i = 1
   div = 2; next = 3; rest = n
   
   do while ( rest /= 1 )
      do while ( mod(rest, div) == 0 ) 
         d(i) = div
         i = i + 1
         rest = rest / div
      end do
      div = next
      next = next + 2
   end do

 end subroutine find_factors

end module PrimeDecompose</lang>

<lang fortran>program Primes 

 use PrimeDecompose
 implicit none

 integer, dimension(100) :: outprimes
 integer i

 outprimes = 0

 call find_factors(12345649494449_huge, outprimes)

 do i = 1, 100
    if ( outprimes(i) == 0 ) exit
    print *, outprimes(i)
 end do

end program Primes</lang>

Haskell

<lang 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</lang>

J

<lang j>q:</lang>

Example use: <lang j>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</lang>

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: <lang java>public boolean prime(BigInteger i);</lang> You will need to import java.util.List, java.util.LinkedList, and java.math.BigInteger. <lang java>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;

}</lang>

Alternate version, optimised to be faster. Uses java.util.List, java.util.LinkedList and java.math.BigInteger. <lang java>public boolean isPrime(BigInteger i);</lang> <lang java>private static final BigInteger two = BigInteger.valueOf(2); public List<BigInteger> primeDecomp(BigInteger a) { 

   //Is it even possible
   if(a.compareTo(two)<0) {
       return null; //impossible for values lower than 2
   }
   //quickly handle even values
   List<BigInteger> result = new LinkedList<BigInteger>();
   while(a.and(BigInteger.ONE).equals(BigInteger.ONE)) {
       a = a.shiftRight(1);
       result.add(two);
   }
   //left with odd values
   if(!a.equals(BigInteger.ONE)) {
       BigInteger b = BigInteger.valueOf(3);
       while(b.compareTo(a)<0) {
           if(isPrime(b)) {
               BigInteger[] dr = a.divideAndRemainder(b);
               if(dr[1].equals(BigInteger.ZERO)) {
                   result.add(b);
                   a = dr[0];
               }
           }
           b = b.add(two);
       }
       result.add(b); //b will always be prime here...
   }
   return result;

}</lang>

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 <lang javascript>function run_factorize(input, output) { 

   var n = new BigInteger(input.value, 10);
   var TWO = new BigInteger("2", 10);
   var divisor = new BigInteger("3", 10);
   var prod = false;

   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); 
   }

}</lang>

<lang logo>to decompose :n [:p 2] 

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

end</lang>

Mathematica

Bare built-in function does: <lang Mathematica> FactorInteger[2016] => {{2, 5}, {3, 2}, {7, 1}}</lang>

Read as: 2 to the power 5 times 3 squared times 7 (to the power 1). To show them nicely we could use the following functions: <lang Mathematica>supscript[x_,y_]:=If[y==1,x,Superscript[x,y]] ShowPrimeDecomposition[input_Integer]:=Print@@{input," = ",Sequence@@Riffle[supscript@@@FactorInteger[input]," "]}</lang>

Example for small prime: <lang Mathematica> ShowPrimeDecomposition[1337]</lang> gives: <lang Mathematica> 1337 = 7 191</lang> Examples for large primes: <lang Mathematica> Table[AbsoluteTiming[ShowPrimeDecomposition[2^a-1]]//Print[#1," sec"]&,{a,50,150,10}];</lang> gives back: <lang Mathematica>1125899906842623 = 3 11 31 251 601 1801 4051 0.000231 sec 1152921504606846975 = 3^2 5^2 7 11 13 31 41 61 151 331 1321 0.000146 sec 1180591620717411303423 = 3 11 31 43 71 127 281 86171 122921 0.001008 sec 1208925819614629174706175 = 3 5^2 11 17 31 41 257 61681 4278255361 0.000340 sec 1237940039285380274899124223 = 3^3 7 11 19 31 73 151 331 631 23311 18837001 0.000192 sec 1267650600228229401496703205375 = 3 5^3 11 31 41 101 251 601 1801 4051 8101 268501 0.000156 sec 1298074214633706907132624082305023 = 3 11^2 23 31 89 683 881 2971 3191 201961 48912491 0.001389 sec 1329227995784915872903807060280344575 = 3^2 5^2 7 11 13 17 31 41 61 151 241 331 1321 61681 4562284561 0.000374 sec 1361129467683753853853498429727072845823 = 3 11 31 131 2731 8191 409891 7623851 145295143558111 0.024249 sec 1393796574908163946345982392040522594123775 = 3 5^2 11 29 31 41 43 71 113 127 281 86171 122921 7416361 47392381 0.009419 sec 1427247692705959881058285969449495136382746623 = 3^2 7 11 31 151 251 331 601 1801 4051 100801 10567201 1133836730401 0.007705 sec</lang>

MATLAB

<lang Matlab>function [outputPrimeDecomposition] = primedecomposition(inputValue) 

  outputPrimeDecomposition = factor(inputValue);</lang>

Maxima

This example was written by a novice in Maxima. If you are familiar with Maxima, please review and edit this example and remove this message. If the example does not work and you cannot fix it, replace this message with {{incorrect|Maxima|description of problem as you see it}}. If the code is correct but unidiomatic and you cannot fix it, replace this message with {{improve|Maxima|description of how it should be improved}}.


<lang maxima>prime_decomposition(x) := map(first, ifactors(x))</lang>

The builtin function factor(integer) also returns the prime decomposition of an integer, but it returns it as a product expression rather than a collection.

<lang maxima>(%i62) prime_decomposition(10); (%o62) [2,5]</lang>

OCaml

<lang 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;;</lang>

Octave

<lang octave>r = factor(120202039393); disp(r);</lang>

PicoLisp

<lang PicoLisp>(de factor (N) 

  (let (D 2  L (1 2 2)  M (sqrt N)  Res)
     (while (>= M D)
        (cond
           ((=0 (% N D)) 
              (push 'Res D) 
              (setq M (sqrt (setq N (/ N D)))) )
           (L (inc 'D (pop 'L)))
           (T (inc 'D 4) 
              (setq L (2 4 2 4 6 2 6)) ) ) )
     (if (= 1 N) Res (cons N Res)) ) )

(factor 1361129467683753853853498429727072845823)</lang> Output: 

-> (145295143558111 7623851 409891 8191 2731 131 31 11 3)

PureBasic

Works with: PureBasic version 4.41

<lang PureBasic> CompilerIf #PB_Compiler_Debugger 

 CompilerError "Turn off the debugger if you want reasonable speed in this example."

CompilerEndIf

Define.q

Procedure Factor(Number, List Factors()) 

 Protected I = 3
 While Number % 2 = 0
   AddElement(Factors())
   Factors() = 2
   Number / 2
 Wend
 Protected Max = Number
 While I <= Max And Number > 1
   While Number % I = 0
     AddElement(Factors())
     Factors() = I
     Number/I
   Wend
   I + 2
 Wend

EndProcedure

Number = 9007199254740991 NewList Factors() time = ElapsedMilliseconds() Factor(Number, Factors()) time = ElapsedMilliseconds()-time S.s = "Factored " + Str(Number) + " in " + StrD(time/1000, 2) + " seconds." ForEach Factors() 

 S + #CRLF$ + Str(Factors())

Next MessageRequester("", S) </lang> 

Factored 9007199254740991 in 0.27 seconds.
6361
69431
20394401

Python

Works with: Python version 2.5.1

<lang python>import sys, time

def is_prime(n): 

   return zip((True, False), decompose(n))[-1][0]

class IsPrimeCached(dict): 

   def __missing__(self, n):
       r = is_prime(n)
       self[n] = r
       return r

is_prime_cached = IsPrimeCached()

def primes(): 

   yield 2
   n = 3
   while n < sys.maxint - 2:
       yield n
       n += 2
       while n < sys.maxint - 2 and not is_prime_cached[n]:
           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
  1. Example: calculate factors of Mersenne numbers to M59 #

for m in primes(): 

   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)
   if m >= 59: break</lang>

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.00s
2**19-1 = 524287, with factors: 524287 => 0.01s
2**23-1 = 8388607, with factors: 47 178481 => 0.00s
2**29-1 = 536870911, with factors: 233 1103 2089 => 0.01s
2**31-1 = 2147483647, with factors: 2147483647 => 1.67s
2**37-1 = 137438953471, with factors: 223 616318177 => 0.02s
2**41-1 = 2199023255551, with factors: 13367 164511353 => 0.01s
2**43-1 = 8796093022207, with factors: 431 9719 2099863 => 0.01s
2**47-1 = 140737488355327, with factors: 2351 4513 13264529 => 0.00s
2**53-1 = 9007199254740991, with factors: 6361 69431 20394401 => 1.17s
2**59-1 = 576460752303423487, with factors: 179951 3203431780337 => 211.07s

Note: Python took 740,238 BogoMI to complete the example.

Modifying the primes() and is_prime() functions as below increases performance. <lang python>def is_prime(n): 

  return all( y==n for y in decompose(n))

primelist = [2,3] max_tested = 3

def primes(): 

  global max_tested
  for n in primelist:
     yield n
  n = max_tested
  pmax = sys.maxint-2
  while n < pmax:
     n += 2
     while not is_prime(n) and n < pmax:
        n += 2
     if n < pmax:
        primelist.append(n)
        max_tested = n
        yield n</lang>

R

<lang R>findfactors <- function(n) { 

 d <- c()
 div <- 2; nxt <- 3; rest <- n
 while( rest != 1 ) {
   while( rest%%div == 0 ) {
     d <- c(d, div)
     rest <- floor(rest / div)
   }
   div <- nxt
   nxt <- nxt + 2
 }
 d

}

print(findfactors(1005025))</lang>


Ruby

<lang ruby># get prime decomposition of integer i

  1. this routine is terribly inefficient, but elegance rules :-)

def prime_factors(i) 

 v = (2..i-1).detect{|j| i % j == 0} 
 v ? ([v] + prime_factors(i/v)) : [i]

end

  1. example: decompose all possible Mersenne primes up to 2**31-1

(2..31).each do |i| 

 puts "prime_factors(#{2**i-1}): #{prime_factors(2**i-1).join(' ')}"

end</lang>

A more efficient version, and quite similar to the Integer#prime_division method added by the

Library: mathn.rb

package in the Ruby stdlib:

<lang ruby>require 'mathn' def prime_factors(n) 

 factors = []
 prime_number_generator = Prime.new
 p = prime_number_generator.next
 while p <= n
   q, r = n.divmod(p)
   if r == 0
     factors << p
     n = q
   elsif p**2 >= n
     break
   else
     p = prime_number_generator.next
   end
 end
 factors << n if n > 1
 factors

end

  1. example: decompose all possible Mersenne primes up to 2**31-1

results = [] png = Prime.new

require 'benchmark' Benchmark.bm(7) do |x| 

 begin
   i = png.next
   n  = 2**i- 1
   f = prime_factors(n)
   results << "%2d : %-20d : %s\n" % [i, n, f.inspect]
   x.report("new-#{i}") {prime_factors(n)}
   x.report("ruby-#{i}") {n.prime_division}
 end while i < 53

end puts results</lang> 

             user     system      total        real
new-2    0.000000   0.000000   0.000000 (  0.000000)
ruby-2   0.000000   0.000000   0.000000 (  0.000000)
new-3    0.000000   0.000000   0.000000 (  0.000000)
ruby-3   0.000000   0.000000   0.000000 (  0.000000)
new-5    0.000000   0.000000   0.000000 (  0.000000)
ruby-5   0.000000   0.000000   0.000000 (  0.000000)
new-7    0.000000   0.000000   0.000000 (  0.000000)
ruby-7   0.000000   0.000000   0.000000 (  0.000000)
new-11   0.000000   0.000000   0.000000 (  0.000000)
ruby-11  0.000000   0.000000   0.000000 (  0.000000)
new-13   0.000000   0.000000   0.000000 (  0.000000)
ruby-13  0.000000   0.000000   0.000000 (  0.002000)
new-17   0.000000   0.000000   0.000000 (  0.005000)
ruby-17  0.015000   0.000000   0.015000 (  0.007000)
new-19   0.016000   0.000000   0.016000 (  0.013000)
ruby-19  0.015000   0.000000   0.015000 (  0.013000)
new-23   0.000000   0.000000   0.000000 (  0.006000)
ruby-23  0.000000   0.000000   0.000000 (  0.006000)
new-29   0.047000   0.000000   0.047000 (  0.024000)
ruby-29  0.031000   0.000000   0.031000 (  0.025000)
new-31  13.016000   0.000000  13.016000 ( 13.390000)
ruby-31 13.437000   0.000000  13.437000 ( 13.578000)
new-37   5.031000   0.000000   5.031000 (  4.639000)
ruby-37  4.750000   0.000000   4.750000 (  4.648000)
new-41   1.594000   0.000000   1.594000 (  1.668000)
ruby-41  1.546000   0.000000   1.546000 (  1.568000)
new-43   0.844000   0.000000   0.844000 (  0.879000)
ruby-43  0.906000   0.000000   0.906000 (  0.914000)
new-47   0.265000   0.000000   0.265000 (  0.256000)
ruby-47  0.266000   0.000000   0.266000 (  0.240000)
new-53  27.938000   0.000000  27.938000 ( 28.369000)
ruby-53 28.562000   0.000000  28.562000 ( 28.227000)
 2 : 3                    : [3]
 3 : 7                    : [7]
 5 : 31                   : [31]
 7 : 127                  : [127]
11 : 2047                 : [23, 89]
13 : 8191                 : [8191]
17 : 131071               : [131071]
19 : 524287               : [524287]
23 : 8388607              : [47, 178481]
29 : 536870911            : [233, 1103, 2089]
31 : 2147483647           : [2147483647]
37 : 137438953471         : [223, 616318177]
41 : 2199023255551        : [13367, 164511353]
43 : 8796093022207        : [431, 9719, 2099863]
47 : 140737488355327      : [2351, 4513, 13264529]
53 : 9007199254740991     : [6361, 69431, 20394401]

Scala

Getting the prime factors does not require identifying prime numbers. Since the problems seems to ask for it, here is one version that does it:

<lang scala>class PrimeFactors(n: BigInt) extends Iterator[BigInt] { 

 val zero = BigInt(0)
 val one = BigInt(1)
 val two = BigInt(2)
 def isPrime(n: BigInt) = n.isProbablePrime(10)
 var currentN = n
 var prime = two

 def nextPrime =
   if (prime == two) {
     prime += one
   } else {
     prime += two
     while (!isPrime(prime)) {
       prime += two
       if (prime * prime > currentN)
         prime = currentN
     }
   }

 def next = {
   if (!hasNext)
     throw new NoSuchElementException("next on empty iterator")
     
   while(currentN % prime != zero) {
     nextPrime
   }
   currentN /= prime
   prime
 }

 def hasNext = currentN != one && currentN > zero

}</lang>

The method isProbablePrime(n) has a chance of 1 - 1/(2^n) of correctly identifying a prime. Next is a version that does not depend on identifying primes, and works with arbitrary integral numbers:

<lang scala>class PrimeFactors[N](n: N)(implicit num: Integral[N]) extends Iterator[N] { 

 import num._
 val two = one + one
 var currentN = n
 var divisor = two

 def next = {
   if (!hasNext)
     throw new NoSuchElementException("next on empty iterator")
     
   while(currentN % divisor != zero) {
     if (divisor == two)
       divisor += one
     else
       divisor += two
       
     if (divisor * divisor > currentN)
       divisor = currentN
   }
   currentN /= divisor
   divisor
 }

 def hasNext = currentN != one && currentN > zero

}</lang>

Both versions can be rather slow, as they accept arbitrarily big numbers, as requested. Test: 

scala> BigInt(2) to BigInt(30) filter (_ isProbablePrime 10) map (p => (p, BigInt(2).pow(p.toInt) - 1)) foreach {
     |   case (prime, n) => println("2**"+prime+"-1 = "+n+", with factors: "+new PrimeFactors(n).mkString(", "))
     | }
2**2-1 = 3, with factors: 3
2**3-1 = 7, with factors: 7
2**5-1 = 31, with factors: 31
2**7-1 = 127, with factors: 127
2**11-1 = 2047, with factors: 23, 89
2**13-1 = 8191, with factors: 8191
2**17-1 = 131071, with factors: 131071
2**19-1 = 524287, with factors: 524287
2**23-1 = 8388607, with factors: 47, 178481
2**29-1 = 536870911, with factors: 233, 1103, 2089
2**31-1 = 2147483647, with factors: 2147483647
2**37-1 = 137438953471, with factors: 223, 616318177
2**41-1 = 2199023255551, with factors: 13367, 164511353
2**43-1 = 8796093022207, with factors: 431, 9719, 2099863
2**47-1 = 140737488355327, with factors: 2351, 4513, 13264529
2**53-1 = 9007199254740991, with factors: 6361, 69431, 20394401
2**59-1 = 576460752303423487, with factors: 179951, 3203431780337

Scheme

<lang scheme>(define (factor number) 

 (define (*factor divisor number)
   (if (> (* divisor divisor) number)
       (list number)
       (if (= (modulo number divisor) 0)
           (cons divisor (*factor divisor (/ number divisor)))
           (*factor (+ divisor 1) number))))
 (*factor 2 number))

(display (factor 111111111111)) (newline)</lang> Output: 

(3 7 11 13 37 101 9901)

Seed7

<lang seed7>const func array integer: factorise (in var integer: number) is func 

 result
   var array integer: result is 0 times 0;
 local
   var integer: checker is 2;
 begin
   while checker * checker <= number do
     if number rem checker = 0 then
       result &:= [](checker);
       number := number div checker;
     else
       incr(checker);
     end if;
   end while;
   if number <> 1 then
     result &:= [](number);
   end if;
 end func;</lang>

Original source: [1]

Slate

Admittedly, this is just based on the Smalltalk entry below: <lang slate>n@(Integer traits) primesDo: block "Decomposes the Integer into primes, applying the block to each (in increasing order)." [| div next remaining | 

 div: 2.
 next: 3.
 remaining: n.
 [[(remaining \\ div) isZero]
    whileTrue:
      [block applyTo: {div}.

remaining: remaining // div]. 

  remaining = 1] whileFalse:
    [div: next.
     next: next + 2] "Just look at the next odd integer."

].</lang>

Smalltalk

<lang 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 ]</lang>

Tcl

<lang tcl>namespace eval primes {}

proc primes::reset {} { 

   variable list [list]
   variable current_index end

}

namespace eval primes {reset}

proc primes::restart {} { 

   variable list
   variable current_index
   if {[llength $list] > 0} {
       set current_index 0
   }

}

proc primes::is_prime {candidate} { 

   variable list

   if {$candidate in $list} {return true}
   foreach prime $list {
       if {$candidate % $prime == 0} {
           return false
       }
       if {$prime * $prime > $candidate} {
           return true
       }
   }
   while true {
       set largest [get_next_prime]
       if {$largest * $largest >= $candidate} {
           return [is_prime $candidate]
       }
   }

}

proc primes::get_next_prime {} { 

   variable list
   variable current_index
   
   if {$current_index ne "end"} {
       set p [lindex $list $current_index]
       if {[incr current_index] == [llength $list]} {
           set current_index end
       }
       return $p
   }
   
   switch -exact -- [llength $list] {
       0 {set candidate 2}
       1 {set candidate 3}
       default {
           set candidate [lindex $list end]
           while true {
               incr candidate 2
               if {[is_prime $candidate]} break
           }
       }
   }
   lappend list $candidate
   return $candidate

}

  1. return the prime factors of a number in a dictionary.
  2. The keys will be the factors, the value will be the number
  3. of times the factor divides the given number
  5. example: 120 = 2**3 * 3 * 5, so
  6. [primes::factors 120] returns 2 3 3 1 5 1
  7. so: set prod 1
  8. dict for {p e} [primes::factors 120] {
  9. set prod [expr {$prod * $p**$e}]
  10. }
  11. expr {$prod == 120} ;# ==> true

proc primes::factors {num} { 

   restart
   set factors [dict create]
   for {set i [get_next_prime]} {$i <= $num} {} {
       if {$num % $i == 0} {
           dict incr factors $i
           set num [expr {$num / $i}]
           continue
       } elseif {$i*$i > $num} {
           dict incr factors $num
           break
       } else {
           set i [get_next_prime]
       }
   }
   return $factors

}</lang> Testing <lang tcl>primes::reset foreach m {2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59} { 

   set n [expr {2**$m - 1}]
   catch {time {set f [dict create {*}[primes::factors $n]]} 1} tm
   set primes [list]
   dict for {p e} $f {lappend primes {*}[lrepeat $e $p]}
   puts [format "2**%02d-1 = %-18s = %-22s => %s" $m $n [join $primes *] $tm]

}</lang> Outputs 

2**02-1 = 3                  = 3                      => 20 microseconds per iteration
2**03-1 = 7                  = 7                      => 16 microseconds per iteration
2**05-1 = 31                 = 31                     => 27 microseconds per iteration
2**07-1 = 127                = 127                    => 33 microseconds per iteration
2**11-1 = 2047               = 23*89                  => 43 microseconds per iteration
2**13-1 = 8191               = 8191                   => 159 microseconds per iteration
2**17-1 = 131071             = 131071                 => 535 microseconds per iteration
2**19-1 = 524287             = 524287                 => 911 microseconds per iteration
2**23-1 = 8388607            = 47*178481              => 162 microseconds per iteration
2**29-1 = 536870911          = 233*1103*2089          => 982 microseconds per iteration
2**31-1 = 2147483647         = 2147483647             => 138831 microseconds per iteration
2**37-1 = 137438953471       = 223*616318177          => 5154 microseconds per iteration
2**41-1 = 2199023255551      = 13367*164511353        => 2901 microseconds per iteration
2**43-1 = 8796093022207      = 431*9719*2099863       => 2141 microseconds per iteration
2**47-1 = 140737488355327    = 2351*4513*13264529     => 1102 microseconds per iteration
2**53-1 = 9007199254740991   = 6361*69431*20394401    => 97472 microseconds per iteration
2**59-1 = 576460752303423487 = 179951*3203431780337   => 12664437 microseconds per iteration

V

like in scheme (using variables) <lang v>[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].</lang>

(mostly) the same thing using stack (with out variables) <lang v>[prime-decomposition 

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

Using it <lang v>|1221 prime-decomposition puts</lang> 

=[3 11 37]
Retrieved from "https://rosettacode.org/wiki/Prime_decomposition?oldid=55743"