Prime decomposition: Difference between revisions

From Rosetta Code
Content added Content deleted
Line 4,926: Line 4,926:
=={{header|S-BASIC}}==
=={{header|S-BASIC}}==
<lang S-BASIC>
<lang S-BASIC>
rem - compute p mod q
rem - return p mod q
function mod(p, q = integer) = integer
function mod(p, q = integer) = integer
end = p - q * (p/q)
end = p - q * (p/q)
Line 4,935: Line 4,935:
Find the prime factors of n and store in global array factors
Find the prime factors of n and store in global array factors
(arrays cannot be passed as parameters) and return the number
(arrays cannot be passed as parameters) and return the number
found. If n is prime, it will be stored as the one and only
found. If n is prime, it will be stored as the only factor.
factor.
end
end
function primefactors(n = integer) = integer
function primefactors(n = integer) = integer

Revision as of 17:15, 8 August 2022

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}


Task

Write a function which returns an array or collection which contains the prime decomposition of a given number     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).


Related tasks



11l

Translation of: D

<lang 11l>F decompose(BigInt number)

  [BigInt] result
  V n = number
  BigInt i = 2
  L n % i == 0
     result.append(i)
     n I/= i
  i = 3
  L n >= i * i
     L n % i == 0
        result.append(i)
        n I/= i
     i += 2
  I n != 1
     result.append(n)
  R result

L(i) 2..9

  print(decompose(i))

print(decompose(1023 * 1024)) print(decompose(2 * 3 * 5 * 7 * 11 * 11 * 13 * 17)) print(decompose(BigInt(16860167264933) * 179951))</lang>

Output:
[2]
[3]
[2, 2]
[5]
[2, 3]
[7]
[2, 2, 2]
[3, 3]
[2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 11, 31]
[2, 3, 5, 7, 11, 11, 13, 17]
[179951, 16860167264933]

360 Assembly

For maximum compatibility, this program uses only the basic instruction set. <lang 360asm>PRIMEDE CSECT

        USING  PRIMEDE,R13
        B      80(R15)            skip savearea
        DC     17F'0'             savearea
        DC     CL8'PRIMEDE'
        STM    R14,R12,12(R13)
        ST     R13,4(R15)
        ST     R15,8(R13)
        LR     R13,R15            end prolog
        LA     R2,0
        LA     R3,1023
        LA     R4,1024
        MR     R2,R4
        ST     R3,N               n=1023*1024
        LA     R5,WBUFFER
        LA     R6,0
        L      R1,N               n
        XDECO  R1,0(R5)
        LA     R5,12(R5)
        MVC    0(3,R5),=C' : '
        LA     R5,3(R5)
        LA     R0,2
        ST     R0,I               i=2

WHILE1 EQU * do while(i<=n/2)

        L      R2,N
        SRA    R2,1
        L      R4,I
        CR     R4,R2              i<=n/2
        BH     EWHILE1

WHILE2 EQU * do while(n//i=0)

        L      R3,N
        LA     R2,0
        D      R2,I
        LTR    R2,R2              n//i=0
        BNZ    EWHILE2
        ST     R3,N               n=n/i
        ST     R3,M               m=n
        L      R1,I               i
        XDECO  R1,WDECO
        MVC    0(5,R5),WDECO+7
        LA     R5,5(R5)
        MVI    OK,X'01'           ok
        B      WHILE2

EWHILE2 EQU *

        L      R4,I
        CH     R4,=H'2'           if i=2 then
        BNE    NE2
        LA     R0,3
        ST     R0,I               i=3
        B      EIFNE2

NE2 L R2,I else

        LA     R2,2(R2)
        ST     R2,I               i=i+2

EIFNE2 B WHILE1 EWHILE1 EQU *

        CLI    OK,X'01'           if ^ok then
        BE     NOTPRIME
        MVC    0(7,R5),=C'[prime]'
        LA     R5,7(R5)
        B      EPRIME

NOTPRIME L R1,M m

        XDECO  R1,WDECO
        MVC    0(5,R5),WDECO+7

EPRIME XPRNT WBUFFER,80 put

        L      R13,4(0,R13)       epilog
        LM     R14,R12,12(R13)
        XR     R15,R15
        BR     R14

N DS F I DS F M DS F OK DC X'00' WBUFFER DC CL80' ' WDECO DS CL16

        YREGS  
        END    PRIMEDE</lang>
Output:
     1047552 :     2    2    2    2    2    2    2    2    2    2    3   11   31

AArch64 Assembly

Works with: as version Raspberry Pi 3B version Buster 64 bits

<lang AArch64 Assembly> /* ARM assembly AARCH64 Raspberry PI 3B */ /* program primeDecomp64.s */

/*******************************************/ /* Constantes file */ /*******************************************/ /* for this file see task include a file in language AArch64 assembly*/ .include "../includeConstantesARM64.inc" .equ NBFACT, 100

/*******************************************/ /* Structures */ /********************************************/ /* structurea area factors */

   .struct  0

fac_value: // factor

   .struct  fac_value + 8

fac_number: // number of identical factors

   .struct  fac_number + 8

fac_end: /*******************************************/ /* Initialized data */ /*******************************************/ .data szMessStartPgm: .asciz "Program start \n" szMessEndPgm: .asciz "Program normal end.\n" szMessNotPrime: .asciz "Not prime.\n" szMessPrime: .asciz "Prime\n" szCarriageReturn: .asciz "\n" szSpaces: .asciz " " szMessNumber: .asciz " The factors of @ are :\n" /*******************************************/ /* UnInitialized data */ /*******************************************/ .bss sZoneConv: .skip 32 .align 4 tbZoneDecom: .skip fac_end * NBFACT /*******************************************/ /* code section */ /*******************************************/ .text .global main main: // program start

   ldr x0,qAdrszMessStartPgm            // display start message
   bl affichageMess
   ldr x20,qVal
   //mov x20,17
   mov x0,x20
   ldr x1,qAdrtbZoneDecom
   bl decompFact                        // decomposition
   cmp x0,#0
   beq 1f
   mov x2,x0
   mov x0,x20
   ldr x1,qAdrtbZoneDecom
   bl displayFactors                    // display factors
   b 2f

1:

   ldr x0,qAdrszMessPrime              // prime
   bl affichageMess
   

2:

   ldr x0,qAdrszMessEndPgm            // display end message
   bl affichageMess

100: // standard end of the program

   mov x0,0                           // return code
   mov x8,EXIT                        // request to exit program
   svc 0                              // perform system call

qAdrszMessStartPgm: .quad szMessStartPgm qAdrszMessEndPgm: .quad szMessEndPgm qAdrszCarriageReturn: .quad szCarriageReturn qAdrszMessNotPrime: .quad szMessNotPrime qAdrszMessPrime: .quad szMessPrime qAdrtbZoneDecom: .quad tbZoneDecom //qVal: .quad 2 <<31 qVal: .quad 1047552 // test not prime //qVal: .quad 1429671721 // test not prime (37811 * 37811)

/******************************************************************/ /* prime decomposition */ /******************************************************************/ /* x0 contains the number */ /* x1 contains address factors array */ /* REMARK no save register x9-x19 */ decompFact:

   stp x1,lr,[sp,-16]!       // save  registers
   mov x12,x0                // save number
   bl isPrime                // prime ?
   cbnz x0,12f               // yes -> no decomposition 
   mov x19,fac_end           // element area size
   mov x18,0                 // raz indice
   mov x16,0                 // prev divisor
   mov x17,0                 // number of identical divisors
   mov x13,2                 // first divisor

2:

   cmp x12,1
   beq 10f
   udiv x14,x12,x13          // division
   msub x15,x14,x13,x12      // remainder = x12 -(x13*x14)
   cbnz x15,5f               // if remainder <> zero  x13 not divisor
   mov x12,x14               // quotient -> new dividende
   cmp x13,x16               // same divisor ?
   beq 4f                    // yes
   cbz x16,3f                // yes it is first divisor ?
   madd x11,x18,x19,x1       // no -> store prev divisor in the area
   str x16,[x11,fac_value]
   str x17,[x11,fac_number]  // and store number of same factor
   add x18,x18,1             // increment indice
   mov x17,0                 // raz number of same factor

3:

   mov x16,x13               // save new divisor

4:

   add x17,x17,1             // increment number of same factor
   mov x0,x12                // the new dividende is prime ?
   bl isPrime
   cbnz x0,10f               // yes
   b 2b                      // else loop

5: // divisor is not a factor

   cmp x13,2                 // begin ?
   cinc x13,x13,ne           // if divisor <> 2 add 1
   add x13,x13,1
   b 2b                      // and loop

10: // new dividende is prime

   cmp x16,x12               // divisor = dividende ?
   cinc x17,x17,eq           //add 1 if last dividende = diviseur
   madd x11,x18,x19,x1
   str x16,[x11,fac_value]   // store divisor in area
   str x17,[x11,fac_number]  // and store number
   add x18,x18,1             // increment indice
   cmp x16,x12               //store last dividende if <> diviseur
   beq 11f
   madd x11,x18,x19,x1
   str x12,[x11,fac_value]   // sinon stockage dans la table
   mov x17,1
   str x17,[x11,fac_number]  // store 1 in number
   add x18,x18,1

11:

   mov x0,x18                // return nb factors
   b 100f

12:

   mov x0,#0                 // number is prime 
   b 100f

100:

   ldp x1,lr,[sp],16         // restaur des  2 registres
   ret                       // retour adresse lr x30

/******************************************************************/ /* prime decomposition */ /******************************************************************/ /* x0 contains the number */ /* x1 contains address factors array */ /* x2 number of factors */ displayFactors:

   stp x1,lr,[sp,-16]!        // save  registres
   mov x19,fac_end            // element area size
   mov x13,x1                 // save area address
   ldr x1,qAdrsZoneConv       // load zone conversion address
   bl conversion10
   ldr x0,qAdrszMessNumber
   bl strInsertAtCharInc                   // insert result at Second @ character
   bl affichageMess
   mov x9,0                   // indice

1:

   madd x10,x9,x19,x13        // compute address area element
   ldr x0,[x10,fac_value]
   ldr x12,[x10,fac_number]
   bl conversion10            // decimal conversion

2:

   mov x0,x1
   bl affichageMess
   ldr x0,qAdrszSpaces
   bl affichageMess
   subs x12,x12,#1
   bgt 2b
   add x9,x9,1
   cmp x9,x2
   blt 1b
   ldr x0,qAdrszCarriageReturn
   bl affichageMess

100:

   ldp x1,lr,[sp],16          // restaur des  2 registres
   ret                        // retour adresse lr x30

qAdrsZoneConv: .quad sZoneConv qAdrszSpaces: .quad szSpaces qAdrszMessNumber: .quad szMessNumber /******************************************************************/ /* test if number is prime */ /******************************************************************/ /* x0 contains the number */ /* x0 return 1 if prime else return 0 */ isPrime:

   stp x1,lr,[sp,-16]!       // save  registers
   cmp x0,1                  // <= 1 ?
   ble 98f
   cmp x0,3                  // 2 and 3 prime
   ble 97f 
   tst x0,1                  //  even ?
   beq 98f 
   mov x9,3                  // first divisor

1:

   udiv x11,x0,x9
   msub x10,x11,x9,x0        // compute remainder
   cbz x10,98f               // end if zero
   add x9,x9,#2              // increment divisor
   cmp x9,x11                // divisors<=quotient ?
   ble 1b                    // loop

97:

   mov x0,1                  // return prime
   b 100f

98:

   mov x0,0                  // not prime
   b 100f

100:

   ldp x1,lr,[sp],16         // restaur  2 registers
   ret                       // return to address lr x30

/********************************************************/ /* File Include fonctions */ /********************************************************/ /* for this file see task include a file in language AArch64 assembly */ .include "../includeARM64.inc"

</lang>

Output:
Program start
 The factors of 1047552 are :
2  2  2  2  2  2  2  2  2  2  3  11  31
Program normal end.

ABAP

<lang ABAP>class ZMLA_ROSETTA definition

 public
 create public .
 public section.
   types:
     enumber         TYPE          N  LENGTH 60,
     listof_enumber  TYPE TABLE OF enumber .
   class-methods FACTORS
     importing
       value(N) type ENUMBER
     exporting
       value(ORET) type LISTOF_ENUMBER .
 protected section.
 private section.

ENDCLASS.


CLASS ZMLA_ROSETTA IMPLEMENTATION.


  • <SIGNATURE>---------------------------------------------------------------------------------------+
  • | Static Public Method ZMLA_ROSETTA=>FACTORS
  • +-------------------------------------------------------------------------------------------------+
  • | [--->] N TYPE ENUMBER
  • | [<---] ORET TYPE LISTOF_ENUMBER
  • +--------------------------------------------------------------------------------------</SIGNATURE>
 method FACTORS.
   CLEAR oret.
   WHILE n mod 2 = 0.
     n = n / 2.
     APPEND 2 to oret.
   ENDWHILE.
   DATA: lim type enumber,
         i   type enumber.
   lim = sqrt( n ).
   i   = 3.
   WHILE i <= lim.
     WHILE n mod i = 0.
       APPEND i to oret.
       n = n / i.
       lim = sqrt( n ).
     ENDWHILE.
     i = i + 2.
   ENDWHILE.
   IF n > 1.
     APPEND n to oret.
   ENDIF.
 endmethod.

ENDCLASS.</lang>

ACL2

<lang Lisp>(include-book "arithmetic-3/top" :dir :system)

(defun prime-factors-r (n i)

  (declare (xargs :mode :program))
  (cond ((or (zp n) (zp (- n i)) (zp i) (< i 2) (< n 2))
         (list n))
        ((= (mod n i) 0)
         (cons i (prime-factors-r (floor n i) 2)))
        (t (prime-factors-r n (1+ i)))))

(defun prime-factors (n)

  (declare (xargs :mode :program))
  (prime-factors-r n 2))</lang>

Ada

The solution is generic.

The package Prime_Numbers 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 same package is used for Almost prime#Ada, Semiprime#Ada, Count in factors#Ada, Primality by Trial Division#Ada, Sequence of primes by Trial Division#Ada, and Ulam_spiral_(for_primes)#Ada.

This is the specification of the generic package Prime_Numbers:

<lang ada>generic

  type Number is private;
  Zero : Number;
  One  : Number;
  Two  : Number;
  with function "+"   (X, Y : Number) return Number 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;
  function Is_Prime (N : Number) return Boolean;

end Prime_Numbers;</lang>

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.

This is the implementation of the generic package Prime_Numbers:

<lang ada>package body Prime_Numbers is

-- auxiliary (internal) functions
  function First_Factor (N : Number; Start : Number) return Number is
     K    : Number  := Start;
  begin
     while ((N mod K) /= Zero) and then (N > (K*K))  loop 
        K := K + One;
     end loop;
     if (N mod K) = Zero then
        return K;
     else
        return N;
     end if;
  end First_Factor;
  
  function Decompose (N : Number; Start : Number) return Number_List is
     F: Number := First_Factor(N, Start);
     M: Number := N / F;
  begin
     if M = One then -- F is the last factor
        return (1 => F);
     else
        return F & Decompose(M, Start);
     end if;
  end Decompose;
  
-- functions visible from the outside
  function Decompose (N : Number) return Number_List is (Decompose(N, Two));
  function Is_Prime (N : Number) return Boolean is
     (N > One and then First_Factor(N, Two)=N);

end Prime_Numbers;</lang>

In the example provided, the package Prime_Numbers is instantiated with plain integer type:

<lang ada>with Prime_Numbers, Ada.Text_IO;

procedure Test_Prime is

  package Integer_Numbers is new 
    Prime_Numbers (Natural, 0, 1, 2); 
  use Integer_Numbers;
    
  procedure Put (List : Number_List) is
  begin
     for Index in List'Range loop
        Ada.Text_IO.Put (Positive'Image (List (Index)));
     end loop;
  end Put;
    

begin

  Put (Decompose (12));

end Test_Prime;</lang>

Output:

(decomposition of 12)

 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 gen 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 # gen 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 gen 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 # gen decompose := (LINT in n, YIELDLINT yield)VOID: (
   LINT n := in n;
 # FOR LINT p IN # gen 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 # gen 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 # gen 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.

ALGOL-M

Sadly, ALGOL-M does not allow arrays to be passed as parameters to procedures or functions, so the routine must store its results in (and know the name of) the external array used for that purpose. <lang ALGOL> BEGIN

INTEGER I, K, NFOUND; INTEGER ARRAY FACTORS[1:16];

COMMENT - RETURN P MOD Q; INTEGER FUNCTION MOD (P, Q); INTEGER P, Q; BEGIN

   MOD := P - Q * (P / Q);

END;

COMMENT

 FIND THE PRIME FACTORS OF N AND STORE IN THE EXTERNAL
 ARRAY "FACTORS", RETURNING THE NUMBER FOUND. IF N IS
 PRIME, IT WILL BE STORED AS THE FIRST AND ONLY FACTOR;

INTEGER FUNCTION PRIMEFACTORS(N); INTEGER N; BEGIN

 INTEGER P, COUNT;
 P := 2;
 COUNT := 1;
 WHILE N >= P * P DO
   BEGIN
     IF MOD(N, P) = 0 THEN
       BEGIN
         FACTORS[COUNT] := P;
         COUNT := COUNT + 1;
         N := N / P;
       END
     ELSE
       P := P + 1;
   END;
 FACTORS[COUNT] := N;
 PRIMEFACTORS := COUNT;

END;

COMMENT -- EXERCISE THE ROUTINE;

FOR I := 77 STEP 2 UNTIL 99 DO

 BEGIN
   WRITE(I,":");
   NFOUND := PRIMEFACTORS(I);
   COMMENT - PRINT OUT THE FACTORS THAT WERE FOUND;
   FOR K := 1 STEP 1 UNTIL NFOUND DO
     BEGIN
       WRITEON(FACTORS[K]);
     END;
 END;

END </lang>

Output:
   77:     7    11
   79:    79
   81:     3     3     3     3
   83:    83
   85:     5    17
   87:     3    29
   89:    89
   91:     7    13
   93:     3    31
   95:     5    19
   97:    97
   99:     3     3    11

Applesoft BASIC

<lang ApplesoftBasic>9040 PF(0) = 0 : SC = 0 9050 FOR CA = 2 TO INT( SQR(I)) 9060 IF I = 1 THEN RETURN 9070 IF INT(I / CA) * CA = I THEN GOSUB 9200 : GOTO 9060 9080 CA = CA + SC : SC = 1 9090 NEXT CA 9100 IF I = 1 THEN RETURN 9110 CA = I

9200 PF(0) = PF(0) + 1 9210 PF(PF(0)) = CA 9220 I = I / CA 9230 RETURN</lang>

Arturo

<lang rebol>decompose: function [num][

   facts: to [:string] factors.prime num
   print [
       pad.right (to :string num) ++ " = " ++  join.with:" x " facts 30 
       "{"++ (join.with:", " unique facts) ++ "}"
   ]

]

loop 2..40 => decompose</lang>

Output:
2 = 2                          {2} 
3 = 3                          {3} 
4 = 2 x 2                      {2} 
5 = 5                          {5} 
6 = 2 x 3                      {2, 3} 
7 = 7                          {7} 
8 = 2 x 2 x 2                  {2} 
9 = 3 x 3                      {3} 
10 = 2 x 5                     {2, 5} 
11 = 11                        {11} 
12 = 2 x 2 x 3                 {2, 3} 
13 = 13                        {13} 
14 = 2 x 7                     {2, 7} 
15 = 3 x 5                     {3, 5} 
16 = 2 x 2 x 2 x 2             {2} 
17 = 17                        {17} 
18 = 2 x 3 x 3                 {2, 3} 
19 = 19                        {19} 
20 = 2 x 2 x 5                 {2, 5} 
21 = 3 x 7                     {3, 7} 
22 = 2 x 11                    {2, 11} 
23 = 23                        {23} 
24 = 2 x 2 x 2 x 3             {2, 3} 
25 = 5 x 5                     {5} 
26 = 2 x 13                    {2, 13} 
27 = 3 x 3 x 3                 {3} 
28 = 2 x 2 x 7                 {2, 7} 
29 = 29                        {29} 
30 = 2 x 3 x 5                 {2, 3, 5} 
31 = 31                        {31} 
32 = 2 x 2 x 2 x 2 x 2         {2} 
33 = 3 x 11                    {3, 11} 
34 = 2 x 17                    {2, 17} 
35 = 5 x 7                     {5, 7} 
36 = 2 x 2 x 3 x 3             {2, 3} 
37 = 37                        {37} 
38 = 2 x 19                    {2, 19} 
39 = 3 x 13                    {3, 13} 
40 = 2 x 2 x 2 x 5             {2, 5}

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)
           return, % f "`n" next
       }
       f++
   }

}</lang>

Optimized Version

Translation of: JavaScript

<lang AutoHotkey>prime_numbers(n) {

   if (n <= 3)
       return [n]
   ans := []
   done := false
   while !done 
   {
       if !Mod(n,2){
           ans.push(2)
           n /= 2
           continue
       }
       if !Mod(n,3) {
           ans.push(3)
           n /= 3
           continue
       }
       if (n = 1)
           return ans
       
       sr := sqrt(n)
       done := true
       ; try to divide the checked number by all numbers till its square root.
       i := 6
       while (i <= sr+6){
           if !Mod(n, i-1) { ; is n divisible by i-1?
               ans.push(i-1)
               n /= i-1
               done := false
               break
           }
           if !Mod(n, i+1) { ; is n divisible by i+1?
               ans.push(i+1)
               n /= i+1
               done := false
               break
           }
           i += 6
       }
   }
   ans.push(n)
   return ans

}</lang> Examples:<lang AutoHotkey>num := 8388607, output := "" for i, p in prime_numbers(num)

   output .= p " * "

MsgBox % num " = " Trim(output, " * ") return</lang>

Output:
8388607 = 47 * 178481


AWK

As the examples show, pretty large numbers can be factored in tolerable time:

<lang awk># Usage: awk -f primefac.awk function pfac(n, r, f){

   r = ""; f = 2
   while (f <= n) {
       while(!(n % f)) {
           n = n / f
           r = r " " f
       }
       f = f + 2 - (f == 2)
   }
   return r

}

  1. For each line of input, print the prime factors.

{ print pfac($1) } </lang>

Output:

entering input on stdin

$
36
 2 2 3 3
77
 7 11
536870911
 233 1103 2089
8796093022207
 431 9719 2099863

Batch file

Unfortunately Batch does'nt have a BigNum library so the maximum number that can be decomposed is 2^31-1 <lang Batch file> @echo off

usage: cmd /k primefactor.cmd number

setlocal enabledelayedexpansion

set /a compo=%1 if "%compo%"=="" goto:eof set list=%compo%= (

set /a div=2 & call :loopdiv set /a div=3 & call :loopdiv set /a div=5,inc=2

looptest

call :loopdiv set /a div+=inc,inc=6-inc,div2=div*div if %div2% lss %compo% goto looptest if %compo% neq 1 set list= %list% %compo% echo %list%) & goto:eof

loopdiv

set /a "res=compo%%div if %res% neq 0 goto:eof set list=%list% %div%, set/a compo/=div goto:loopdiv </lang>

Befunge

Works with: befungee

Handles safely integers only up to 250 (or ones which don't have prime divisors greater than 250). <lang Befunge>& 211p > : 1 - #v_ 25*, @ > 11g:. / v

               > : 11g %!|
                         > 11g 1+ 11p v
      ^                               <</lang>

BQN

An efficient Factor function using trial division and Pollard's rho algorithm is given in bqn-libs primes.bqn. The following standalone version is based on the trial division there, and builds in the sieve from Extensible prime generator. <lang bqn>Factor ← { 𝕊n:

 # Prime sieve
 primes ← ↕0
 Sieve ← { p 𝕊 a‿b:
   p(⍋↑⊣)↩√b ⋄ l←b-a
   E ← {↕∘⌈⌾(((𝕩|-a)+𝕩×⊢)⁼)l}        # Indices of multiples of 𝕩
   a + / (1⥊˜l) E⊸{0¨⌾(𝕨⊸⊏)𝕩}´ p     # Primes in segment [a,b)
 }
 # Factor by trial division
 r ← ↕0                              # Result list
 Try ← {
   m ← (1+⌊√n) ⌊ 2×𝕩                 # Upper bound for factors this round
   𝕩<m ?                             # Stop if no factors
   primes ∾↩ np ← primes Sieve 𝕩‿m   # New primes
   {0=𝕩|n? r∾↩𝕩 ⋄ n÷↩𝕩 ⋄ 𝕊𝕩 ;@}¨ np  # Try each one
   𝕊 m                               # Next segment
 ;@}
 Try 2
 r ∾ 1⊸<⊸⥊n

}</lang>

Output:

<lang bqn> > ⋈⟜Factor¨ 1232123+↕4 # Some factored numbers ┌─ ╵ 1232123 ⟨ 29 42487 ⟩

 1232124 ⟨ 2 2 3 102677 ⟩   
 1232125 ⟨ 5 5 5 9857 ⟩     
 1232126 ⟨ 2 7 17 31 167 ⟩  
                           ┘</lang>

Burlesque

<lang burlesque> blsq ) 12fC {2 2 3} </lang>

C

Version 1

Relatively sophiscated sieve method based on size 30 prime wheel. The code does not pretend to handle prime factors larger than 64 bits. All 32-bit primes are cached with 137MB data. Cache data takes about a minute to compute the first time the program is run, which is also saved to the current directory, and will be loaded in a second if needed again. <lang c>#include <inttypes.h>

  1. include <stdio.h>
  2. include <stdlib.h>
  3. include <string.h>
  4. include <assert.h>

typedef uint32_t pint; typedef uint64_t xint; typedef unsigned int uint;

  1. define PRIuPINT PRIu32 /* printf macro for pint */
  2. define PRIuXINT PRIu64 /* printf macro for xint */
  3. define MAX_FACTORS 63 /* because 2^64 is too large for xint */

uint8_t *pbits;

  1. define MAX_PRIME (~(pint)0)
  2. define MAX_PRIME_SQ 65535U
  3. define PBITS (MAX_PRIME / 30 + 1)

pint next_prime(pint); int is_prime(xint); void sieve(pint);

uint8_t bit_pos[30] = {

   0, 1<<0, 0, 0, 0,    0,
   0, 1<<1, 0, 0, 0, 1<<2,
   0, 1<<3, 0, 0, 0, 1<<4,
   0, 1<<5, 0, 0, 0, 1<<6,
   0,    0, 0, 0, 0, 1<<7,

};

uint8_t rem_num[] = { 1, 7, 11, 13, 17, 19, 23, 29 };

void init_primes() {

   FILE *fp;
   pint s, tgt = 4;
   if (!(pbits = malloc(PBITS))) {
       perror("malloc");
       exit(1);
   }
   if ((fp = fopen("primebits", "r"))) {
       fread(pbits, 1, PBITS, fp);
       fclose(fp);
       return;
   }
   memset(pbits, 255, PBITS);
   for (s = 7; s <= MAX_PRIME_SQ; s = next_prime(s)) {
       if (s > tgt) {
           tgt *= 2;
           fprintf(stderr, "sieve %"PRIuPINT"\n", s);
       }
       sieve(s);
   }
   fp = fopen("primebits", "w");
   fwrite(pbits, 1, PBITS, fp);
   fclose(fp);

}

int is_prime(xint x) {

   pint p;
   if (x > 5) {
       if (x < MAX_PRIME)
           return pbits[x/30] & bit_pos[x % 30];
       for (p = 2; p && (xint)p * p <= x; p = next_prime(p))
           if (x % p == 0) return 0;
       return 1;
   }
   return x == 2 || x == 3 || x == 5;

}

void sieve(pint p) {

   unsigned char b[8];
   off_t ofs[8];
   int i, q;
   for (i = 0; i < 8; i++) {
       q = rem_num[i] * p;
       b[i] = ~bit_pos[q % 30];
       ofs[i] = q / 30;
   }
   for (q = ofs[1], i = 7; i; i--)
       ofs[i] -= ofs[i-1];
   for (ofs[0] = p, i = 1; i < 8; i++)
       ofs[0] -= ofs[i];
   for (i = 1; q < PBITS; q += ofs[i = (i + 1) & 7])
       pbits[q] &= b[i];

}

pint next_prime(pint p) {

   off_t addr;
   uint8_t bits, rem;
   if (p > 5) {
       addr = p / 30;
       bits = bit_pos[ p % 30 ] << 1;
       for (rem = 0; (1 << rem) < bits; rem++);
       while (pbits[addr] < bits || !bits) {
           if (++addr >= PBITS) return 0;
           bits = 1;
           rem = 0;
       }
       if (addr >= PBITS) return 0;
       while (!(pbits[addr] & bits)) {
           rem++;
           bits <<= 1;
       }
       return p = addr * 30 + rem_num[rem];
   }
   switch(p) {
       case 2: return 3;
       case 3: return 5;
       case 5: return 7;
   }
   return 2;

}

int decompose(xint n, xint *f) {

   pint p = 0;
   int i = 0;
   /* check small primes: not strictly necessary */
   if (n <= MAX_PRIME && is_prime(n)) {
       f[0] = n;
       return 1;
   }
   while (n >= (xint)p * p) {
       if (!(p = next_prime(p))) break;
       while (n % p == 0) {
           n /= p;
           f[i++] = p;
       }
   }
   if (n > 1) f[i++] = n;
   return i;

}

int main() {

   int i, len;
   pint p = 0;
   xint f[MAX_FACTORS], po;
   init_primes();
   for (p = 1; p < 64; p++) {
       po = (1LLU << p) - 1;
       printf("2^%"PRIuPINT" - 1 = %"PRIuXINT, p, po);
       fflush(stdout);
       if ((len = decompose(po, f)) > 1)
           for (i = 0; i < len; i++)
               printf(" %c %"PRIuXINT, i?'x':'=', f[i]);
       putchar('\n');
   }
   return 0;

}</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 */

/* 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.

To understand what was going on with the above code, pass it through cpp and read the outcome. Translated into normal C code sans the function call overhead, it's really this (the following uses a adjustable cache, although setting it beyond a few thousands doesn't gain further benefit):<lang c>#include <stdio.h>

  1. include <stdlib.h>
  2. include <stdint.h>

typedef uint32_t pint; typedef uint64_t xint; typedef unsigned int uint;

int is_prime(xint);

inline int next_prime(pint p) {

   if (p == 2) return 3;
   for (p += 2; p > 1 && !is_prime(p); p += 2);
   if (p == 1) return 0;
   return p;

}

int is_prime(xint n) {

  1. define NCACHE 256
  2. define S (sizeof(uint) * 2)
   static uint cache[NCACHE] = {0};

   pint p = 2;
   int ofs, bit = -1;

   if (n < NCACHE * S) {
       ofs = n / S;
       bit = 1 << ((n & (S - 1)) >> 1);
       if (cache[ofs] & bit) return 1;
   }

   do {
       if (n % p == 0) return 0;
       if (p * p > n) break;
   } while ((p = next_prime(p)));

   if (bit != -1) cache[ofs] |= bit;
   return 1;

}

int decompose(xint n, pint *out) {

   int i = 0;
   pint p = 2;
   while (n > p * p) {
       while (n % p == 0) {
           out[i++] = p;
           n /= p;
       }
       if (!(p = next_prime(p))) break;
   }
   if (n > 1) out[i++] = n;
   return i;

}

int main() {

   int i, j, len;
   xint z;
   pint out[100];
   for (i = 2; i < 64; i = next_prime(i)) {
       z = (1ULL << i) - 1;
       printf("2^%d - 1 = %llu = ", i, z);
       fflush(stdout);
       len = decompose(z, out);
       for (j = 0; j < len; j++)
           printf("%u%s", out[j], j < len - 1 ? " x " : "\n");
   }

   return 0;

}</lang>


Version 2

<lang c> typedef unsigned long long int ulong; // define a type that represent the limit (64-bit)

ulong mod_mul(ulong a, ulong b, const ulong mod) { ulong res = 0, c; // return (a * b) % mod, avoiding overflow errors while doing modular multiplication. for (b %= mod; a; a & 1 ? b >= mod - res ? res -= mod : 0, res += b : 0, a >>= 1, (c = b) >= mod - b ? c -= mod : 0, b += c); return res % mod; }

ulong mod_pow(ulong n, ulong exp, const ulong mod) { ulong res = 1; // return (n ^ exp) % mod for (n %= mod; exp; exp & 1 ? res = mod_mul(res, n, mod) : 0, n = mod_mul(n, n, mod), exp >>= 1); return res; }

ulong square_root(const ulong N) { ulong res = 0, rem = N, c, d; for (c = 1 << 62; c; c >>= 2) { d = res + c; res >>= 1; if (rem >= d) rem -= d, res += c; } // returns the square root of N. return res; }

int is_prime(const ulong N) { ulong i = 1; // return a truthy value about the primality of N. if (N > 1) for (; i < 64 && mod_pow(i, N - 1, N) <= 1; ++i); return i == 64; }

ulong pollard_rho(const ulong N) { // Require : N is a composite number, not a square. // Ensure : res is a non-trivial factor of N. // Option : change the timeout, change the rand function. static const int timeout = 18; static unsigned long long rand_val = 2994439072U; rand_val = (rand_val * 1025416097U + 286824428U) % 4294967291LLU; ulong res = 1, a, b, c, i = 0, j = 1, x = 1, y = 1 + rand_val % (N - 1); for (; res == 1; ++i) { if (i == j) { if (j >> timeout) break; j <<= 1; x = y; } a = y, b = y; // performs y = (y * y) % N for (y = 0; a; a & 1 ? b >= N - y ? y -= N : 0, y += b : 0, a >>= 1, (c = b) >= N - b ? c -= N : 0, b += c); y = (1 + y) % N; for (a = y > x ? y - x : x - y, b = N; (a %= b) && (b %= a);); // compute the gcd(abs(y - x), N); res = a | b; } return res; }

void factor(const ulong N, ulong *array) { // very basic manager that fill the given array (the size of the result is the first array element) // it does not perform initial trial divisions, which is generally highly recommended. if (N < 4 || is_prime(N)) { if (N > 1 || !*array) array[++*array] = N; return; } ulong x = square_root(N); if (x * x != N) x = pollard_rho(N); factor(x, array); factor(N / x, array); }

  1. include <stdio.h>

int main(void) { // simple test. unsigned long long n = 18446744073709551615U; ulong fac[65] = {0}; factor(n, fac); for (ulong i = 1; i <= *fac; ++i) printf("* %llu\n", fac[i]); }

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

Simple trial division

This version a translation from Java of the sample presented by Robert C. Martin during a TDD talk at NDC 2011. Although this three-line algorithm does not mention anything about primes, the fact that factors are taken out of the number n in ascending order garantees the list will only contain primes.

<lang csharp>using System.Collections.Generic;

namespace PrimeDecomposition {

   public class Primes
   {
       public List<int> FactorsOf(int n)
       {
           var factors = new List<int>();
           for (var divisor = 2; n > 1; divisor++)
               for (; n % divisor == 0; n /= divisor)
                   factors.Add(divisor);
           return factors;
       }

}</lang>

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>

Simple trial division

<lang cpp>// Factorization by trial division in C++11

  1. include <iostream>
  2. include <vector>

using long_pair = std::pair<long,long>; using lp_vec = std::vector<long_pair>;

lp_vec factorize(long n) {

   lp_vec fs;
   int cnt = 0;
   for (;n%2==0; n/=2) cnt++;    // optimized by compiler
   if (cnt > 0)
       fs.push_back({2, cnt});
   for (long i=3; i*i<=n; i+=2) {
       cnt = 0;
       for (;n%i==0; n/=i) cnt++;
       if (cnt>0)
           fs.push_back({i, cnt});
   }
   if (n>1)
       fs.push_back({n, 1});
   return fs;

}

int main() {

   long n;
   std::cin >> n;
   auto fs = factorize(n);
   for (auto fp : fs) {
       std::cout << fp.first << "^" << fp.second << "\n";
   }
   return 0;

}</lang>

Clojure

<lang clojure>;;; No stack consuming algorithm (defn factors

 "Return a list of factors of N."
 ([n]
   (factors n 2 ()))
 ([n k acc]
   (if (= 1 n)      
     acc
     (if (= 0 (rem n k))        
       (recur (quot n k) k (cons k acc))
       (recur n (inc k) acc)))))</lang>

Commodore BASIC

Works with: Commodore BASIC version 2.0

It's not easily possible to have arbitrary precision integers in PET basic, so here is at least a version using built-in data types (reals). On return from the subroutine starting at 9000 the global array pf contains the number of factors followed by the factors themselves: <lang zxbasic>9000 REM ----- function generate 9010 REM in ... i ... number 9020 REM out ... pf() ... factors 9030 REM mod ... ca ... pf candidate 9040 pf(0)=0 : ca=2 : REM special case 9050 IF i=1 THEN RETURN 9060 IF INT(i/ca)*ca=i THEN GOSUB 9200 : GOTO 9050 9070 FOR ca=3 TO INT( SQR(i)) STEP 2 9080 IF i=1 THEN RETURN 9090 IF INT(i/ca)*ca=i THEN GOSUB 9200 : GOTO 9080 9100 NEXT 9110 IF i>1 THEN ca=i : GOSUB 9200 9120 RETURN 9200 pf(0)=pf(0)+1 9210 pf(pf(0))=ca 9220 i=i/ca 9230 RETURN</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>

<lang Lisp>;;; Tail-recursive version (defun factor (n &optional (acc '()))

 (when (> n 1) (loop with max-d = (isqrt n)
          for d = 2 then (if (evenp d) (1+ d) (+ d 2)) do
            (cond ((> d max-d) (return (cons (list n 1) acc)))
              ((zerop (rem n d)) 
               (return (factor (truncate n d) (if (eq d (caar acc))
                                  (cons 
                               (list (caar acc) (1+ (cadar acc)))
                               (cdr acc))
                                  (cons (list d 1) acc)))))))))</lang>

D

<lang d>import std.stdio, std.bigint, std.algorithm, std.traits, std.range;

Unqual!T[] decompose(T)(in T number) pure nothrow in {

   assert(number > 1);

} body {

   typeof(return) result;
   Unqual!T n = number;
   for (Unqual!T i = 2; n % i == 0; n /= i)
       result ~= i;
   for (Unqual!T i = 3; n >= i * i; i += 2)
       for (; n % i == 0; n /= i)
           result ~= i;
   if (n != 1)
       result ~= n;
   return result;

}

void main() {

   writefln("%(%s\n%)", iota(2, 10).map!decompose);
   decompose(1023 * 1024).writeln;
   BigInt(2 * 3 * 5 * 7 * 11 * 11 * 13 * 17).decompose.writeln;
   decompose(16860167264933UL.BigInt * 179951).writeln;
   decompose(2.BigInt ^^ 100_000).group.writeln;

}</lang>

Output:
[2]
[3]
[2, 2]
[5]
[2, 3]
[7]
[2, 2, 2]
[3, 3]
[2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 11, 31]
[2, 3, 5, 7, 11, 11, 13, 17]
[179951, 16860167264933]
[Tuple!(BigInt, uint)(2, 100000)]

Delphi

Translation of: C#

<lang Delphi> program Prime_decomposition;

{$APPTYPE CONSOLE}

uses

 System.SysUtils;

function IsPrime(n: UInt64): Boolean; var

 i: Integer;

begin

 if n <= 1 then
   exit(False);
 i := 2;
 while i < Sqrt(n) do
 begin
   if n mod i = 0 then
     exit(False);
   inc(i);
 end;
 Result := True;

end;

function GetPrimes(n: UInt64): TArray<UInt64>; var

 i: Integer;

begin

 while n > 1 do
 begin
   i := 1;
   while True do
   begin
     if IsPrime(i) then
     begin
       if n / i = (round(n / i)) then
       begin
         n := n div i;
         SetLength(Result, Length(Result) + 1);
         Result[High(Result)] := i;
         Break;
       end;
     end;
     inc(i);
   end;
 end;

end;

begin

 for var v in GetPrimes(12) do
   write(v, ' ');
 readln;

end.</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>

EchoLisp

The built-in prime-factors function performs the task. <lang lisp>(prime-factors 1024)

  → (2 2 2 2 2 2 2 2 2 2)

(lib 'bigint)

2^59 - 1

(prime-factors (1- (expt 2 59)))

   → (179951 3203431780337)

(prime-factors 100000000000000000037)

   → (31 821 66590107 59004541)</lang>

Eiffel

Uses the feature prime from the Task Primality by Trial Devision in the contract to check if the Result contains only prime numbers. <lang Eiffel>class

   PRIME_DECOMPOSITION

feature

   factor (p: INTEGER): ARRAY [INTEGER]
           -- Prime decomposition of 'p'.
       require
           p_positive: p > 0
       local
           div, i, next, rest: INTEGER
       do
           create Result.make_empty
           if p = 1 then
               Result.force (1, 1)
           end
           div := 2
           next := 3
           rest := p
           from
               i := 1
           until
               rest = 1
           loop
               from
               until
                   rest \\ div /= 0
               loop
                   Result.force (div, i)
                   rest := (rest / div).floor
                   i := i + 1
               end
               div := next
               next := next + 2
           end
       ensure
           is_divisor: across Result as r all p \\ r.item = 0 end
           is_prime: across Result as r all prime (r.item) end
       end</lang>

The test was done in an application class. (Similar as in other Eiffel examples (ex. Selectionsort).)

factor(5000)

Output:
2x2x2x5x5x5x5

Ela

Translation of: F#

<lang ela>open integer //arbitrary sized integers

decompose_prime n = loop n 2I

 where 
   loop c p | c < (p * p) = [c]
            | c % p == 0I = p :: (loop (c / p) p)
            | else = loop c (p + 1I)

decompose_prime 600851475143I</lang>

Output:
[71,839,1471,6857]

Elixir

<lang elixir>defmodule Prime do

 def decomposition(n), do: decomposition(n, 2, [])
 
 defp decomposition(n, k, acc) when n < k*k, do: Enum.reverse(acc, [n])
 defp decomposition(n, k, acc) when rem(n, k) == 0, do: decomposition(div(n, k), k, [k | acc])
 defp decomposition(n, k, acc), do: decomposition(n, k+1, acc)

end

prime = Stream.iterate(2, &(&1+1)) |>

       Stream.filter(fn n-> length(Prime.decomposition(n)) == 1 end) |>
       Enum.take(17)

mersenne = Enum.map(prime, fn n -> {n, round(:math.pow(2,n)) - 1} end) Enum.each(mersenne, fn {n,m} ->

 :io.format "~3s :~20w = ~s~n", ["M#{n}", m, Prime.decomposition(m) |> Enum.join(" x ")]

end)</lang>

Output:
 M2 :                   3 = 3
 M3 :                   7 = 7
 M5 :                  31 = 31
 M7 :                 127 = 127
M11 :                2047 = 23 x 89
M13 :                8191 = 8191
M17 :              131071 = 131071
M19 :              524287 = 524287
M23 :             8388607 = 47 x 178481
M29 :           536870911 = 233 x 1103 x 2089
M31 :          2147483647 = 2147483647
M37 :        137438953471 = 223 x 616318177
M41 :       2199023255551 = 13367 x 164511353
M43 :       8796093022207 = 431 x 9719 x 2099863
M47 :     140737488355327 = 2351 x 4513 x 13264529
M53 :    9007199254740991 = 6361 x 69431 x 20394401
M59 :  576460752303423487 = 179951 x 3203431780337

Erlang

<lang erlang>% no stack consuming version

factors(N) ->

    factors(N,2,[]).

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

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

factors(N,K,Acc) ->

   factors(N,K+1,Acc).</lang>

ERRE

<lang ERRE> PROGRAM DECOMPOSE


! ! for rosettacode.org !

!VAR NUM,J

DIM PF[100]

PROCEDURE STORE_FACTOR

  PF[0]=PF[0]+1
  PF[PF[0]]=CA
  I=I/CA

END PROCEDURE

PROCEDURE DECOMP(I)

 PF[0]=0  CA=2 ! special case
 LOOP
    IF I=1 THEN EXIT PROCEDURE END IF
    EXIT IF INT(I/CA)*CA<>I
    STORE_FACTOR
 END LOOP
 FOR CA=3 TO INT(SQR(I)) STEP 2 DO
    LOOP
       IF I=1 THEN EXIT PROCEDURE END IF
       EXIT IF INT(I/CA)*CA<>I
       STORE_FACTOR
    END LOOP
 END FOR
 IF I>1 THEN CA=I STORE_FACTOR END IF

END PROCEDURE

BEGIN

! ----- function generate
! in ...  I     ... number
! out ... PF[]  ... factors
!         PF[0] ... # of factors
! mod ... CA    ... pr.fact. candidate
PRINT(CHR$(12);) !CLS
INPUT("Numero ",NUM)
DECOMP(NUM)
PRINT(NUM;"=";)
FOR J=1 TO PF[0] DO
   PRINT(PF[J];)
END FOR
PRINT

END PROGRAM </lang> This version is a translation from Commodore BASIC program.

Ezhil

<lang Ezhil>

    1. இந்த நிரல் தரப்பட்ட எண்ணின் பகாஎண் கூறுகளைக் கண்டறியும்

நிரல்பாகம் பகாஎண்ணா(எண்1)

 ## இந்த நிரல்பாகம் தரப்பட்ட எண் பகு எண்ணா அல்லது பகா எண்ணா என்று கண்டறிந்து சொல்லும்
 ## பகுஎண் என்றால் 0 திரும்பத் தரப்படும்
 ## பகாஎண் என்றால் 1 திரும்பத் தரப்படும்
 @(எண்1 < 0) ஆனால்
  ## எதிர்மறை எண்களை நேராக்குதல்
   எண்1 = எண்1 * (-1)
 முடி
 @(எண்1 < 2) ஆனால்
  ## பூஜ்ஜியம், ஒன்று ஆகியவை பகா எண்கள் அல்ல
   பின்கொடு 0
 முடி
 @(எண்1 == 2) ஆனால்
   ## இரண்டு என்ற எண் ஒரு பகா எண்
   பின்கொடு 1
 முடி
 மீதம் = எண்1%2
 @(மீதம் == 0) ஆனால்
   ## இரட்டைப்படை எண், ஆகவே, இது பகா எண் அல்ல
   பின்கொடு 0
 முடி
   எண்1வர்க்கமூலம் = எண்1^0.5
   @(எண்2 = 3, எண்2 <= எண்1வர்க்கமூலம், எண்2 = எண்2 + 2) ஆக
     மீதம்1 = எண்1%எண்2
     @(மீதம்1 == 0) ஆனால்
       ## ஏதேனும் ஓர் எண்ணால் முழுமையாக வகுபட்டுவிட்டது, ஆகவே அது பகா எண் அல்ல
       பின்கொடு 0
     முடி
   முடி
   பின்கொடு 1

முடி

நிரல்பாகம் பகுத்தெடு(எண்1)

 ## இந்த எண் தரப்பட்ட எண்ணின் பகா எண் கூறுகளைக் கண்டறிந்து பட்டியல் இடும்
 கூறுகள் = பட்டியல்()
 @(எண்1 < 0) ஆனால்
   ## எதிர்மறை எண்களை நேராக்குதல்
   எண்1 = எண்1 * (-1)
 முடி
 @(எண்1 <= 1) ஆனால்
   ## ஒன்று அல்லது அதற்குக் குறைவான எண்களுக்குப் பகா எண் விகிதம் கண்டறியமுடியாது
   பின்கொடு கூறுகள்
 முடி
 
 @(பகாஎண்ணா(எண்1) == 1) ஆனால்
   ## தரப்பட்ட எண்ணே பகா எண்ணாக அமைந்துவிட்டால், அதற்கு அதுவே பகாஎண் கூறு ஆகும்
   பின்இணை(கூறுகள், எண்1)
   பின்கொடு கூறுகள்
 முடி
 தாற்காலிகஎண் = எண்1
 எண்2 = 2
 @(எண்2 <= தாற்காலிகஎண்) வரை
   விடை1 = பகாஎண்ணா(எண்2)
   மீண்டும்தொடங்கு = 0
   @(விடை1 == 1) ஆனால்
     விடை2 = தாற்காலிகஎண்%எண்2
     
     @(விடை2 == 0) ஆனால்
       ## பகா எண்ணால் முழுமையாக வகுபட்டுள்ளது, அதனைப் பட்டியலில் இணைக்கிறோம்
       பின்இணை(கூறுகள், எண்2)
       தாற்காலிகஎண் = தாற்காலிகஎண்/எண்2
       ## மீண்டும் இரண்டில் தொடங்கி இதே கணக்கிடுதலைத் தொடரவேண்டும்
       எண்2 = 2
       மீண்டும்தொடங்கு = 1
     முடி
     
   முடி
   @(மீண்டும்தொடங்கு == 0) ஆனால்
     ## அடுத்த எண்ணைத் தேர்ந்தெடுத்துக் கணக்கிடுதலைத் தொடரவேண்டும்
     எண்2 = எண்2 + 1
   முடி
 முடி
 பின்கொடு கூறுகள்

முடி

அ = int(உள்ளீடு("உங்களுக்குப் பிடித்த ஓர் எண்ணைத் தாருங்கள்: "))

பகாஎண்கூறுகள் = பட்டியல்()

பகாஎண்கூறுகள் = பகுத்தெடு(அ)

பதிப்பி "நீங்கள் தந்த எண்ணின் பகா எண் கூறுகள் இவை: ", பகாஎண்கூறுகள் </lang>

F#

<lang Fsharp>let decompose_prime n =

 let rec loop c p =
   if c < (p * p) then [c]
   elif c % p = 0I then p :: (loop (c/p) p)
   else loop c (p + 1I)

 loop n 2I
 

printfn "%A" (decompose_prime 600851475143I)</lang>

Output:
[71; 839; 1471; 6857]

Factor

factors from the math.primes.factors vocabulary converts a number into a sequence of its prime divisors; the rest of the code prints this sequence. <lang factor>USING: io kernel math math.parser math.primes.factors sequences ;

27720 factors [ number>string ] map " " join print ;</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>

FreeBASIC

<lang freebasic>' FB 1.05.0 Win64

Function isPrime(n As Integer) As Boolean

 If n Mod 2 = 0 Then Return n = 2
 If n Mod 3 = 0 Then Return n = 3
 Dim d As Integer = 5
 While d * d <= n
   If n Mod d = 0 Then Return False
   d += 2
   If n Mod d = 0 Then Return False
   d += 4
 Wend
 Return True

End Function

Sub getPrimeFactors(factors() As UInteger, n As UInteger)

 If n < 2 Then Return
 If isPrime(n) Then
   Redim factors(0 To 0)
   factors(0) = n
   Return
 End If
 Dim factor As UInteger = 2
 Do
   If n Mod factor = 0 Then
     Redim Preserve factors(0 To UBound(factors) + 1)
     factors(UBound(factors)) = factor
     n \= factor     
     If n = 1 Then Return
     If isPrime(n) Then factor = n
   Else
     factor += 1  
   End If    
 Loop

End Sub

Dim factors() As UInteger Dim primes(1 To 17) As UInteger = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59} Dim n As UInteger For i As UInteger = 1 To 17

 Erase factors
 n = 1 Shl primes(i) - 1
 getPrimeFactors factors(), n
 Print "2^";Str(primes(i)); Tab(5); " - 1 = "; Str(n); Tab(30);" => ";
 For j As UInteger = LBound(factors) To UBound(factors)
    Print factors(j);
    If j < UBound(factors) Then Print " x ";
 Next j
 Print

Next i Print Print "Press any key to quit" Sleep</lang>

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

Frink

Frink has a built-in factoring function which uses wheel factoring, trial division, Pollard p-1 factoring, and Pollard rho factoring. It also recognizes some special forms (e.g. Mersenne numbers) and handles them efficiently. <lang frink>println[factor[2^508-1]]</lang>

Output:

(total process time including JVM startup = 1.515 s)

[[3, 1], [5, 1], [509, 1], [18797, 1], [26417, 1], [72118729, 1], [140385293, 1], [2792688414613, 1], [8988357880501, 1], [90133566917913517709497, 1], [56713727820156410577229101238628035243, 1], [170141183460469231731687303715884105727, 1]]

Note that this means 31 * 51 * ...

GAP

Built-in function : <lang gap>FactorsInt(2^67-1);

  1. [ 193707721, 761838257287 ]</lang>

Or using the FactInt package : <lang gap>FactInt(2^67-1);

  1. [ [ 193707721, 761838257287 ], [ ] ]</lang>

Go

<lang go>package main

import (

   "fmt"
   "math/big"

)

var (

   ZERO = big.NewInt(0)
   ONE  = big.NewInt(1)

)

func Primes(n *big.Int) []*big.Int {

   res := []*big.Int{}
   mod, div := new(big.Int), new(big.Int)
   for i := big.NewInt(2); i.Cmp(n) != 1; {
       div.DivMod(n, i, mod)
       for mod.Cmp(ZERO) == 0 {
           res = append(res, new(big.Int).Set(i))
           n.Set(div)
           div.DivMod(n, i, mod)
       }
       i.Add(i, ONE)
   }
   return res

}

func main() {

   vals := []int64{
       1 << 31,
       1234567,
       333333,
       987653,
       2 * 3 * 5 * 7 * 11 * 13 * 17,
   }
   for _, v := range vals {
       fmt.Println(v, "->", Primes(big.NewInt(v)))
   }

}</lang>

Output:
2147483648 -> [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2]
1234567 -> [127 9721]
333333 -> [3 3 7 11 13 37]
987653 -> [29 34057]
510510 -> [2 3 5 7 11 13 17]

Groovy

This solution uses the fact that a given factor must be prime if no smaller factor divides it evenly, so it does not require an "isPrime-like function", assumed or otherwise. <lang groovy>def factorize = { long target ->

   if (target == 1) return [1L]

   if (target < 4) return [1L, target]

   def targetSqrt = Math.sqrt(target)
   def lowfactors = (2L..targetSqrt).findAll { (target % it) == 0 }
   if (lowfactors == []) return [1L, target]
   def nhalf = lowfactors.size() - ((lowfactors[-1]**2 == target) ? 1 : 0)

   [1] + lowfactors + (0..<nhalf).collect { target.intdiv(lowfactors[it]) }.reverse() + [target]

}

def decomposePrimes = { target ->

   def factors = factorize(target) - [1]
   def primeFactors = []
   factors.eachWithIndex { f, i ->
       if (i==0 || factors[0..<i].every {f % it != 0}) {
           primeFactors << f
           def pfPower = f*f
           while (target % pfPower == 0) {
               primeFactors << f
               pfPower *= f
           }
       }
   }
   primeFactors

}</lang>

Test #1:

<lang groovy>((1..30) + [97*4, 1000, 1024, 333333]).each { println ([number:it, primes:decomposePrimes(it)]) }</lang>

Output #1:
[number:1, primes:[]]
[number:2, primes:[2]]
[number:3, primes:[3]]
[number:4, primes:[2, 2]]
[number:5, primes:[5]]
[number:6, primes:[2, 3]]
[number:7, primes:[7]]
[number:8, primes:[2, 2, 2]]
[number:9, primes:[3, 3]]
[number:10, primes:[2, 5]]
[number:11, primes:[11]]
[number:12, primes:[2, 2, 3]]
[number:13, primes:[13]]
[number:14, primes:[2, 7]]
[number:15, primes:[3, 5]]
[number:16, primes:[2, 2, 2, 2]]
[number:17, primes:[17]]
[number:18, primes:[2, 3, 3]]
[number:19, primes:[19]]
[number:20, primes:[2, 2, 5]]
[number:21, primes:[3, 7]]
[number:22, primes:[2, 11]]
[number:23, primes:[23]]
[number:24, primes:[2, 2, 2, 3]]
[number:25, primes:[5, 5]]
[number:26, primes:[2, 13]]
[number:27, primes:[3, 3, 3]]
[number:28, primes:[2, 2, 7]]
[number:29, primes:[29]]
[number:30, primes:[2, 3, 5]]
[number:388, primes:[2, 2, 97]]
[number:1000, primes:[2, 2, 2, 5, 5, 5]]
[number:1024, primes:[2, 2, 2, 2, 2, 2, 2, 2, 2, 2]]
[number:333333, primes:[3, 3, 7, 11, 13, 37]]
Test #2:

<lang groovy>def isPrime = {factorize(it).size() == 2} (1..60).step(2).findAll(isPrime).each { println ([number:"2**${it}-1", value:2**it-1, primes:decomposePrimes(2**it-1)]) }</lang>

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

Perhaps a more sophisticated algorithm is in order. It took well over 1 hour to calculate the last three decompositions using this solution.

Haskell

The task description hints at using the isPrime function from the trial division task:

<lang haskell>factorize n = [ d | p <- [2..n], isPrime p, d <- divs n p ]

          -- [2..n] >>= (\p-> [p|isPrime p]) >>= divs n
   where
   divs n p | rem n p == 0 = p : divs (quot n p) p 
            | otherwise    = []</lang>

but it is not very efficient, to put it mildly. Inlining and fusing gets us the progressively more optimized <lang haskell>import Data.Maybe (listToMaybe) import Data.List (unfoldr)

factorize :: Integer -> [Integer] factorize n

 = unfoldr (\n     -> listToMaybe [(x, div n x)      | x <- [2..n], mod n x==0]) n
 = unfoldr (\(d,n) -> listToMaybe [(x, (x, div n x)) | x <- [d..n], mod n x==0]) (2,n)
 = unfoldr (\(d,n) -> listToMaybe [(x, (x, div n x)) | x <- 
                   takeWhile ((<=n).(^2)) [d..] ++ [n|n>1], mod n x==0]) (2,n)
 = unfoldr (\(ds,n) -> listToMaybe [(x, (dropWhile (< x) ds, div n x)) | n>1, x <-
                   takeWhile ((<=n).(^2)) ds ++ [n|n>1], mod n x==0]) (primesList,n)</lang>

The library function listToMaybe gets at most one element from its list argument. The last variant can be written as the optimal

<lang haskell>factorize n = divs n primesList

    where
    divs n ds@(d:t) | d*d > n    = [n | n > 1]
                    | r == 0     =  d : divs q ds
                    | otherwise  =      divs n t
           where  (q,r) = quotRem n d</lang>

See Sieve of Eratosthenes or Primality by trial division for a source of primes to use with this function. Actually as some other entries notice, with any ascending order list containing all primes (e.g. 2:[3,5..]) used in place of primesList, the factors found by this function are guaranteed to be prime, so no separate testing for primality is strictly needed; however using just primes is more efficient, if we already have them.

Output:
λ> mapM_ (print . factorize) $ take 11 [123123451..]
[11,41,273001]
[2,2,17,53,127,269]
[3,229,277,647]
[2,61561727]
[5,7,13,270601]
[2,2,2,2,2,2,2,2,3,3,3,47,379]
[37,109,30529]
[2,19,97,33403]
[3,3167,12959]
[2,2,5,6156173]
[123123461]

Icon and Unicon

<lang Icon>procedure main() factors := primedecomp(2^43-1) # a big int end

procedure primedecomp(n) #: return a list of factors local F,o,x F := []

every writes(o,n|(x := genfactors(n))) do {

  \o := "*"
  /o := "="
  put(F,x)   # build a list of factors to satisfy the task
  }

write() return F end

link factors</lang>

Uses genfactors and prime from factors

Sample Output showing factors of a large integer:

8796093022207=431*9719*2099863

J

<lang j>q:</lang>

Example use:

<lang j> q: 3684 2 2 3 307</lang>

and, more elaborately:

<lang j> _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. <lang java>private static final BigInteger two = BigInteger.valueOf(2);

public List<BigInteger> primeDecomp(BigInteger a) {

   // impossible for values lower than 2
   if (a.compareTo(two) < 0) {
       return null; 
   }
   //quickly handle even values
   List<BigInteger> result = new ArrayList<BigInteger>();
   while (a.and(BigInteger.ONE).equals(BigInteger.ZERO)) {
       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 (b.isProbablePrime(10)) {
               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>

Another alternate version designed to make fewer modular calculations: <lang java> private static final BigInteger TWO = BigInteger.valueOf(2); private static final BigInteger THREE = BigInteger.valueOf(3); private static final BigInteger FIVE = BigInteger.valueOf(5);

public static ArrayList<BigInteger> primeDecomp(BigInteger n){

   if(n.compareTo(TWO) < 0) return null;
   ArrayList<BigInteger> factors = new ArrayList<BigInteger>();
   
   // handle even values
   while(n.and(BigInteger.ONE).equals(BigInteger.ZERO)){
       n = n.shiftRight(1);
       factors.add(TWO);
   }
   
   // handle values divisible by three
   while(n.mod(THREE).equals(BigInteger.ZERO)){
       factors.add(THREE);
       n = n.divide(THREE);
   }
   
   // handle values divisible by five
   while(n.mod(FIVE).equals(BigInteger.ZERO)){
       factors.add(FIVE);
       n = n.divide(FIVE);
   }
   
   // much like how we can skip multiples of two, we can also skip
   // multiples of three and multiples of five. This increment array
   // helps us to accomplish that
   int[] pattern = {4,2,4,2,4,6,2,6};
   int pattern_index = 0;
   BigInteger current_test = BigInteger.valueOf(7);
   while(!n.equals(BigInteger.ONE)){
       while(n.mod(current_test).equals(BigInteger.ZERO)){
           factors.add(current_test);
           n = n.divide(current_test);
       }
       current_test = current_test.add(BigInteger.valueOf(pattern[pattern_index]));
       pattern_index = (pattern_index + 1) & 7;
   }
   
   return factors;

} </lang>

Translation of: C#

Simple but very inefficient method, because it will test divisibility of all numbers from 2 to max prime factor. When decomposing a large prime number this will take O(n) trial divisions instead of more common O(log n). <lang java>public static List<BigInteger> primeFactorBig(BigInteger a){

   List<BigInteger> ans = new LinkedList<BigInteger>();
   for(BigInteger divisor = BigInteger.valueOf(2);
       a.compareTo(ONE) > 0; divisor = divisor.add(ONE))
       while(a.mod(divisor).equals(ZERO)){
            ans.add(divisor);
            a = a.divide(divisor);
       }
   return ans;

}</lang>

JavaScript

This code uses the BigInteger Library jsbn and jsbn2 <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].equals(BigInteger.ZERO)) {
           if (prod) 
               output.value += "*"; 
           else 
               prod = true; 
           output.value += "2";
           n = qr[0];
       }
       else 
           break; 
   }
   while (!n.equals(BigInteger.ONE)) {
       var qr = n.divideAndRemainder(divisor);
       if (qr[1].equals(BigInteger.ZERO)) {
           if (prod) 
               output.value += "*"; 
           else 
               prod = true; 
           output.value += divisor;
           n = qr[0];
       }
       else 
           divisor = divisor.add(TWO); 
   }

}</lang>

Without any library. <lang javascript>function run_factorize(n) {

   if (n <= 3)
       return [n];
   var ans = [];
   var done = false;
   while (!done) {
       if (n % 2 === 0) {
           ans.push(2);
           n /= 2;
           continue;
       }
       if (n % 3 === 0) {
           ans.push(3);
           n /= 3;
           continue;
       }
       if (n === 1)
           return ans;
       var sr = Math.sqrt(n);
       done = true;
       // try to divide the checked number by all numbers till its square root.
       for (var i = 6; i <= (sr + 6); i += 6) {
           if (n % (i - 1) === 0) { // is n divisible by i-1?
               ans.push((i - 1));
               n /= (i - 1);
               done = false;
               break;
           }
           if (n % (i + 1) === 0) { // is n divisible by i+1?
               ans.push((i + 1));
               n /= (i + 1);
               done = false;
               break;
           }
       }
   }
   ans.push(n);
   return ans;

}</lang>

TDD using Jasmine

PrimeFactors.js <lang javascript>function factors(n) {

 if (!n || n < 2)
   return [];
 var f = [];
 for (var i = 2; i <= n; i++){
   while (n % i === 0){
     f.push(i);
     n /= i;
   }
 }
 return f;

}; </lang>

SpecPrimeFactors.js (with tag for Chutzpah) <lang javascript>/// <reference path="PrimeFactors.js" />

describe("Prime Factors", function() {

 it("Given nothing, empty is returned", function() {
   expect(factors()).toEqual([]);
 });
 it("Given 1, empty is returned", function() {
   expect(factors(1)).toEqual([]);
 });
 it("Given 2, 2 is returned", function() {
   expect(factors(2)).toEqual([2]);
 });
 it("Given 3, 3 is returned", function() {
   expect(factors(3)).toEqual([3]);
 });
 it("Given 4, 2 and 2 is returned", function() {
   expect(factors(4)).toEqual([2, 2]);
 });
 it("Given 5, 5 is returned", function() {
   expect(factors(5)).toEqual([5]);
 });
 it("Given 6, 2 and 3 is returned", function() {
   expect(factors(6)).toEqual([2, 3]);
 });
 it("Given 7, 7 is returned", function() {
   expect(factors(7)).toEqual([7]);
 });
 it("Given 8; 2, 2, and 2 is returned", function() {
   expect(factors(8)).toEqual([2, 2, 2]);
 });
 it("Given a large number, many primes factors are returned", function() {
   expect(factors(2*2*2*3*3*7*11*17))
     .toEqual([2, 2, 2, 3, 3, 7, 11, 17]);
 });
 it("Given a large prime number, that number is returned", function() {
   expect(factors(997)).toEqual([997]);
 });

}); </lang>

jq

Works with: jq version 1.5

Works with gojq, the Go implementation of jq

`factors` as defined below emits a stream of all the prime factors of the input integer. The implementation is compact, fast and space-efficient: no space is required to store the primes or factors already computed, there is no reliance on an "is_prime" function, and square roots are only computed if needed.

The economy comes about through the use of the builtin filter recurse/1, and the use of the state vector: [p, n, valid, sqrt], where p is the candidate factor, n is the number still to be factored, valid is a flag, and sqrt is either null or the square root of n.

gojq supports unlimited-precision integer arithmetic, but the C implementation of jq currently uses IEEE 754 64-bit numbers, so using the latter, the following program will only be reliable for integers up to and including 9,007,199,254,740,992 (2^53). However, "factors" could be easily modified to work with a "BigInt" library for jq, such as BigInt.jq. <lang jq>def factors:

 . as $in 
 | [2, $in, false]
 | recurse(
     . as [$p, $q, $valid, $s]
     | if $q == 1        then empty
       elif $q % $p == 0 then [$p, $q/$p, true]
       elif $p == 2      then [3, $q, false, $s]
       else ($s // ($q | sqrt)) as $s
       | if $p + 2 <= $s then [$p + 2, $q, false, $s]
         else [$q, 1, true]
         end
       end )
  | if .[2] then .[0] else empty end ;</lang> 

Examples: <lang jq> 24 | factors

  1. => 2 2 2 3

[9007199254740992 | factors] | length

  1. => 53
  1. 2**29-1 is 536870911

[ 536870911 | factors ]

  1. => [233,1103,2089]</lang>

Julia

using package Primes.jl: <lang julia> julia> Pkg.add("Primes") julia> factor(8796093022207) [9719=>1,431=>1,2099863=>1] </lang> (The factor function returns a dictionary whose keys are the factors and whose values are the multiplicity of each factor.)

Kotlin

<lang scala>// version 1.0.6

import java.math.BigInteger

val bigTwo = BigInteger.valueOf(2L) val bigThree = BigInteger.valueOf(3L)

fun getPrimeFactors(n: BigInteger): MutableList<BigInteger> {

   val factors = mutableListOf<BigInteger>()
   if (n < bigTwo) return factors
   if (n.isProbablePrime(20)) {
       factors.add(n)
       return factors
   }
   var factor = bigTwo
   var nn = n
   while (true) {
       if (nn % factor == BigInteger.ZERO) {
           factors.add(factor)
           nn /= factor
           if (nn == BigInteger.ONE) return factors
           if (nn.isProbablePrime(20)) factor = nn
       }
       else if (factor >= bigThree) factor += bigTwo 
       else factor = bigThree
   }

}

fun main(args: Array<String>) {

   val primes = intArrayOf(2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97)
   for (prime in primes) {
       val bigPow2 = bigTwo.pow(prime) - BigInteger.ONE
       println("2^${"%2d".format(prime)} - 1 = ${bigPow2.toString().padEnd(30)} => ${getPrimeFactors(bigPow2)}")
   }

}</lang>

Output:
2^ 2 - 1 = 3                              => [3]
2^ 3 - 1 = 7                              => [7]
2^ 5 - 1 = 31                             => [31]
2^ 7 - 1 = 127                            => [127]
2^11 - 1 = 2047                           => [23, 89]
2^13 - 1 = 8191                           => [8191]
2^17 - 1 = 131071                         => [131071]
2^19 - 1 = 524287                         => [524287]
2^23 - 1 = 8388607                        => [47, 178481]
2^29 - 1 = 536870911                      => [233, 1103, 2089]
2^31 - 1 = 2147483647                     => [2147483647]
2^37 - 1 = 137438953471                   => [223, 616318177]
2^41 - 1 = 2199023255551                  => [13367, 164511353]
2^43 - 1 = 8796093022207                  => [431, 9719, 2099863]
2^47 - 1 = 140737488355327                => [2351, 4513, 13264529]
2^53 - 1 = 9007199254740991               => [6361, 69431, 20394401]
2^59 - 1 = 576460752303423487             => [179951, 3203431780337]
2^61 - 1 = 2305843009213693951            => [2305843009213693951]
2^67 - 1 = 147573952589676412927          => [193707721, 761838257287]
2^71 - 1 = 2361183241434822606847         => [228479, 48544121, 212885833]
2^73 - 1 = 9444732965739290427391         => [439, 2298041, 9361973132609]
2^79 - 1 = 604462909807314587353087       => [2687, 202029703, 1113491139767]
2^83 - 1 = 9671406556917033397649407      => [167, 57912614113275649087721]
2^89 - 1 = 618970019642690137449562111    => [618970019642690137449562111]
2^97 - 1 = 158456325028528675187087900671 => [11447, 13842607235828485645766393]

Lambdatalk

<lang scheme> {def prime_fact.smallest

 {def prime_fact.smallest.r
 {lambda {:q :r :i}
  {if {and {> :r 0} {< :i :q}}
   then {prime_fact.smallest.r :q {% :q {+ :i 1}} {+ :i 1}}
   else :i}}}
 {lambda {:q} {prime_fact.smallest.r :q {% :q 2} 2}}}

{def prime_fact

{def prime_fact.r
 {lambda {:q :d}
  {if {> :q 1}
   then {let { {:q :q} {:d :d} 
               {:i {prime_fact.smallest :q}}}
             {prime_fact.r {floor {/ :q :i}} {#.push! :d :i}} }
   else {if {= {#.length :d} 1} then {b :d} else :d}}}}
{lambda {:n} :n:{prime_fact.r :n {#.new}}}}

{prime_fact {* 2 3 3 3 31 47 173}} -> 13611294:[2,3,3,3,31,47,173]

{map prime_fact {serie 2 101}} -> 2:[2] 3:[3] 4:[2,2] 5:[5] 6:[2,3] 7:[7] 8:[2,2,2] 9:[3,3] 10:[2,5] 11:[11] 12:[2,2,3] 13:[13] 14:[2,7] 15:[3,5] 16:[2,2,2,2] 17:[17] 18:[2,3,3] 19:[19] 20:[2,2,5] 21:[3,7] 22:[2,11] 23:[23] 24:[2,2,2,3] 25:[5,5] 26:[2,13] 27:[3,3,3] 28:[2,2,7] 29:[29] 30:[2,3,5] 31:[31] 32:[2,2,2,2,2] 33:[3,11] 34:[2,17] 35:[5,7] 36:[2,2,3,3] 37:[37] 38:[2,19] 39:[3,13] 40:[2,2,2,5] 41:[41] 42:[2,3,7] 43:[43] 44:[2,2,11] 45:[3,3,5] 46:[2,23] 47:[47] 48:[2,2,2,2,3] 49:[7,7] 50:[2,5,5] 51:[3,17] 52:[2,2,13] 53:[53] 54:[2,3,3,3] 55:[5,11] 56:[2,2,2,7] 57:[3,19] 58:[2,29] 59:[59] 60:[2,2,3,5] 61:[61] 62:[2,31] 63:[3,3,7] 64:[2,2,2,2,2,2] 65:[5,13] 66:[2,3,11] 67:[67] 68:[2,2,17] 69:[3,23] 70:[2,5,7] 71:[71] 72:[2,2,2,3,3] 73:[73] 74:[2,37] 75:[3,5,5] 76:[2,2,19] 77:[7,11] 78:[2,3,13] 79:[79] 80:[2,2,2,2,5] 81:[3,3,3,3] 82:[2,41] 83:[83] 84:[2,2,3,7] 85:[5,17] 86:[2,43] 87:[3,29] 88:[2,2,2,11] 89:[89] 90:[2,3,3,5] 91:[7,13] 92:[2,2,23] 93:[3,31] 94:[2,47] 95:[5,19] 96:[2,2,2,2,2,3] 97:[97] 98:[2,7,7] 99:[3,3,11] 100:[2,2,5,5] 101:[101] </lang>

LFE

<lang lisp> (defun factors (n)

 (factors n 2 '()))

(defun factors

 ((1 _ acc)
   acc)
 ((n k acc) (when (== 0 (rem n k)))
   (factors (div n k) k (cons k acc)))
 ((n k acc)
   (factors n (+ k 1) acc)))

</lang>

Lingo

<lang lingo>-- Returns list of prime factors for given number. -- To overcome the limits of integers (signed 32-bit in Lingo), -- the number can be specified as float (which works up to 2^53). -- For the same reason, values in returned list are floats, not integers. on getPrimeFactors (n)

 f = []
 f.sort()
 c = sqrt(n)
 i = 1.0
 repeat while TRUE
   i=i+1
   if i>c then exit repeat
   check = n/i
   if bitOr(check,0)=check then
     f.add(i)
     n = check
     c = sqrt(n)
     i = 1.0
   end if
 end repeat
 f.add(n)
 return f

end</lang> <lang lingo>put getPrimeFactors(12) -- [2.0000, 2.0000, 3.0000]

-- print floats without fractional digits the floatPrecision=0

put getPrimeFactors(12) -- [2, 2, 3]

put getPrimeFactors(1125899906842623.0) -- [3, 251, 601, 4051, 614141]</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>

Lua

The code of the used auxiliary function "IsPrime(n)" is located at Primality by trial division#Lua

<lang lua>function PrimeDecomposition( n )

   local f = {}
   
   if IsPrime( n ) then
       f[1] = n
       return f
   end
   local i = 2
   repeat
       while n % i == 0 do
           f[#f+1] = i
           n = n / i
       end
       
       repeat
           i = i + 1
       until IsPrime( i )       
   until n == 1
   
   return f

end</lang>

M2000 Interpreter

<lang M2000 Interpreter> Module Prime_decomposition {

     Inventory Known1=2@, 3@
     IsPrime=lambda  Known1 (x as decimal) -> {
                 =0=1
                 if exist(Known1, x) then =1=1 : exit
                 if x<=5 OR frac(x) then {if x == 2 OR x == 3 OR x == 5 then Append Known1, x  : =1=1
                 Break}
                 if frac(x/2) else exit
                 if frac(x/3) else exit
                 x1=sqrt(x):d = 5@
                 {if frac(x/d ) else exit
                       d += 2: if d>x1 then Append Known1, x : =1=1 : exit
                       if frac(x/d) else exit
                       d += 4: if d<= x1 else Append Known1, x :  =1=1: exit
                  loop}
           }
     decompose=lambda IsPrime (n as decimal) -> {
           Inventory queue Factors
           {
                k=2@
                While frac(n/k)=0 {
                n/=k
                     Append Factors, k
                }
                if n=1 then exit
                k++ 
                While frac(n/k)=0 {
                n/=k
                       Append Factors, k
                }   
                if n=1 then exit
                {
                k+=2
                while not isprime(k) {k+=2}
                While frac(n/k)=0 {
                n/=k
                       Append Factors, k
                }
                     if n=1 then exit
                     loop
                }             
           }
           =Factors
     }
     Data 10, 100, 12, 144, 496, 1212454
     while not empty {
       Print Decompose(Number)
     }

} Prime_decomposition </lang>

Maple

Maple has two commands for integer factorization: ifactor, which returns results in a form resembling textbook presentation and ifactors, which returns a list of two-element lists of prime factors and their multiplicities:

<lang Maple>> ifactor(1337);

                                  (7)  (191)

</lang> <lang Maple>> ifactors(1337);

                           [1, [[7, 1], [191, 1]]]

</lang>

Mathematica/Wolfram Language

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

Using the built-in function: <lang maxima>(%i1) display2d: false$ /* disable rendering exponents as superscripts */ (%i2) factor(2016); (%o2) 2^5*3^2*7 </lang> Using the underlying language: <lang maxima>prime_dec(n) := flatten(create_list(makelist(first(a), second(a)), a, ifactors(n)))$

/* or, slighlty more "functional" */ prime_dec(n) := flatten(map(lambda([a], apply(makelist, a)), ifactors(n)))$

prime_dec(2^4*3^5*5*7^2); /* [2, 2, 2, 2, 3, 3, 3, 3, 3, 5, 7, 7] */</lang>

MUMPS

<lang MUMPS>ERATO1(HI)

SET HI=HI\1
KILL ERATO1 ;Don't make it new - we want it to remain after the quit
NEW I,J,P
FOR I=2:1:(HI**.5)\1 DO
.FOR J=I*I:I:HI DO
..SET P(J)=1 ;$SELECT($DATA(P(J))#10:P(J)+1,1:1)
;WRITE !,"Prime numbers between 2 and ",HI,": "
FOR I=2:1:HI DO
.S:'$DATA(P(I)) ERATO1(I)=I ;WRITE $SELECT((I<3):"",1:", "),I
KILL I,J,P
QUIT

PRIMDECO(N)

;Returns its results in the string PRIMDECO
;Kill that before the first call to this recursive function
QUIT:N<=1
IF $D(PRIMDECO)=1 SET PRIMDECO="" D ERATO1(N)
SET N=N\1,I=0
FOR  SET I=$O(ERATO1(I)) Q:+I<1  Q:'(N#I)
IF I>1 SET PRIMDECO=$S($L(PRIMDECO)>0:PRIMDECO_"^",1:"")_I D PRIMDECO(N/I)
;that is, if I is a factor of N, add it to the string
QUIT</lang>
Usage:
USER>K ERATO1,PRIMDECO D PRIMDECO^ROSETTA(31415) W PRIMDECO
5^61^103
USER>K ERATO,PRIMDECO D PRIMDECO^ROSETTA(31318) W PRIMDECO
2^7^2237
USER>K ERATO,PRIMDECO D PRIMDECO^ROSETTA(34) W PRIMDECO
2^17
USER>K ERATO,PRIMDECO D PRIMDECO^ROSETTA(68) W PRIMDECO
2^2^17
USER>K ERATO,PRIMDECO D PRIMDECO^ROSETTA(7) W PRIMDECO
7
USER>K ERATO,PRIMDECO D PRIMDECO^ROSETTA(777) W PRIMDECO
3^7^37

Nim

Based on python floating point solution, but using integers rather than floats. <lang nim>import math, sequtils, strformat, strutils, times

proc getStep(n: int64): int64 {.inline.} =

 result = 1 + n shl 2 - n shr 1 shl 1

proc primeFac(n: int64): seq[int64] =

 var maxq = int64(sqrt(float(n)))
 var d = 1
 var q: int64 = 2 + (n and 1)   # Start with 2 or 3 according to oddity.
 while q <= maxq and n %% q != 0:
   q = getStep(d)
   inc d
 if q <= maxq:
   let q1 = primeFac(n /% q)
   let q2 = primeFac(q)
   result = concat(q2, q1, result)
 else:
   result.add(n)

iterator primes(limit: int): int =

 var isPrime = newSeq[bool](limit + 1)
 for n in 2..limit: isPrime[n] = true
 for n in 2..limit:
   if isPrime[n]:
     yield n
     for i in countup(n *% n, limit, n):
       isPrime[i] = false

when isMainModule:

 # Example: calculate factors of Mersenne numbers from M2 to M59.
 for m in primes(59):
   let p = 2i64^m - 1
   let s = &"2^{m}-1"
   stdout.write &"{s:<6} = {p} with factors: "
   let start = cpuTime()
   stdout.write primeFac(p).join(", ")
   echo &" => {(1000 * (cpuTime() - start)).toInt} ms"</lang>
Output:

Compiled with option -d:release

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

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

Oforth

Oforth handles aribitrary precision integers.

<lang Oforth>: factors(n) // ( aInteger -- aList ) | k p |

  ListBuffer new
  2 ->k
  n nsqrt ->p
  while( k p <= ) [
     n k /mod swap ifZero: [ 
        dup ->n nsqrt ->p 
        k over add continue
        ]
     drop k 1+ ->k
     ]
  n 1 > ifTrue: [ n over add ] 
  dup freeze ;</lang>
Output:
>2 128 pow 1 - dup println factors println
340282366920938463463374607431768211455
[3, 5, 17, 257, 641, 65537, 274177, 6700417, 67280421310721]
ok

PARI/GP

GP normally returns factored integers as a matrix with the first column representing the primes and the second their exponents. Thus factor(12)==[2,2;3,1] is true. But it's simple enough to convert this to a vector with repetition: <lang parigp>pd(n)={

 my(f=factor(n),v=f[,1]~);
 for(i=1,#v,
   while(f[i,2]--,
     v=concat(v,f[i,1])
   )
 );
 vecsort(v)

};</lang>

Pascal

<lang pascal>Program PrimeDecomposition(output);

type

 DynArray = array of integer;

procedure findFactors(n: Int64; var d: DynArray);

 var
   divisor, next, rest: Int64;
   i: integer;
begin
   i := 0;
   divisor := 2;
   next := 3;
   rest := n;
   while (rest <> 1) do
   begin
     while (rest mod divisor = 0) do
     begin
       setlength(d, i+1);
       d[i] := divisor;
       inc(i);
       rest := rest div divisor;
     end;
     divisor := next;
     next := next + 2;
   end;
 end;

var

 factors: DynArray;
 j: integer;

begin

 setlength(factors, 1);
 findFactors(1023*1024, factors);
 for j := low(factors) to high(factors) do
   writeln (factors[j]);

end.</lang>

Output:
% ./PrimeDecomposition
2
2
2
2
2
2
2
2
2
2
3
11
31

Optimization:

<lang pascal>Program PrimeDecomposition(output);

type

 DynArray = array of integer;

procedure findFactors(n: Int64; var d: DynArray);

 var
   divisor, next, rest: Int64;
   i: integer;
begin
   i := 0;
   divisor := 2;
   next := 3;
   rest := n;
   while (rest <> 1) do
   begin
     while (rest mod divisor = 0) do
     begin
       setlength(d, i+1);
       d[i] := divisor;
       inc(i);
       rest := rest div divisor;
     end;
     divisor := next;
     next := next + 2;  // try only odd numbers
     // cut condition: avoid many useless iterations
     if (rest < divisor * divisor) then
       begin
         setlength(d, i+1);
         d[i] := rest;
         rest := 1;
       end;
   end;
 end;

var

 factors: DynArray;
 j: integer;

begin

 setlength(factors, 1);
 findFactors(1023*1024, factors);
 for j := low(factors) to high(factors) do
   writeln (factors[j]);
 readln;

end.</lang>

Perl

These will work for large integers by adding the use bigint; clause.

Trivial trial division (very slow)

<lang perl>sub prime_factors {

   my ($n, $d, @out) = (shift, 1);
   while ($n > 1 && $d++) {
       $n /= $d, push @out, $d until $n % $d;
   }
   @out

}

print "@{[prime_factors(1001)]}\n";</lang>

Better trial division

This is much faster than the trivial version above. <lang perl>sub prime_factors {

 my($n, $p, @out) = (shift, 3);
 return if $n < 1;
 while (!($n&1)) { $n >>= 1; push @out, 2; }
 while ($n > 1 && $p*$p <= $n) {
   while ( ($n % $p) == 0) {
     $n /= $p;
     push @out, $p;
   }
   $p += 2;
 }
 push @out, $n if $n > 1;
 @out;

}</lang>

Modules

As usual, there are CPAN modules for this that will be much faster. These both take about 1 second to factor all Mersenne numbers from M_1 to M_150.

Library: ntheory

<lang perl>use ntheory qw/factor forprimes/; use bigint;

forprimes {

 my $p = 2 ** $_ - 1;
 print "2**$_-1: ", join(" ", factor($p)), "\n";

} 100, 150;</lang>

Output:
2^101-1: 7432339208719 341117531003194129
2^103-1: 2550183799 3976656429941438590393
2^107-1: 162259276829213363391578010288127
2^109-1: 745988807 870035986098720987332873
2^113-1: 3391 23279 65993 1868569 1066818132868207
2^127-1: 170141183460469231731687303715884105727
2^131-1: 263 10350794431055162386718619237468234569
2^137-1: 32032215596496435569 5439042183600204290159
2^139-1: 5625767248687 123876132205208335762278423601
2^149-1: 86656268566282183151 8235109336690846723986161

<lang perl>use Math::Pari qw/:int factorint isprime/;

  1. Convert Math::Pari's format into simple vector

sub factor {

 my ($pn,$pc) = @{Math::Pari::factorint(shift)};
 map { ($pn->[$_]) x $pc->[$_] } 0 .. $#$pn;

}

for (100 .. 150) {

 next unless isprime($_);
 my $p = 2 ** $_ - 1;
 print "2^$_-1: ", join(" ", factor($p)), "\n";

}</lang> With the same output.

Phix

For small numbers less than 253 on 32bit and 264 on 64bit just use prime_factors().

Library: Phix/mpfr
with javascript_semantics
requires("1.0.0")
include mpfr.e
atom t0 = time()
mpz z = mpz_init()
for i=1 to 25 do
    integer pi = get_prime(i)
    mpz_ui_pow_ui(z,2,pi)
    mpz_sub_ui(z,z,1)
    string zs = mpz_get_str(z),
           fs = mpz_factorstring(mpz_pollard_rho(z))
    if fs!=zs then zs &= " = "&fs end if
    printf(1,"2^%d-1 = %s\n",{pi,zs})
end for
string s = "600851475143"
for i=1 to 2 do
    mpz_set_str(z,s)
    printf(1,"%s = %s\n",{s,mpz_factorstring(mpz_pollard_rho(z))})
    s = "100000000000000000037"
end for
?elapsed(time()-t0)
Output:
2^2-1 = 3
2^3-1 = 7
2^5-1 = 31
2^7-1 = 127
2^11-1 = 2047 = 23*89
2^13-1 = 8191
2^17-1 = 131071
2^19-1 = 524287
2^23-1 = 8388607 = 47*178481
2^29-1 = 536870911 = 233*1103*2089
2^31-1 = 2147483647
2^37-1 = 137438953471 = 223*616318177
2^41-1 = 2199023255551 = 13367*164511353
2^43-1 = 8796093022207 = 431*9719*2099863
2^47-1 = 140737488355327 = 2351*4513*13264529
2^53-1 = 9007199254740991 = 6361*69431*20394401
2^59-1 = 576460752303423487 = 179951*3203431780337
2^61-1 = 2305843009213693951
2^67-1 = 147573952589676412927 = 193707721*761838257287
2^71-1 = 2361183241434822606847 = 228479*48544121*212885833
2^73-1 = 9444732965739290427391 = 439*2298041*9361973132609
2^79-1 = 604462909807314587353087 = 2687*202029703*1113491139767
2^83-1 = 9671406556917033397649407 = 167*57912614113275649087721
2^89-1 = 618970019642690137449562111
2^97-1 = 158456325028528675187087900671 = 11447*13842607235828485645766393
600851475143 = 71*839*1471*6857
100000000000000000037 = 31*821*59004541*66590107
"1.1s"

Picat

<lang Picat>go =>

 % Checking 2**prime-1
 foreach(P in primes(60))
    Factors = factors(2**P-1),
    println([n=2**P-1,factors=Factors])
 end,
 nl,
 % Testing a larger number
 println(factors(1361129467683753853853498429727072845823)),
 nl.

% % factors of N % factors(N) = Factors =>

  Factors = [],
  M = N,    
  while (M mod 2 == 0) 
     Factors := Factors ++ [2],  
     M := M div 2 
  end,
  T = 3,
  while (M > 1, T < 1+(sqrt(M)))
     if M mod T == 0 then
        [Divisors, NewM] = alldivisorsM(M, T),
        Factors := Factors ++ Divisors,
        M := NewM
     end,
     T := T + 2
  end,
  if M > 1 then Factors := Factors ++ [M] end.

alldivisorsM(N,Div) = [Divisors,M] =>

  M = N,
  Divisors = [],
  while (M mod Div == 0) 
     Divisors := Divisors ++ [Div],
     M := M div Div
  end.</lang>
Output:
[n = 3,factors = [3]]
[n = 7,factors = [7]]
[n = 31,factors = [31]]
[n = 127,factors = [127]]
[n = 2047,factors = [23,89]]
[n = 8191,factors = [8191]]
[n = 131071,factors = [131071]]
[n = 524287,factors = [524287]]
[n = 8388607,factors = [47,178481]]
[n = 536870911,factors = [233,1103,2089]]
[n = 2147483647,factors = [2147483647]]
[n = 137438953471,factors = [223,616318177]]
[n = 2199023255551,factors = [13367,164511353]]
[n = 8796093022207,factors = [431,9719,2099863]]
[n = 140737488355327,factors = [2351,4513,13264529]]
[n = 9007199254740991,factors = [6361,69431,20394401]]
[n = 576460752303423487,factors = [179951,3203431780337]]

[3,11,31,131,2731,8191,409891,7623851,145295143558111]

PicoLisp

The following solution generates a sequence of "trial divisors" (2 3 5 7 11 13 17 19 23 29 31 37 ..), as described by Donald E. Knuth, "The Art of Computer Programming", Vol.2, p.365. <lang PicoLisp>(de factor (N)

  (make
     (let (D 2  L (1 2 2 . (4 2 4 2 4 6 2 6 .))  M (sqrt N))
        (while (>= M D)
           (if (=0 (% N D))
              (setq M (sqrt (setq N (/ N (link D)))))
              (inc 'D (pop 'L)) ) )
        (link N) ) ) )

(factor 1361129467683753853853498429727072845823)</lang>

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

PL/I

<lang pli> test: procedure options (main, reorder);

  declare (n, i) fixed binary (31);
  get list (n);
  put edit ( n, '[' ) (x(1), a);

restart:

  if is_prime(n) then
     do;
        put edit (trim(n), ']' ) (x(1), a);
        stop;
     end;
  do i = n/2 to 2 by -1;
     if is_prime(i) then
        if (mod(n, i) = 0) then
           do;
              put edit ( trim(i) ) (x(1), a);
              n = n / i;
              go to restart;
           end;
  end;
  put edit ( ' ]' ) (a);

is_prime: procedure (n) options (reorder) returns (bit(1));

  declare n fixed binary (31);
  declare i fixed binary (31);
  if n < 2 then return ('0'b);
  if n = 2 then return ('1'b);
  if mod(n, 2) = 0 then return ('0'b);
  do i = 3 to sqrt(n) by 2;
     if mod(n, i) = 0 then return ('0'b);
  end;
  return ('1'b);

end is_prime;

end test; </lang>

Results from various runs:
        1234567 [ 9721 127 ]
          32768 [ 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ]
             99 [ 11 3 3 ]
        9876543 [ 14503 227 3 ]
            100 [ 5 5 2 2 ]
        9999999 [ 4649 239 3 3 ]
           5040 [ 7 5 3 3 2 2 2 2 ]

PowerShell

<lang PowerShell> function eratosthenes ($n) {

   if($n -gt 1){
       $prime = @(1..($n+1) | foreach{$true})
       $prime[1] = $false
       $m = [Math]::Floor([Math]::Sqrt($n))
       function multiple($i) {
           for($j = $i*$i; $j -le $n; $j += $i) {
               $prime[$j] = $false
           }
       }
       multiple 2
       for($i = 3; $i -le $m; $i += 2) {
           if($prime[$i]) {multiple $i}
       }
       1..$n | where{$prime[$_]}
   } else {
       Write-Error "$n is not greater than 1"
   }

} function prime-decomposition ($n) {

   $array = eratosthenes $n
   $prime = @()
   foreach($p in $array) {
       while($n%$p -eq 0) {
           $n /= $p
           $prime += @($p)
       }
   }
   $prime

} "$(prime-decomposition 12)" "$(prime-decomposition 100)" </lang> Output:

2 2 3
2 2 5 5

Prolog

<lang Prolog>prime_decomp(N, L) :-

   SN is sqrt(N),
   prime_decomp_1(N, SN, 2, [], L).


prime_decomp_1(1, _, _, L, L) :- !.

% Special case for 2, increment 1 prime_decomp_1(N, SN, D, L, LF) :-

   (   0 is N mod D ->
       Q is N / D,
       SQ is sqrt(Q),
       prime_decomp_1(Q, SQ, D, [D |L], LF)
   ;
       D1 is D+1,
       (   D1 > SN ->
           LF = [N |L]
       ;
           prime_decomp_2(N, SN, D1, L, LF)
       )
   ).

% General case, increment 2 prime_decomp_2(1, _, _, L, L) :- !.

prime_decomp_2(N, SN, D, L, LF) :-

   (   0 is N mod D ->
       Q is N / D,
       SQ is sqrt(Q),
       prime_decomp_2(Q, SQ, D, [D |L], LF);
       D1 is D+2,
       (   D1 > SN ->
           LF = [N |L]
       ;
           prime_decomp_2(N, SN, D1, L, LF)
       )
   ).</lang>
Output:

<lang Prolog> ?- time(prime_decomp(9007199254740991, L)). % 138,882 inferences, 0.344 CPU in 0.357 seconds (96% CPU, 404020 Lips) L = [20394401,69431,6361].

?- time(prime_decomp(576460752303423487, L)).

% 2,684,734 inferences, 0.672 CPU in 0.671 seconds (100% CPU, 3995883 Lips) L = [3203431780337,179951].

?- time(prime_decomp(1361129467683753853853498429727072845823, L)).

% 18,080,807 inferences, 7.953 CPU in 7.973 seconds (100% CPU, 2273422 Lips) L = [145295143558111,7623851,409891,8191,2731,131,31,11,3].</lang>

Simple version

Translation of: Erlang

Optimized to stop on square root, and count by +2 on odds, above 2.

<lang Prolog>factors( N, FS):-

   factors2( N, FS).

factors2( N, FS):-

   ( N < 2        -> FS = [] 
   ; 4 > N        -> FS = [N] 
   ; 0 is N rem 2 -> FS = [K|FS2], N2 is N div 2, factors2( N2, FS2)
   ;                 factors( N, 3, FS)
   ).

factors( N, K, FS):-

   ( N < 2        -> FS = [] 
   ; K*K > N      -> FS = [N] 
   ; 0 is N rem K -> FS = [K|FS2], N2 is N div K, factors( N2, K, FS2)
   ;                 K2 is K+2, factors( N, K2, FS)
   ).</lang>

Expression Tree version

Uses a 2*3*5*7 factor wheel, but the main feature is that it returns the decomposition as a fully simplified expression tree. <lang Prolog> wheel2357(L) :-

   W = [2,  4,  2,  4,  6,  2,  6,  4,
        2,  4,  6,  6,  2,  6,  4,  2,
        6,  4,  6,  8,  4,  2,  4,  2,
        4,  8,  6,  4,  6,  2,  4,  6,
        2,  6,  6,  4,  2,  4,  6,  2,
        6,  4,  2,  4,  2, 10,  2, 10 | W],
   L = [1, 2, 2, 4 | W].

factor(1, 1) :- !. factor(N, Fac) :-

   N > 1,
   wheel2357(W),
   factor(N, 2, W, 1, Fac0),
   reverse_factors(Fac0, Fac).

factor(N, F, _, Fac1, Fac2) :- F*F > N, !, add_factor(N, Fac1, Fac2). factor(N, F, W, Fac1, Fac) :-

   divmod(N, F, Q, 0), !,
   add_factor(F, Fac1, Fac2),
   factor(Q, F, W, Fac2, Fac).

factor(N, F1, [A|As], Fac1, Fac) :-

   F2 is F1 + A,
   factor(N, F2, As, Fac1, Fac).

add_factor(F, 1, F) :- !. add_factor(F, F, F**2) :- !. add_factor(F, F**Ex1, F**Ex2) :- succ(Ex1, Ex2), !.

add_factor(F, F*A, F**2*A) :- !. add_factor(F, F**Ex1*Rest, F**Ex2*Rest) :- succ(Ex1, Ex2), !. add_factor(F, Fac, F*Fac).

reverse_factors(A*B, C*A) :- reverse_factors(B, C), !. reverse_factors(A, A). </lang>

Output:
?- factor(277,X).
X = 277.

?- factor(1003,X).
X = 17*59.

?- factor(1024,X).
X = 2**10.

?- factor(768,X).
X = 2**8*3.

?- factor(1361129467683753853853498429727072845823,X).
X = 3*11*31*131*2731*8191*409891*7623851*145295143558111.

?- factor(360,X).
X = 2**3*3**2*5.

Pure

<lang pure>factor n = factor 2 n with

 factor k n = k : factor k (n div k) if n mod k == 0;
        = if n>1 then [n] else [] if k*k>n;
        = factor (k+1) n if k==2;
        = factor (k+2) n otherwise;

end;</lang>

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>

Output:
Factored 9007199254740991 in 0.27 seconds.
6361
69431
20394401

Python

Python: Using Croft Spiral sieve

Note: the program below is saved to file prime_decomposition.py and imported as a library here, here, here, here and here.

<lang python>from __future__ import print_function

import sys from itertools import islice, cycle, count

try:

   from itertools import compress

except ImportError:

   def compress(data, selectors):
       """compress('ABCDEF', [1,0,1,0,1,1]) --> A C E F"""
       return (d for d, s in zip(data, selectors) if s)


def is_prime(n):

   return list(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 croft():

   """Yield prime integers using the Croft Spiral sieve.
   This is a variant of wheel factorisation modulo 30.
   """
   # Copied from:
   #   https://code.google.com/p/pyprimes/source/browse/src/pyprimes.py
   # Implementation is based on erat3 from here:
   #   http://stackoverflow.com/q/2211990
   # and this website:
   #   http://www.primesdemystified.com/
   # Memory usage increases roughly linearly with the number of primes seen.
   # dict ``roots`` stores an entry x:p for every prime p.
   for p in (2, 3, 5):
       yield p
   roots = {9: 3, 25: 5}  # Map d**2 -> d.
   primeroots = frozenset((1, 7, 11, 13, 17, 19, 23, 29))
   selectors = (1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0)
   for q in compress(
           # Iterate over prime candidates 7, 9, 11, 13, ...
           islice(count(7), 0, None, 2),
           # Mask out those that can't possibly be prime.
           cycle(selectors)
           ):
       # Using dict membership testing instead of pop gives a
       # 5-10% speedup over the first three million primes.
       if q in roots:
           p = roots[q]
           del roots[q]
           x = q + 2*p
           while x in roots or (x % 30) not in primeroots:
               x += 2*p
           roots[x] = p
       else:
           roots[q*q] = q
           yield q

primes = croft

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


if __name__ == '__main__':

   # Example: calculate factors of Mersenne numbers to M59 #

   import time

   for m in primes():
       p = 2 ** m - 1
       print( "2**{0:d}-1 = {1:d}, with factors:".format(m, p) )
       start = time.time()
       for factor in decompose(p):
           print(factor, end=' ')
           sys.stdout.flush()

       print( "=> {0:.2f}s".format( 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.01s
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.00s
2**23-1 = 8388607, with factors:
47 178481 => 0.01s
2**29-1 = 536870911, with factors:
233 1103 2089 => 0.01s
2**31-1 = 2147483647, with factors:
2147483647 => 0.03s
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.01s
2**53-1 = 9007199254740991, with factors:
6361 69431 20394401 => 0.04s
2**59-1 = 576460752303423487, with factors:
179951 3203431780337 => 1.22s

Python: Using floating point

Here a shorter and marginally faster algorithm:

<lang python>from math import floor, sqrt try:

   long

except NameError:

   long = int

def fac(n):

   step = lambda x: 1 + (x<<2) - ((x>>1)<<1)
   maxq = long(floor(sqrt(n)))
   d = 1
   q = 2 if n % 2 == 0 else 3 
   while q <= maxq and n % q != 0:
       q = step(d)
       d += 1
   return [q] + fac(n // q) if q <= maxq else [n]

if __name__ == '__main__':

   import time
   start = time.time()
   tocalc =  2**59-1
   print("%s = %s" % (tocalc, fac(tocalc)))
   print("Needed %ss" % (time.time() - start))</lang>
Output:
576460752303423487 = [3203431780337, 179951]
Needed 0.9240529537200928s

Quackery

<lang Quackery> [ [] swap

   dup times
   [ [ dup i^ 2 + /mod
       0 = while
       nip dip 
         [ i^ 2 + join ]
       again ]
     drop
     dup 1 = if conclude ] 
   drop ]                  is primefactors ( n --> [ )
1047552 primefactors echo</lang>
Output:
[ 2 2 2 2 2 2 2 2 2 2 3 11 31 ]

R

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

 x <- NULL
 firstprime<- 2; secondprime <- 3; everyprime <- num
 while( everyprime != 1 ) {
   while( everyprime%%firstprime == 0 ) {
     x <- c(x, firstprime)
     everyprime <- floor(everyprime/ firstprime)
   }
   firstprime <- secondprime
   secondprime <- secondprime + 2
 }
 x

}

print(findfactors(1027*4))</lang>

Or a more explicit (but less efficient) recursive approach:

Recursive Approach (Less efficient for large numbers)

<lang R> primes <- as.integer(c())

max_prime_checker <- function(n){

 divisor <<- NULL
 primes <- primes[primes <= n]
 for(i in 1:length(primes)){
   if((n/primes[i]) %% 1 == 0){
     divisor[i]<<-1
   } else {
     divisor[i]<<-0
   } 
 }
 num_find <<- primes*as.integer(divisor)
 
 return(max(num_find))

}

  1. recursive prime finder

prime_factors <- function(n){

 factors <- NULL
 
 large <- max_prime_checker(n)
 n1 <- n/large 
 
 if(max_prime_checker(n1) == n1){
   factors <- c(large,n1)
   return(factors)
 } else {
   factors <- c(large, prime_factors(n1))
   return(factors)
 }

} </lang>

Alternate solution

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

 a <- NULL
 if (n > 1) {
   while (n %% 2 == 0) {
     a <- c(a, 2)
     n <- n %/% 2
   }
   k <- 3
   while (k * k <= n) {
     while (n %% k == 0) {
       a <- c(a, k)
       n <- n %/% k
     }
     k <- k + 2
   }
   if (n > 1) a <- c(a, n)
 }
 a

}</lang>

Racket

<lang Racket>

  1. lang racket

(require math) (define (factors n)

 (append-map (λ (x) (make-list (cadr x) (car x))) (factorize n)))

</lang>

Or, an explicit (and less efficient) computation: <lang Racket>

  1. lang racket

(define (factors number)

 (let loop ([n number] [i 2])
   (if (= n 1)
     '()
     (let-values ([(q r) (quotient/remainder n i)])
       (if (zero? r) (cons i (loop q i)) (loop n (add1 i)))))))

</lang>

Raku

(formerly Perl 6)

Pure Raku

This is a pure Raku version that uses no outside libraries. It uses a variant of Pollard's rho factoring algorithm and is fairly performent when factoring numbers < 2⁸⁰; typically taking well under a second on an i7. It starts to slow down with larger numbers, but really bogs down factoring numbers that have more than 1 factor larger than about 2⁴⁰.

<lang perl6>sub prime-factors ( Int $n where * > 0 ) {

   return $n if $n.is-prime;
   return () if $n == 1;
   my $factor = find-factor( $n );
   sort flat ( $factor, $n div $factor ).map: &prime-factors;

}

sub find-factor ( Int $n, $constant = 1 ) {

   return 2 unless $n +& 1;
   if (my $gcd = $n gcd 6541380665835015) > 1 { # magic number: [*] primes 3 .. 43
       return $gcd if $gcd != $n
   }
   my $x      = 2;
   my $rho    = 1;
   my $factor = 1;
   while $factor == 1 {
       $rho = $rho +< 1;
       my $fixed = $x;
       my int $i = 0;
       while $i < $rho {
           $x = ( $x * $x + $constant ) % $n;
           $factor = ( $x - $fixed ) gcd $n;
           last if 1 < $factor;
           $i = $i + 1;
       }
   }
   $factor = find-factor( $n, $constant + 1 ) if $n == $factor;
   $factor;

}

.put for (2²⁹-1, 2⁴¹-1, 2⁵⁹-1, 2⁷¹-1, 2⁷⁹-1, 2⁹⁷-1, 2¹¹⁷-1, 2²⁴¹-1, 5465610891074107968111136514192945634873647594456118359804135903459867604844945580205745718497)\ .hyper(:1batch).map: -> $n {

   my $start = now;
  "factors of $n: ",
   prime-factors($n).join(' × '), " \t in ", (now - $start).fmt("%0.3f"), " sec."

}</lang>

Output:
factors of 536870911:  233 × 1103 × 2089    in  0.004  sec.
factors of 2199023255551:  13367 × 164511353     in  0.011  sec.
factors of 576460752303423487:  179951 × 3203431780337       in  0.023  sec.
factors of 2361183241434822606847:  228479 × 48544121 × 212885833    in  0.190  sec.
factors of 604462909807314587353087:  2687 × 202029703 × 1113491139767       in  0.294  sec.
factors of 158456325028528675187087900671:  11447 × 13842607235828485645766393       in  0.005  sec.
factors of 166153499473114484112975882535043071:  7 × 73 × 79 × 937 × 6553 × 8191 × 86113 × 121369 × 7830118297      in  0.022  sec.
factors of 3533694129556768659166595001485837031654967793751237916243212402585239551:  22000409 × 160619474372352289412737508720216839225805656328990879953332340439     in  0.085  sec.
factors of 5465610891074107968111136514192945634873647594456118359804135903459867604844945580205745718497:  165901 × 10424087 × 18830281 × 53204737 × 56402249 × 59663291 × 91931221 × 95174413 × 305293727939 × 444161842339 × 790130065009     in  28.427  sec.

There is a Raku module available: Prime::Factor, that uses essentially this algorithm with some minor performance tweaks.

External library

If you really need a speed boost, load the highly optimized Perl 5 ntheory module. It needs a little extra plumbing to deal with the lack of built-in big integer support, but for large number factoring the interface overhead is worth it. <lang perl6>use Inline::Perl5; my $p5 = Inline::Perl5.new(); $p5.use( 'ntheory' );

sub prime-factors ($i) {

   my &primes = $p5.run('sub { map { ntheory::todigitstring $_ } sort {$a <=> $b} ntheory::factor $_[0] }');
   primes("$i");

}

for 2²⁹-1, 2⁴¹-1, 2⁵⁹-1, 2⁷¹-1, 2⁷⁹-1, 2⁹⁷-1, 2¹¹⁷-1, 5465610891074107968111136514192945634873647594456118359804135903459867604844945580205745718497

->  $n {
   my $start = now;
   say "factors of $n: ",
   prime-factors($n).join(' × '), " \t in ", (now - $start).fmt("%0.3f"), " sec."

}</lang>

Output:
factors of 536870911: 233 × 1103 × 2089     in 0.001 sec.
factors of 2199023255551: 13367 × 164511353      in 0.001 sec.
factors of 576460752303423487: 179951 × 3203431780337    in 0.001 sec.
factors of 2361183241434822606847: 228479 × 48544121 × 212885833     in 0.012 sec.
factors of 604462909807314587353087: 2687 × 202029703 × 1113491139767    in 0.003 sec.
factors of 158456325028528675187087900671: 11447 × 13842607235828485645766393    in 0.001 sec.
factors of 166153499473114484112975882535043071: 7 × 73 × 79 × 937 × 6553 × 8191 × 86113 × 121369 × 7830118297   in 0.001 sec.
factors of 5465610891074107968111136514192945634873647594456118359804135903459867604844945580205745718497: 165901 × 10424087 × 18830281 × 53204737 × 56402249 × 59663291 × 91931221 × 95174413 × 305293727939 × 444161842339 × 790130065009      in 0.064 sec.

REXX

optimized slightly

No (error) checking was done for the input arguments to test their validity.

The number of decimal digits is adjusted to match the size of the top-of-the-range (top).

Also, a count of primes found is shown.

If the   top   number is negative, only the number of primes up to   abs(top)   is shown.

A method exists in this REXX program to also test Mersenne-type numbers   (2n - 1).

Since the majority of computing time is spent looking for primes, that part of the program was
optimized somewhat (but could be extended if more optimization is wanted). <lang rexx>/*REXX pgm does prime decomposition of a range of positive integers (with a prime count)*/ numeric digits 1000 /*handle thousand digits for the powers*/ parse arg bot top step base add /*get optional arguments from the C.L. */ if bot== then do; bot=1; top=100; end /*no BOT given? Then use the default.*/ if top== then top=bot /* " TOP? " " " " " */ if step== then step= 1 /* " STEP? " " " " " */ if add == then add= -1 /* " ADD? " " " " " */ tell= top>0; top=abs(top) /*if TOP is negative, suppress displays*/ w=length(top) /*get maximum width for aligned display*/ if base\== then w=length(base**top) /*will be testing powers of two later? */ @.=left(, 7); @.0="{unity}"; @.1='[prime]' /*some literals: pad; prime (or not).*/ numeric digits max(9, w+1) /*maybe increase the digits precision. */

  1. =0 /*#: is the number of primes found. */
       do n=bot  to top  by step                /*process a single number  or  a range.*/
       ?=n;  if base\==  then ?=base**n + add /*should we perform a "Mercenne" test? */
       pf=factr(?);      f=words(pf)            /*get prime factors; number of factors.*/
       if f==1  then #=#+1                      /*Is N prime?  Then bump prime counter.*/
       if tell  then say right(?,w)   right('('f")",9)   'prime factors: '     @.f     pf
       end   /*n*/

say ps= 'primes'; if p==1 then ps= "prime" /*setup for proper English in sentence.*/ say right(#, w+9+1) ps 'found.' /*display the number of primes found. */ exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ factr: procedure; parse arg x 1 d,$ /*set X, D to argument 1; $ to null.*/ if x==1 then return /*handle the special case of X = 1. */

      do  while x//2==0;  $=$ 2;  x=x%2;  end   /*append all the  2  factors of new  X.*/
      do  while x//3==0;  $=$ 3;  x=x%3;  end   /*   "    "   "   3     "     "  "   " */
      do  while x//5==0;  $=$ 5;  x=x%5;  end   /*   "    "   "   5     "     "  "   " */
      do  while x//7==0;  $=$ 7;  x=x%7;  end   /*   "    "   "   7     "     "  "   " */
                                                /*                                  ___*/

q=1; do while q<=x; q=q*4; end /*these two lines compute integer √ X */ r=0; do while q>1; q=q%4; _=d-r-q; r=r%2; if _>=0 then do; d=_; r=r+q; end; end

      do j=11  by 6  to r                       /*insure that  J  isn't divisible by 3.*/
      parse var j    -1  _                    /*obtain the last decimal digit of  J. */
      if _\==5  then  do  while x//j==0;  $=$ j;  x=x%j;  end     /*maybe reduce by J. */
      if _ ==3  then iterate                    /*Is next  Y  is divisible by 5?  Skip.*/
      y=j+2;          do  while x//y==0;  $=$ y;  x=x%y;  end     /*maybe reduce by J. */
      end   /*j*/
                                                /* [↓]  The $ list has a leading blank.*/

if x==1 then return $ /*Is residual=unity? Then don't append.*/

             return $ x                         /*return   $   with appended residual. */</lang>

output   when using the default input of:   1   100

(Shown at three-quarter size.)

  1       (0) prime factors:  {unity}
  2       (1) prime factors:  [prime]  2
  3       (1) prime factors:  [prime]  3
  4       (2) prime factors:           2 2
  5       (1) prime factors:  [prime]  5
  6       (2) prime factors:           2 3
  7       (1) prime factors:  [prime]  7
  8       (3) prime factors:           2 2 2
  9       (2) prime factors:           3 3
 10       (2) prime factors:           2 5
 11       (1) prime factors:  [prime]  11
 12       (3) prime factors:           2 2 3
 13       (1) prime factors:  [prime]  13
 14       (2) prime factors:           2 7
 15       (2) prime factors:           3 5
 16       (4) prime factors:           2 2 2 2
 17       (1) prime factors:  [prime]  17
 18       (3) prime factors:           2 3 3
 19       (1) prime factors:  [prime]  19
 20       (3) prime factors:           2 2 5
 21       (2) prime factors:           3 7
 22       (2) prime factors:           2 11
 23       (1) prime factors:  [prime]  23
 24       (4) prime factors:           2 2 2 3
 25       (2) prime factors:           5 5
 26       (2) prime factors:           2 13
 27       (3) prime factors:           3 3 3
 28       (3) prime factors:           2 2 7
 29       (1) prime factors:  [prime]  29
 30       (3) prime factors:           2 3 5
 31       (1) prime factors:  [prime]  31
 32       (5) prime factors:           2 2 2 2 2
 33       (2) prime factors:           3 11
 34       (2) prime factors:           2 17
 35       (2) prime factors:           5 7
 36       (4) prime factors:           2 2 3 3
 37       (1) prime factors:  [prime]  37
 38       (2) prime factors:           2 19
 39       (2) prime factors:           3 13
 40       (4) prime factors:           2 2 2 5
 41       (1) prime factors:  [prime]  41
 42       (3) prime factors:           2 3 7
 43       (1) prime factors:  [prime]  43
 44       (3) prime factors:           2 2 11
 45       (3) prime factors:           3 3 5
 46       (2) prime factors:           2 23
 47       (1) prime factors:  [prime]  47
 48       (5) prime factors:           2 2 2 2 3
 49       (2) prime factors:           7 7
 50       (3) prime factors:           2 5 5
 51       (2) prime factors:           3 17
 52       (3) prime factors:           2 2 13
 53       (1) prime factors:  [prime]  53
 54       (4) prime factors:           2 3 3 3
 55       (2) prime factors:           5 11
 56       (4) prime factors:           2 2 2 7
 57       (2) prime factors:           3 19
 58       (2) prime factors:           2 29
 59       (1) prime factors:  [prime]  59
 60       (4) prime factors:           2 2 3 5
 61       (1) prime factors:  [prime]  61
 62       (2) prime factors:           2 31
 63       (3) prime factors:           3 3 7
 64       (6) prime factors:           2 2 2 2 2 2
 65       (2) prime factors:           5 13
 66       (3) prime factors:           2 3 11
 67       (1) prime factors:  [prime]  67
 68       (3) prime factors:           2 2 17
 69       (2) prime factors:           3 23
 70       (3) prime factors:           2 5 7
 71       (1) prime factors:  [prime]  71
 72       (5) prime factors:           2 2 2 3 3
 73       (1) prime factors:  [prime]  73
 74       (2) prime factors:           2 37
 75       (3) prime factors:           3 5 5
 76       (3) prime factors:           2 2 19
 77       (2) prime factors:           7 11
 78       (3) prime factors:           2 3 13
 79       (1) prime factors:  [prime]  79
 80       (5) prime factors:           2 2 2 2 5
 81       (4) prime factors:           3 3 3 3
 82       (2) prime factors:           2 41
 83       (1) prime factors:  [prime]  83
 84       (4) prime factors:           2 2 3 7
 85       (2) prime factors:           5 17
 86       (2) prime factors:           2 43
 87       (2) prime factors:           3 29
 88       (4) prime factors:           2 2 2 11
 89       (1) prime factors:  [prime]  89
 90       (4) prime factors:           2 3 3 5
 91       (2) prime factors:           7 13
 92       (3) prime factors:           2 2 23
 93       (2) prime factors:           3 31
 94       (2) prime factors:           2 47
 95       (2) prime factors:           5 19
 96       (6) prime factors:           2 2 2 2 2 3
 97       (1) prime factors:  [prime]  97
 98       (3) prime factors:           2 7 7
 99       (3) prime factors:           3 3 11
100       (4) prime factors:           2 2 5 5

           25 primes found.

output   when using the input of:   9007199254740991

9007199254740991       (3) prime factors:           6361 69431 20394401

              0 primes found.

output   when using the input of:   2543821448263974486045199

2543821448263974486045199       (6) prime factors:           701 1123 1123 2411 1092461 1092461

              0 primes found.

output   when using the input of:   1   -1000000

            78498 primes found.

output   when using the input of:   2   50   1   2   -1

(essentially testing for Mersenne primes:   2n -1)

(Shown at three-quarter size.)

               3       (1) prime factors:  [prime]  3
               7       (1) prime factors:  [prime]  7
              15       (2) prime factors:           3 5
              31       (1) prime factors:  [prime]  31
              63       (3) prime factors:           3 3 7
             127       (1) prime factors:  [prime]  127
             255       (3) prime factors:           3 5 17
             511       (2) prime factors:           7 73
            1023       (3) prime factors:           3 11 31
            2047       (2) prime factors:           23 89
            4095       (5) prime factors:           3 3 5 7 13
            8191       (1) prime factors:  [prime]  8191
           16383       (2) prime factors:           3 5461
           32767       (2) prime factors:           7 4681
           65535       (4) prime factors:           3 5 17 257
          131071       (1) prime factors:  [prime]  131071
          262143       (5) prime factors:           3 3 3 7 1387
          524287       (1) prime factors:  [prime]  524287
         1048575       (6) prime factors:           3 5 5 11 41 31
         2097151       (3) prime factors:           7 7 42799
         4194303       (4) prime factors:           3 23 89 683
         8388607       (2) prime factors:           47 178481
        16777215       (7) prime factors:           3 3 5 7 13 17 241
        33554431       (1) prime factors:  [prime]  33554431
        67108863       (2) prime factors:           3 22369621
       134217727       (2) prime factors:           7 19173961
       268435455       (5) prime factors:           3 5 29 113 5461
       536870911       (3) prime factors:           233 1103 2089
      1073741823       (5) prime factors:           3 3 7 11 1549411
      2147483647       (1) prime factors:  [prime]  2147483647
      4294967295       (5) prime factors:           3 5 17 257 65537
      8589934591       (4) prime factors:           7 23 89 599479
     17179869183       (3) prime factors:           3 43691 131071
     34359738367       (3) prime factors:           71 122921 3937
     68719476735       (7) prime factors:           3 3 3 5 7 13 5593771
    137438953471       (1) prime factors:  [prime]  137438953471
    274877906943       (2) prime factors:           3 91625968981
    549755813887       (2) prime factors:           7 78536544841
   1099511627775       (7) prime factors:           3 5 5 11 17 41 1912111
   2199023255551       (2) prime factors:           13367 164511353
   4398046511103       (5) prime factors:           3 3 7 7 9972894583
   8796093022207       (3) prime factors:           431 9719 2099863
  17592186044415       (6) prime factors:           3 5 23 89 683 838861
  35184372088831       (2) prime factors:           7 5026338869833
  70368744177663       (4) prime factors:           3 47 178481 2796203
 140737488355327       (2) prime factors:           2351 59862819377
 281474976710655       (8) prime factors:           3 3 5 7 13 17 257 15732721
 562949953421311       (1) prime factors:  [prime]  562949953421311
1125899906842623       (4) prime factors:           3 11 251 135928999981

                        11 primes found.

output   when using the input of:   1   50   1   2   +1

(essentially testing for   2n +1)

(Shown at three-quarter size.)

               3       (1) prime factors:  [prime]  3
               5       (1) prime factors:  [prime]  5
               9       (2) prime factors:           3 3
              17       (1) prime factors:  [prime]  17
              33       (2) prime factors:           3 11
              65       (2) prime factors:           5 13
             129       (2) prime factors:           3 43
             257       (1) prime factors:  [prime]  257
             513       (4) prime factors:           3 3 3 19
            1025       (3) prime factors:           5 5 41
            2049       (2) prime factors:           3 683
            4097       (2) prime factors:           17 241
            8193       (2) prime factors:           3 2731
           16385       (3) prime factors:           5 29 113
           32769       (4) prime factors:           3 3 11 331
           65537       (1) prime factors:  [prime]  65537
          131073       (2) prime factors:           3 43691
          262145       (3) prime factors:           5 13 4033
          524289       (2) prime factors:           3 174763
         1048577       (2) prime factors:           17 61681
         2097153       (3) prime factors:           3 3 233017
         4194305       (2) prime factors:           5 838861
         8388609       (2) prime factors:           3 2796203
        16777217       (2) prime factors:           257 65281
        33554433       (4) prime factors:           3 11 251 4051
        67108865       (4) prime factors:           5 53 1613 157
       134217729       (5) prime factors:           3 3 3 3 1657009
       268435457       (2) prime factors:           17 15790321
       536870913       (3) prime factors:           3 59 3033169
      1073741825       (5) prime factors:           5 5 13 41 80581
      2147483649       (2) prime factors:           3 715827883
      4294967297       (2) prime factors:           641 6700417
      8589934593       (4) prime factors:           3 3 683 1397419
     17179869185       (4) prime factors:           5 137 953 26317
     34359738369       (5) prime factors:           3 11 281 86171 43
     68719476737       (2) prime factors:           17 4042322161
    137438953473       (2) prime factors:           3 45812984491
    274877906945       (2) prime factors:           5 54975581389
    549755813889       (3) prime factors:           3 3 61083979321
   1099511627777       (2) prime factors:           257 4278255361
   2199023255553       (3) prime factors:           3 83 8831418697
   4398046511105       (5) prime factors:           5 13 29 113 20647621
   8796093022209       (2) prime factors:           3 2932031007403
  17592186044417       (3) prime factors:           17 353 2931542417
  35184372088833       (5) prime factors:           3 3 3 11 118465899289
  70368744177665       (4) prime factors:           5 1013 30269 458989
 140737488355329       (2) prime factors:           3 46912496118443
 281474976710657       (2) prime factors:           65537 4294901761
 562949953421313       (2) prime factors:           3 187649984473771
1125899906842625       (6) prime factors:           5 5 5 41 101 2175126601

                         5 primes found.

optimized more

This REXX version is about   20%   faster than the 1st REXX version when factoring one million numbers. <lang rexx>/*REXX pgm does prime decomposition of a range of positive integers (with a prime count)*/ numeric digits 1000 /*handle thousand digits for the powers*/ parse arg bot top step base add /*get optional arguments from the C.L. */ if bot== then do; bot=1; top=100; end /*no BOT given? Then use the default.*/ if top== then top=bot /* " TOP? " " " " " */ if step== then step= 1 /* " STEP? " " " " " */ if add == then add= -1 /* " ADD? " " " " " */ tell= top>0; top=abs(top) /*if TOP is negative, suppress displays*/ w=length(top) /*get maximum width for aligned display*/ if base\== then w=length(base**top) /*will be testing powers of two later? */ @.=left(, 7); @.0="{unity}"; @.1='[prime]' /*some literals: pad; prime (or not).*/ numeric digits max(9, w+1) /*maybe increase the digits precision. */

  1. =0 /*#: is the number of primes found. */
       do n=bot  to top  by step                /*process a single number  or  a range.*/
       ?=n;  if base\==  then ?=base**n + add /*should we perform a "Mercenne" test? */
       pf=factr(?);      f=words(pf)            /*get prime factors; number of factors.*/
       if f==1  then #=#+1                      /*Is N prime?  Then bump prime counter.*/
       if tell  then say right(?,w)   right('('f")",9)   'prime factors: '     @.f     pf
       end   /*n*/

say ps= 'primes'; if p==1 then ps= "prime" /*setup for proper English in sentence.*/ say right(#, w+9+1) ps 'found.' /*display the number of primes found. */ exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ factr: procedure; parse arg x 1 d,$ /*set X, D to argument 1; $ to null.*/ if x==1 then return /*handle the special case of X = 1. */

      do  while x// 2==0;  $=$  2;  x=x%2;  end /*append all the  2  factors of new  X.*/
      do  while x// 3==0;  $=$  3;  x=x%3;  end /*   "    "   "   3     "     "  "   " */
      do  while x// 5==0;  $=$  5;  x=x%5;  end /*   "    "   "   5     "     "  "   " */
      do  while x// 7==0;  $=$  7;  x=x%7;  end /*   "    "   "   7     "     "  "   " */
      do  while x//11==0;  $=$ 11;  x=x%11; end /*   "    "   "  11     "     "  "   " */    /* ◄■■■■ added.*/
      do  while x//13==0;  $=$ 13;  x=x%13; end /*   "    "   "  13     "     "  "   " */    /* ◄■■■■ added.*/
      do  while x//17==0;  $=$ 17;  x=x%17; end /*   "    "   "  17     "     "  "   " */    /* ◄■■■■ added.*/
      do  while x//19==0;  $=$ 19;  x=x%19; end /*   "    "   "  19     "     "  "   " */    /* ◄■■■■ added.*/
      do  while x//23==0;  $=$ 23;  x=x%23; end /*   "    "   "  23     "     "  "   " */    /* ◄■■■■ added.*/
                                                /*                                  ___*/

q=1; do while q<=x; q=q*4; end /*these two lines compute integer √ X */ r=0; do while q>1; q=q%4; _=d-r-q; r=r%2; if _>=0 then do; d=_; r=r+q; end; end

      do j=29  by 6  to r                       /*insure that  J  isn't divisible by 3.*/    /* ◄■■■■ changed.*/
      parse var j    -1  _                    /*obtain the last decimal digit of  J. */
      if _\==5  then  do  while x//j==0;  $=$ j;  x=x%j;  end     /*maybe reduce by J. */
      if _ ==3  then iterate                    /*Is next  Y  is divisible by 5?  Skip.*/
      y=j+2;          do  while x//y==0;  $=$ y;  x=x%y;  end     /*maybe reduce by J. */
      end   /*j*/
                                                /* [↓]  The $ list has a leading blank.*/

if x==1 then return $ /*Is residual=unity? Then don't append.*/

             return $ x                         /*return   $   with appended residual. */</lang>

output   is identical to the 1st REXX version.

Ring

<lang ring> prime = 18705 decomp(prime)

func decomp nr x = "" for i = 1 to nr

   if isPrime(i) and nr % i = 0
      x = x + string(i) + " * " ok
   if i = nr
      x2 = substr(x,1,(len(x)-2))
      see string(nr) + " = " + x2 + nl ok

next

func isPrime num

    if (num <= 1) return 0 ok
    if (num % 2 = 0) and num != 2 return 0 ok
    for i = 3 to floor(num / 2) -1 step 2
        if (num % i = 0) return 0 ok
    next
    return 1

</lang>

Ruby

Built in

<lang ruby>irb(main):001:0> require 'prime' => true irb(main):003:0> 2543821448263974486045199.prime_division => [[701, 1], [1123, 2], [2411, 1], [1092461, 2]]</lang>

Simple algorithm

<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. This may take several minutes to show that 2**31-1 is prime.

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

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

end</lang>

Output:
...
2**28-1 = 268435455 = 3 * 5 * 29 * 43 * 113 * 127
2**29-1 = 536870911 = 233 * 1103 * 2089
2**30-1 = 1073741823 = 3 * 3 * 7 * 11 * 31 * 151 * 331
2**31-1 = 2147483647 = 2147483647

Faster algorithm

<lang ruby># Get prime decomposition of integer _i_.

  1. This routine is more efficient than prime_factors,
  2. and quite similar to Integer#prime_division of MRI 1.9.

def prime_factors_faster(i)

 factors = []
 check = proc do |p|
   while(q, r = i.divmod(p)
         r.zero?)
     factors << p
     i = q
   end
 end
 check[2]
 check[3]
 p = 5
 while p * p <= i
   check[p]
   p += 2
   check[p]
   p += 4    # skip multiples of 2 and 3
 end
 factors << i if i > 1
 factors

end

  1. Example: Decompose all possible Mersenne primes up to 2**70-1.
  2. This may take several minutes to show that 2**61-1 is prime,
  3. but 2**62-1 and 2**67-1 are not prime.

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

 factors = prime_factors_faster(2**i-1)
 puts "2**#{i}-1 = #{2**i-1} = #{factors.join(' * ')}"

end</lang>

Output:
...
2**67-1 = 147573952589676412927 = 193707721 * 761838257287
2**68-1 = 295147905179352825855 = 3 * 5 * 137 * 953 * 26317 * 43691 * 131071
2**69-1 = 590295810358705651711 = 7 * 47 * 178481 * 10052678938039
2**70-1 = 1180591620717411303423 = 3 * 11 * 31 * 43 * 71 * 127 * 281 * 86171 * 122921

This benchmark compares the different implementations.

<lang ruby>require 'benchmark' require 'mathn' Benchmark.bm(24) do |x|

 [2**25 - 6, 2**35 - 7].each do |i|
   puts "#{i} = #{prime_factors_faster(i).join(' * ')}"
   x.report("  prime_factors") { prime_factors(i) }
   x.report("  prime_factors_faster") { prime_factors_faster(i) }
   x.report("  Integer#prime_division") { i.prime_division }
 end

end</lang>

With MRI 1.8, prime_factors is slow, Integer#prime_division is fast, and prime_factors_faster is very fast. With MRI 1.9, Integer#prime_division is also very fast.

Rust

Rust's largest built-in integer type is u128 (128-bit unsigned integer) which is pretty large, but not unlimited. The solution therefore uses external crates for big integers. Add the dependencies in Cargo.toml:

[package]
name = "prime_decomposition"
version = "0.1.1"
edition = "2018"

[dependencies]
num-bigint = "0.3.0"
num-traits = "0.2.12"

The implementation:

<lang Rust>use num_bigint::BigUint; use num_traits::{One, Zero}; use std::fmt::{Display, Formatter};

  1. [derive(Clone, Debug)]

pub struct Factors {

   pub number: BigUint,
   pub result: Vec<BigUint>,

}

impl Factors {

   pub fn of(number: BigUint) -> Factors {
       let mut factors = Self {
           number: number.clone(),
           result: Vec::new(),
       };
       let big_2 = BigUint::from(2u8);
       let big_4 = BigUint::from(4u8);
       factors.check(&big_2);
       factors.check(&BigUint::from(3u8));
       let mut divisor = BigUint::from(5u8);
       while &divisor * &divisor <= factors.number {
           factors.check(&divisor);
           divisor += &big_2;
           factors.check(&divisor);
           divisor += &big_4;
       }
       if factors.number > BigUint::one() {
           factors.result.push(factors.number);
       }
       factors.number = number; // Restore the number
       factors
   }
   pub fn is_prime(&self) -> bool {
       self.result.len() == 1
   }
   fn check(&mut self, divisor: &BigUint) {
       while (&self.number % divisor).is_zero() {
           self.result.push(divisor.clone());
           self.number /= divisor;
       }
   }

}

impl Display for Factors {

   fn fmt(&self, f: &mut Formatter) -> std::fmt::Result {
       let mut iter = self.result.iter();
       match iter.next() {
           None => write!(f, "[]"),
           Some(first) => {
               write!(f, "[{}", first)?;
               for next in iter {
                   write!(f, ", {}", next)?;
               }
               write!(f, "]")
           }
       }
   }

}

fn print_factors(number: BigUint) {

   let factors = Factors::of(number);
   if factors.is_prime() {
       println!("{} -> {} (prime)", factors.number, factors);
   } else {
       println!("{} -> {}", factors.number, factors);
   }

}

fn main() {

   print_factors(24u32.into());
   print_factors(32u32.into());
   print_factors(37u32.into());
   // Find Mersenne primes
   for n in 2..70 {
       print!("2**{} - 1: ", n);
       print_factors((BigUint::from(2u8) << n) - BigUint::one());
   }

} </lang>

S-BASIC

<lang S-BASIC> rem - return p mod q function mod(p, q = integer) = integer end = p - q * (p/q)

dim integer factors(16) rem log2(maxint) is sufficiently large

comment

 Find the prime factors of n and store in global array factors
 (arrays cannot be passed as parameters) and return the number
 found. If n is prime, it will be stored as the only factor.

end function primefactors(n = integer) = integer

 var p, count = integer
 p = 2
 count = 1
 while n >= (p * p) do
    begin
      if mod(n, p) = 0 then
        begin
          factors(count) = p
          count = count + 1
          n = n / p
        end
      else
        p = p + 1
    end
  factors(count) = n

end = count

rem -- exercise the routine by checking odd numbers from 77 to 99

var i, k, nfound = integer

for i = 77 to 99 step 2

 nfound = primefactors(i)
 print i;"; ";
 for k = 1 to nfound
   print factors(k);
 next k
 print

next i

end </lang>

Output:
77: 7 11
79: 79
81: 3 3 3 3
83: 83
85: 5 17
87: 3 29
89: 89
91: 7 13
93: 3 31
95: 5 19
97: 97
99: 3 3 11

Scala

Library: Scala

<lang Scala>import annotation.tailrec import collection.parallel.mutable.ParSeq

object PrimeFactors extends App {

 def factorize(n: Long): List[Long] = {
   @tailrec
   def factors(tuple: (Long, Long, List[Long], Int)): List[Long] = {
     tuple match {
       case (1, _, acc, _)                 => acc
       case (n, k, acc, _) if (n % k == 0) => factors((n / k, k, acc ++ ParSeq(k), Math.sqrt(n / k).toInt))
       case (n, k, acc, sqr) if (k < sqr)  => factors(n, k + 1, acc, sqr)
       case (n, k, acc, sqr) if (k >= sqr) => factors((1, k, acc ++ ParSeq(n), 0))
     }
   }
   factors((n, 2, List[Long](), Math.sqrt(n).toInt))
 }
 def mersenne(p: Int): BigInt = (BigInt(2) pow p) - 1
 def sieve(nums: Stream[Int]): Stream[Int] =
   Stream.cons(nums.head, sieve((nums.tail) filter (_ % nums.head != 0)))
 // An infinite stream of primes, lazy evaluation and memo-ized
 val oddPrimes = sieve(Stream.from(3, 2))
 def primes = sieve(2 #:: oddPrimes)
 oddPrimes takeWhile (_ <= 59) foreach { p =>
   { // Needs some intermediate results for nice formatting
     val numM = s"M${p}"
     val nMersenne = mersenne(p).toLong
     val lit = f"${nMersenne}%30d"
     val datum = System.nanoTime
     val result = factorize(nMersenne)
     val mSec = ((System.nanoTime - datum) / 1.0e+6).round
     def decStr = { if (lit.length > 30) f"(M has ${lit.length}%3d dec)" else "" }
     def sPrime = { if (result.isEmpty) " is a prime number." else "" }
     println(
       f"$numM%4s = 2^$p%03d - 1 = ${lit}%s${sPrime} ($mSec%,4d msec) composed of ${result.mkString(" × ")}")
   }
 }

}</lang>

Output:
  M3 = 2^003 - 1 =                              7 (  23 msec) composed of 7
  M5 = 2^005 - 1 =                             31 (   0 msec) composed of 31
  M7 = 2^007 - 1 =                            127 (   0 msec) composed of 127
 M11 = 2^011 - 1 =                           2047 (   0 msec) composed of 23 × 89
 M13 = 2^013 - 1 =                           8191 (   0 msec) composed of 8191
 M17 = 2^017 - 1 =                         131071 (   1 msec) composed of 131071
 M19 = 2^019 - 1 =                         524287 (   1 msec) composed of 524287
 M23 = 2^023 - 1 =                        8388607 (   1 msec) composed of 47 × 178481
 M29 = 2^029 - 1 =                      536870911 (   2 msec) composed of 233 × 1103 × 2089
 M31 = 2^031 - 1 =                     2147483647 (  39 msec) composed of 2147483647
 M37 = 2^037 - 1 =                   137438953471 (   8 msec) composed of 223 × 616318177
 M41 = 2^041 - 1 =                  2199023255551 (   2 msec) composed of 13367 × 164511353
 M43 = 2^043 - 1 =                  8796093022207 (   2 msec) composed of 431 × 9719 × 2099863
 M47 = 2^047 - 1 =                140737488355327 (   2 msec) composed of 2351 × 4513 × 13264529
 M53 = 2^053 - 1 =               9007199254740991 (   7 msec) composed of 6361 × 69431 × 20394401
 M59 = 2^059 - 1 =             576460752303423487 ( 152 msec) composed of 179951 × 3203431780337

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>

Output:

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

Alternatively, Scala LazyLists and Iterators support quite elegant one-line encodings of iterative/recursive algorithms, allowing us to to define the prime factorization like so: <lang scala>import spire.math.SafeLong import spire.implicits._ def pFactors(num: SafeLong): Vector[SafeLong] = Iterator.iterate((Vector[SafeLong](), num, SafeLong(2))){case (ac, n, f) => if(n%f == 0) (ac :+ f, n/f, f) else (ac, n, f + 1)}.dropWhile(_._2 != 1).next._1</lang>

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]

SequenceL

Recursive Using isPrime

<lang sequencel>isPrime(n) := n = 2 or (n > 1 and none(n mod ([2]++((1...floor(sqrt(n)/2))*2+1)) = 0));

primeFactorization(num) := primeFactorizationHelp(num, []);

primeFactorizationHelp(num, current(1)) :=

    let
       primeFactors[i] := i when num mod i = 0 and isPrime(i) foreach i within 2 ... num;
    in
           current when size(primeFactors) = 0
       else
           primeFactorizationHelp(num / product(primeFactors), current ++ primeFactors);</lang>

Using isPrime Based On: [2]

Recursive Trial Division

<lang sequencel>primeFactorization(num) := primeFactorizationHelp(num, 2, []);

primeFactorizationHelp(num, divisor, factors(1)) :=

       factors when num <= 1
   else
       primeFactorizationHelp(num, divisor + 1, factors) when num mod divisor /= 0
   else
       primeFactorizationHelp(num / divisor, divisor, factors ++ [divisor]);</lang>

Sidef

Built-in: <lang ruby>say factor(536870911) #=> [233, 1103, 2089] say factor_exp(536870911) #=> [[233, 1], [1103, 1], [2089, 1]]</lang>

Trial division: <lang ruby>func prime_factors(n) {

   return [] if (n < 1)
   gather {
       while (!(n & 1)) {
           n >>= 1
           take(2)
       }
       var p = 3
       while ((n > 1) && (p*p <= n)) {
           while (n %% p) {
               n //= p
               take(p)
           }
           p += 2
       }
       take(n) if (n > 1)
   }

}</lang>

Calling the function: <lang ruby>say prime_factors(536870911) #=> [233, 1103, 2089]</lang>

Simula

Simula has no built-in function to test for prime numbers.
Code for class bignum can be found here: https://rosettacode.org/wiki/Pi#Simula <lang simula> EXTERNAL CLASS BIGNUM; BIGNUM BEGIN

   CLASS TEXTLIST;
   BEGIN
       CLASS TEXTARRAY(N); INTEGER N;
       BEGIN
           TEXT ARRAY DATA(1:N);
       END TEXTARRAY;
       PROCEDURE EXPAND(N); INTEGER N;
       BEGIN
           REF(TEXTARRAY) NEWARR;
           INTEGER I;
           NEWARR :- NEW TEXTARRAY(N);
           FOR I := 1 STEP 1 UNTIL SIZE DO BEGIN
               NEWARR.DATA(I) :- ARR.DATA(I);
           END;
           ARR :- NEWARR;
       END EXPAND;
       PROCEDURE APPEND(T); TEXT T;
       BEGIN
           IF SIZE = ARR.N THEN
               EXPAND(2*ARR.N);
           SIZE := SIZE+1;
           ARR.DATA(SIZE) :- T;
       END EXPAND;
       TEXT PROCEDURE GET(I); INTEGER I;
           GET :- ARR.DATA(I);
       REF(TEXTARRAY) ARR;
       INTEGER SIZE;
       EXPAND(20);
   END TEXTLIST;
   REF(TEXTLIST) PROCEDURE PRIME_FACTORS(N); TEXT N;
   BEGIN
       REF(TEXTLIST) FACTORS;
       REF(DIVMOD) DM;
       TEXT P;
       FACTORS :- NEW TEXTLIST;
       IF TCMP(N, "1") < 0 THEN
           GOTO RETURN;
       P :- "2";
       FOR DM :- TDIVMOD(N,P) WHILE TISZERO(DM.MOD) DO BEGIN
           N :- DM.DIV;
           FACTORS.APPEND(P);
       END;
       P :- "3";
       WHILE TCMP(N,"1") > 0 AND THEN TCMP(TMUL(P,P),N) <= 0 DO BEGIN
           FOR DM :- TDIVMOD(N, P) WHILE TISZERO(DM.MOD) DO BEGIN
               N :- DM.DIV;
               FACTORS.APPEND(P);
           END;
           P :- TADD(P,"2");
       END;
       IF TCMP(N,"1") > 0 THEN
           FACTORS.APPEND(N);
   RETURN:
       PRIME_FACTORS :- FACTORS;
   END PRIME_FACTORS;
   REF(TEXTLIST) FACTORS;
   TEXT INP;
   INTEGER I;
   FOR INP :- "536870911", "6768768", "1957", "64865899369365843" DO BEGIN
       FACTORS :- PRIME_FACTORS(INP);
       OUTTEXT("PRIME FACTORS OF ");
       OUTTEXT(INP);
       OUTTEXT(" => [");
       FOR I := 1 STEP 1 UNTIL FACTORS.SIZE DO BEGIN
           IF I > 1 THEN
               OUTTEXT(", ");
           OUTTEXT(FACTORS.GET(I));
       END;
       OUTTEXT("]");
       OUTIMAGE;
   END;

END; </lang>

Output:
PRIME FACTORS OF 536870911 => [233, 1103, 2089]
PRIME FACTORS OF 6768768 => [2, 2, 2, 2, 2, 2, 2, 3, 17627]
PRIME FACTORS OF 1957 => [19, 103]
PRIME FACTORS OF 64865899369365843 => [3, 7, 397, 276229, 28166791]

5320 garbage collection(s) in 1.9 seconds.

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>

SPAD

<lang SPAD>

(1) -> factor 102400

        12 2
  (1)  2  5
                                                     Type: Factored(Integer)

(2) -> factor 23193931893819371

  (2)  83 3469 71341 1129153
                                                     Type: Factored(Integer)

</lang>

Domain:Factored(R)

Standard ML

Trial division <lang SML> val factor = fn n :IntInf.int => let

val unfactored  = fn (u,_,_)   => u
val factors     = fn (_,f,_)   => f
val try         = fn (_,_,i)   => i
fun getresult t = unfactored t::(factors t)

fun until done change x =
   if done x
     then    getresult x
     else    until done change (change x);       (* iteration *)
fun lastprime t = unfactored t  <  (try t)*(try t)
fun trymore t   = if unfactored t mod (try t) = 0
          then (unfactored t div (try t) , try t::(factors t) , try t    )
      else (unfactored t             ,         factors t  , try t + 1)

in

until lastprime trymore (n,[],2)

end; </lang>

- Array.fromList(factor 122489234920000001278234798233);;
val it = fromList[658601127263, 41259943, 34942753, 43, 3]: IntInf.int array

Stata

The following Mata function will factor any representable positive integer (that is, between 1 and 2^53).

<lang stata>function factor(n_) {

   n = n_
   a = J(0,2,.)
   if (n<2) {
       return(a)
   }
   else if (n<4) {
       return((n,1))
   }
   else {
       if (mod(n,2)==0) {
           for (i=0; mod(n,2)==0; i++) n = floor(n/2)
           a = a\(2,i)
       }
           
       for (k=3; k*k<=n; k=k+2) {
           if (mod(n,k)==0) {
               for (i=0; mod(n,k)==0; i++) n = floor(n/k)
               a = a\(k,i)
           }
       }
       
       if (n>1) a = a\(n,1)
       return(a)
   }

}</lang>

Swift

Translation of: Python

Uses the sieve of Eratosthenes. This is generic on any type that conforms to BinaryInteger. So in theory any BigInteger library should work with it.

<lang swift>func primeDecomposition<T: BinaryInteger>(of n: T) -> [T] {

 guard n > 2 else { return [] }
 func step(_ x: T) -> T {
   return 1 + (x << 2) - ((x >> 1) << 1)
 }
 let maxQ = T(Double(n).squareRoot())
 var d: T = 1
 var q: T = n % 2 == 0 ? 2 : 3
 while q <= maxQ && n % q != 0 {
   q = step(d)
   d += 1
 }
 return q <= maxQ ? [q] + primeDecomposition(of: n / q) : [n]

}

for prime in Eratosthenes(upTo: 60) {

 let m = Int(pow(2, Double(prime))) - 1
 let decom = primeDecomposition(of: m)
 print("2^\(prime) - 1 = \(m) => \(decom)")

}</lang>

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

Tcl

<lang tcl>proc factors {x} {

   # list the prime factors of x in ascending order
   set result [list]
   while {$x % 2 == 0} {
       lappend result 2
       set x [expr {$x / 2}]
   }
   for {set i 3} {$i*$i <= $x} {incr i 2} {
       while {$x % $i == 0} {
           lappend result $i
           set x [expr {$x / $i}]
       }
   }
   if {$x != 1} {lappend result $x}
   return $result

} </lang> Testing <lang tcl>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 primes [factors $n]} 1} tm
   puts [format "2**%02d-1 = %-18s = %-22s => %s" $m $n [join $primes *] $tm]

}</lang>

Output:
2**02-1 = 3                  = 3                      => 184 microseconds per iteration
2**03-1 = 7                  = 7                      => 8 microseconds per iteration
2**05-1 = 31                 = 31                     => 8 microseconds per iteration
2**07-1 = 127                = 127                    => 23 microseconds per iteration
2**11-1 = 2047               = 23*89                  => 12 microseconds per iteration
2**13-1 = 8191               = 8191                   => 22 microseconds per iteration
2**17-1 = 131071             = 131071                 => 69 microseconds per iteration
2**19-1 = 524287             = 524287                 => 131 microseconds per iteration
2**23-1 = 8388607            = 47*178481              => 81 microseconds per iteration
2**29-1 = 536870911          = 233*1103*2089          => 199 microseconds per iteration
2**31-1 = 2147483647         = 2147483647             => 9509 microseconds per iteration
2**37-1 = 137438953471       = 223*616318177          => 4377 microseconds per iteration
2**41-1 = 2199023255551      = 13367*164511353        => 2389 microseconds per iteration
2**43-1 = 8796093022207      = 431*9719*2099863       => 1711 microseconds per iteration
2**47-1 = 140737488355327    = 2351*4513*13264529     => 802 microseconds per iteration
2**53-1 = 9007199254740991   = 6361*69431*20394401    => 13109 microseconds per iteration
2**59-1 = 576460752303423487 = 179951*3203431780337   => 316009 microseconds per iteration

TI-83 BASIC

<lang ti83b>::prgmPREMIER Disp "FACTEURS PREMIER" Prompt N If N<1:Stop ClrList L1�,L2 0→K iPart(√(N))→L N→M For(I,2,L) 0→J While fPart(M/I)=0 J+1→J M/I→M End If J≠0 Then K+1→K I→L�1(K) J→L2(K) I→Z:prgmVSTR " "+Str0→Str1 If J≠1 Then J→Z:prgmVSTR Str1+"^"+Str0→Str1 End Disp Str1 End If M=1:Stop End If M≠1 Then If M≠N Then M→Z:prgmVSTR " "+Str0→Str1 Disp Str1 Else Disp "PREMIER" End End

prgmVSTR

{Z,Z}→L5 {1,2}→L6 LinReg(ax+b)L6,L5,Y�₀ Equ►String(Y₀,Str0) length(Str0)→O sub(Str0,4,O-3)→Str0 ClrList L5,L6 DelVar Y�</lang>

Output:
FACTEURS PREMIER
N=?1047552
   2^10
      3
     11
     31

Tiny BASIC

<lang Tiny BASIC>10 PRINT "Enter a number." 20 INPUT N 25 PRINT "------" 30 IF N<0 THEN LET N = -N 40 IF N<2 THEN END 50 LET I = 2 60 IF I*I > N THEN GOTO 200 70 IF (N/I)*I = N THEN GOTO 300 80 LET I = I + 1 90 GOTO 60 200 PRINT N 210 END 300 LET N = N / I 310 PRINT I 320 GOTO 50</lang>

Output:
Enter a number.
32
------
2
2
2
2
2

Enter a number.
2520
------
2
2
2
3
3
5
7

Enter a number.
13
------
13

TXR

Translation of: Common Lisp

<lang txr>@(next :args) @(do

 (defun factor (n)
   (if (> n 1)
     (for ((max-d (isqrt n))
           (d 2))
          ()
          ((inc d (if (evenp d) 1 2)))
       (cond ((> d max-d) (return (list n)))
             ((zerop (mod n d)) 
              (return (cons d (factor (trunc n d))))))))))

@{num /[0-9]+/} @(bind factors @(factor (int-str num 10))) @(output) @num -> {@(rep)@factors, @(last)@factors@(end)} @(end)</lang>

Output:
$ txr factor.txr 1139423842450982345
1139423842450982345 -> {5, 19, 37, 12782467, 25359769}
$ txr factor.txr 1
1 -> {}
$ txr factor.txr 2
2 -> {2}
$ txr factor.txr 3
3 -> {3}
$ txr factor.txr 2
2 -> {2}
$ txr factor.txr 3
3 -> {3}
$ txr factor.txr 4
4 -> {2, 2}
$ txr factor.txr 5
5 -> {5}
$ txr factor.txr 6
6 -> {2, 3}

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]

VBScript

<lang vb>Function PrimeFactors(n)

   arrP = Split(ListPrimes(n)," ")
   divnum = n
   Do Until divnum = 1
       'The -1 is to account for the null element of arrP
       For i = 0 To UBound(arrP)-1
           If divnum = 1 Then
               Exit For
           ElseIf divnum Mod arrP(i) = 0 Then
               divnum = divnum/arrP(i)
               PrimeFactors = PrimeFactors & arrP(i) & " "
           End If
       Next
   Loop

End Function

Function IsPrime(n)

   If n = 2 Then
       IsPrime = True
   ElseIf n <= 1 Or n Mod 2 = 0 Then
       IsPrime = False
   Else
       IsPrime = True
       For i = 3 To Int(Sqr(n)) Step 2
           If n Mod i = 0 Then
               IsPrime = False
               Exit For
           End If
       Next
   End If

End Function

Function ListPrimes(n)

   ListPrimes = ""
   For i = 1 To n
       If IsPrime(i) Then
           ListPrimes = ListPrimes & i & " "
       End If
   Next

End Function

WScript.StdOut.Write PrimeFactors(CInt(WScript.Arguments(0))) WScript.StdOut.WriteLine</lang>

Output:
C:\>cscript /nologo primefactors.vbs 12
2 3 2

C:\>cscript /nologo primefactors.vbs 50
2 5 5

Wren

Library: Wren-big
Library: Wren-fmt

The examples are borrowed from the Go solution. <lang ecmascript>import "/big" for BigInt import "/fmt" for Fmt

var vals = [1 << 31, 1234567, 333333, 987653, 2 * 3 * 5 * 7 * 11 * 13 * 17] for (val in vals) {

   Fmt.print("$10d -> $n", val, BigInt.primeFactors(val))

}</lang>

Output:
2147483648 -> [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
   1234567 -> [127, 9721]
    333333 -> [3, 3, 7, 11, 13, 37]
    987653 -> [29, 34057]
    510510 -> [2, 3, 5, 7, 11, 13, 17]

XSLT

Let's assume that in XSLT the application of a template is similar to the invocation of a function. So when the following template <lang xml><xsl:stylesheet xmlns:xsl="http://www.w3.org/1999/XSL/Transform" version="1.0">

   <xsl:template match="/numbers">
       <html>
           <body>
    <xsl:apply-templates />
           </body>
       </html>
   </xsl:template>
   <xsl:template match="number">
  • Number: <xsl:apply-templates mode="value" /> Factors: <xsl:apply-templates mode="factors" />
  • </xsl:template> <xsl:template match="value" mode="value"> <xsl:apply-templates /> </xsl:template> <xsl:template match="value" mode="factors"> <xsl:call-template name="generate"> <xsl:with-param name="number" select="number(current())" /> <xsl:with-param name="candidate" select="number(2)" /> </xsl:call-template> </xsl:template> <xsl:template name="generate"> <xsl:param name="number" /> <xsl:param name="candidate" /> <xsl:choose> <xsl:when test="$number = 1"></xsl:when> <xsl:when test="$candidate * $candidate > $number"> <xsl:value-of select="$number" /> </xsl:when> <xsl:when test="$number mod $candidate = 0"> <xsl:value-of select="$candidate" /> <xsl:text> </xsl:text> <xsl:call-template name="generate"> <xsl:with-param name="number" select="$number div $candidate" /> <xsl:with-param name="candidate" select="$candidate" /> </xsl:call-template> </xsl:when> <xsl:otherwise> <xsl:choose> <xsl:when test="$candidate = 2"> <xsl:call-template name="generate"> <xsl:with-param name="number" select="$number" /> <xsl:with-param name="candidate" select="$candidate + 1" /> </xsl:call-template> </xsl:when> <xsl:otherwise> <xsl:call-template name="generate"> <xsl:with-param name="number" select="$number" /> <xsl:with-param name="candidate" select="$candidate + 2" /> </xsl:call-template> </xsl:otherwise> </xsl:choose> </xsl:otherwise> </xsl:choose> </xsl:template> </xsl:stylesheet></lang> is applied against the document <lang xml><numbers> <number><value>1</value></number> <number><value>2</value></number> <number><value>4</value></number> <number><value>8</value></number> <number><value>9</value></number> <number><value>255</value></number> </numbers></lang> then the output contains the prime decomposition of each number: <lang html><html> <body>

    • Number: 1 Factors:
    • Number: 2 Factors: 2
    • Number: 4 Factors: 2 2
    • Number: 8 Factors: 2 2 2
    • Number: 9 Factors: 3 3
    • Number: 255 Factors: 3 5 17

    </body> </html></lang>

    zkl

    With 64 bit ints: <lang zkl>fcn primeFactors(n){ // Return a list of factors of n

      acc:=fcn(n,k,acc,maxD){  // k is 2,3,5,7,9,... not optimum
         if(n==1 or k>maxD) acc.close();
         else{
        q,r:=n.divr(k);   // divr-->(quotient,remainder)
        if(r==0) return(self.fcn(q,k,acc.write(k),q.toFloat().sqrt()));
        return(self.fcn(n,k+1+k.isOdd,acc,maxD))
         }
      }(n,2,Sink(List),n.toFloat().sqrt());
      m:=acc.reduce('*,1);      // mulitply factors
      if(n!=m) acc.append(n/m); // opps, missed last factor
      else acc;
    

    }</lang> <lang zkl>foreach n in (T(5,12, 2147483648, 2199023255551, 8796093022207,

       9007199254740991, 576460752303423487)){
      println(n,": ",primeFactors(n).concat(", ")) 
    

    }</lang>

    Output:
    5: 5
    12: 2, 2, 3
    2147483648: 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2
    2199023255551: 13367, 164511353
    8796093022207: 431, 9719, 2099863
    9007199254740991: 6361, 69431, 20394401
    576460752303423487: 179951, 3203431780337
    

    Unfortunately, big ints (GMP) don't have (quite) the same interface as ints (since there is no big float, BI.toFloat() truncates to a double so BI.toFloat().sqrt() is wrong). So mostly duplicate code is needed: <lang zkl>fcn factorsBI(n){ // Return a list of factors of n

      acc:=fcn(n,k,acc,maxD){  // k is 2,3,5,7,9,... not optimum
         if(n==1 or k>maxD) acc.close();
         else{
        q,r:=n.div2(k);   // divr-->(quotient,remainder)
        if(r==0) return(self.fcn(q,k,acc.write(k),q.root(2)));
        return(self.fcn(n,k+1+k.isOdd,acc,maxD))
         }
      }(n,2,Sink(List),n.root(2));
      m:=acc.reduce('*,BN(1));  // mulitply factors
      if(n!=m) acc.append(n/m); // opps, missed last factor
      else acc;
    

    }</lang> <lang zkl>var BN=Import("zklBigNum"); foreach n in (T(BN("12"),

       BN("340282366920938463463374607431768211455"))){
      println(n,": ",factorsBI(n).concat(", ")) 
    

    }</lang>

    Output:
    12: 2, 2, 3
    340282366920938463463374607431768211455: 3, 5, 17, 257, 641, 65537, 274177, 6700417, 67280421310721