Prime decomposition

You are encouraged to 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 n {\displaystyle n} greater than 1.

If your language does not have an isPrime-like function available, you may assume that you have a function which determines whether a number is prime (note its name before your code).

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

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





Related tasks







For maximum compatibility, this program uses only the basic instruction set.

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

Output:

1047552 : 2 2 2 2 2 2 2 2 2 2 3 11 31

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



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

"

szMessEndPgm: .asciz "Program normal end.

"

szMessNotPrime: .asciz "Not prime.

"

szMessPrime: .asciz "Prime

"

szCarriageReturn: .asciz "

"

szSpaces: .asciz " "

szMessNumber: .asciz " The factors of @ are :

"

/*******************************************/

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





Output:

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

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 .

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

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:

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;

The function Decompose first estimates the maximal result length as log 2 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:

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;

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

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;

Output:

2 2 3

Translation of: Python

Works with: ALGOL 68 version Revision 1 - no extensions to language used

#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 ;



#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

) ;



# 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 : (

# FOR LINT m IN # gen primes ( # ) DO ( #

## ( 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

# OD # ) ) ;

done : EMPTY

)

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

(decomposition of 12):- note: This specimen retains the original Python coding style.

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

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

loop [filter 2..60 => isPrime] @(num){

n: 2^num-1

print "2^" + num + "-1 = " + n + " => prime decomposition: " + [primeFactors n]

}



Output:

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

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

}

}

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

# 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

}



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

{ print pfac ( $1 ) }



Output:

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

entering input on stdin:

Unfortunately Batch does'nt have a BigNum library so the maximum number that can be decomposed is 2^31-1



@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



Works with: befungee

Handles safely integers only up to 250 (or ones which don't have prime divisors greater than 250).

& 211p > : 1 - #v_ 25*, @ > 11g:. / v

> : 11g %!|

> 11g 1+ 11p v

^ <



blsq ) 12fC

{2 2 3}



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.

#include <inttypes.h>

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <assert.h>



typedef uint32_t pint ;

typedef uint64_t xint ;

typedef unsigned int uint ;

#define PRIuPINT PRIu32 /* printf macro for pint */

#define PRIuXINT PRIu64 /* printf macro for xint */

#define MAX_FACTORS 63 /* because 2^64 is too large for xint */



uint8_t * pbits ;



#define MAX_PRIME (~(pint)0)

#define MAX_PRIME_SQ 65535U

#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 "

" , 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 ( '

' ) ;

}



return 0 ;

}

Using GNU Compiler Collection gcc extensions [ edit ]

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.

#include <limits.h>

#include <stdio.h>

#include <math.h>



typedef enum { false = 0 , true = 1 } bool ;

const int max_lint = LONG_MAX ;



typedef long long int lint ;

#assert sizeof_long_long_int (LONG_MAX>=8) /* XXX */



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

#define FOR(i,iterator) auto bool lambda(i); yield_init = (void *)λ iterator; bool lambda(i)

#define DO {

#define YIELD(x) if(!yield(x))return

#define BREAK return false

#define CONTINUE return true

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

#define YIELDS(type) bool (*yield)(type) = yield_init



typedef unsigned int bits ;

#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 ( "

" , m ) ;

if ( m >= 59 ) BREAK ;

CONTINUE ;

}

}

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.

cpp

#include <stdio.h>

#include <stdlib.h>

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

{

# define NCACHE 256

# 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 " : "

" ) ;

}



return 0 ;

}

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 ;

}

}

}

Simple trial division [ edit ]

To understand what was going on with the above code, pass it throughand 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):

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.

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 ;

}

}

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

#include <iostream>

#include <gmpxx.h>



// This function template works for any type representing integers or

// nonnegative integers, and has the standard operator overloads for

// arithmetic and comparison operators, as well as explicit conversion

// from int.

//

// OutputIterator must be an output iterator with value_type Integer.

// It receives the prime factors.

template < typename Integer, typename OutputIterator >

void decompose ( Integer n, OutputIterator out )

{

Integer i ( 2 ) ;



while ( n ! = 1 )

{

while ( n % i == Integer ( 0 ) )

{

* out ++ = i ;

n / = i ;

}

++ i ;

}

}



// this is an output iterator similar to std::ostream_iterator, except

// that it outputs the separation string *before* the value, but not

// before the first value (i.e. it produces an infix notation).

template < typename T > class infix_ostream_iterator :

public std :: iterator < T, std :: output_iterator_tag >

{

class Proxy ;

friend class Proxy ;

class Proxy

{

public :

Proxy ( infix_ostream_iterator & iter ) : iterator ( iter ) { }

Proxy & operator = ( T const & value )

{

if ( ! iterator. first )

{

iterator. stream << iterator. infix ;

}

iterator. stream << value ;

}

private :

infix_ostream_iterator & iterator ;

} ;

public :

infix_ostream_iterator ( std :: ostream & os, char const * inf ) :

stream ( os ) ,

first ( true ) ,

infix ( inf )

{

}

infix_ostream_iterator & operator ++ ( ) { first = false ; return * this ; }

infix_ostream_iterator operator ++ ( int )

{

infix_ostream_iterator prev ( * this ) ;

++* this ;

return prev ;

}

Proxy operator * ( ) { return Proxy ( * this ) ; }

private :

std :: ostream & stream ;

bool first ;

char const * infix ;

} ;



int main ( )

{

std :: cout << "please enter a positive number: " ;

mpz_class number ;

std :: cin >> number ;



if ( number <= 0 )

std :: cout << "this number is not positive!

;" ;

else

{

std :: cout << "decomposition: " ;

decompose ( number, infix_ostream_iterator < mpz_class > ( std :: cout , " * " ) ) ;

std :: cout << "

" ;

}

}

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

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:

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

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

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

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

%)" , 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 ;

}

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

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.

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

}

The built-in prime-factors function performs the task.

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

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

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

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

factor(5000)

Output:

2x2x2x5x5x5x5

Translation of: F#

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

Output:

[71,839,1471,6857]

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)

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

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



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



This version is a translation from Commodore BASIC program.



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



நிரல்பாகம் பகாஎண்ணா(எண்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(உள்ளீடு("உங்களுக்குப் பிடித்த ஓர் எண்ணைத் தாருங்கள்: "))



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



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



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



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 )

Output:

[71; 839; 1471; 6857]

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.

USING: io kernel math math.parser math.primes.factors sequences ;



27720 factors

[ number>string ] map

" " join print ;

: 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 . ;

Works with: Fortran version 90 and later

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

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

' 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

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

println[factor[2^508-1]]

Output:

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

(total process time including JVM startup = 1.515 s):

Note that this means 31 * 51 * ...

Built-in function :

FactorsInt(2^67-1);

# [ 193707721, 761838257287 ]

Or using the FactInt package :

FactInt(2^67-1);

# [ [ 193707721, 761838257287 ], [ ] ]

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

}

}

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]

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.

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

}

Test #1:

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

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:

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

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.

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

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 = [ ]

but it is not very efficient, to put it mildly. Inlining and fusing gets us the progressively more optimized

import Data . Maybe ( listToMaybe )

import Data . List ( unfoldr )



factorize :: Integer -> [ Integer ]

factorize n

= unfoldr (

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

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

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

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]

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

Sample Output showing factors of a large integer:

8796093022207=431*9719*2099863

q:

Example use:

q: 3684

2 2 3 307

and, more elaborately:

_1 + 2 ^ 128x

340282366920938463463374607431768211455

q: _1 + 2 ^ 128x

3 5 17 257 641 65537 274177 6700417 67280421310721

*/ q: _1 + 2 ^ 128x

340282366920938463463374607431768211455

Works with: Java version 1.5+

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

public boolean prime ( BigInteger i ) ;

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

public static List < BigInteger > primeFactorBig ( BigInteger a ) {

List < BigInteger > ans = new LinkedList < BigInteger > ( ) ;

//loop until we test the number itself or the number is 1

for ( BigInteger i = BigInteger . valueOf ( 2 ) ; i. compareTo ( a ) <= 0 && ! a. equals ( BigInteger . ONE ) ;

i = i. add ( BigInteger . ONE ) ) {

while ( a. remainder ( i ) . equals ( BigInteger . ZERO ) && prime ( i ) ) { //if we have a prime factor

ans. add ( i ) ; //put it in the list

a = a. divide ( i ) ; //factor it out of the number

}

}

return ans ;

}

Alternate version, optimised to be faster.

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 ;

}

Another alternate version designed to make fewer modular calculations:



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 ;

}



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

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 ;

}

This code uses the BigInteger Library jsbn and jsbn2

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

}

}

Without any library.

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 ;

}

TDD using Jasmine

PrimeFactors.js

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 ;

} ;



SpecPrimeFactors.js (with tag for Chutzpah)

/// <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 ] ) ;

} ) ;

} ) ;



Works with: jq version 1.4

"factors" as defined below emits a stream of all the prime factors of the input integer. The implementation is compact, fast and highly 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.

The caveat is that the program uses jq's builtin arithmetic operations. Since jq currently uses IEEE 754 64-bit numbers, 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.

def factors:

. as $in

| [2, $in, false]

| recurse( .[0] as $p |.[1] as $q | .[2] as $valid | .[3] as $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 ;

Examples:

[9007199254740992 | factors] | length

#=> 53



# 2**29-1 = 536870911

[ 536870911 | factors ]



#=> [233,1103,2089]

using package Primes.jl:



julia> Pkg.add("Primes")

julia> factor(8796093022207)

[9719=>1,431=>1,2099863=>1]



(The factor function returns a dictionary whose keys are the factors and whose values are the multiplicity of each factor.)

// 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)}" )

}

}

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]



{ 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 ]





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



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

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]

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

The code of the used auxiliary function "IsPrime(n)" is located at Primality by trial division#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



Module Prime_decomposition {

Inventory [email protected], [email protected]

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 = [email protected]

{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

{

[email protected]

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



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:

> ifactor(1337);

(7) (191)



> ifactors(1337);

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



Bare built-in function does:

FactorInteger[2016] => {{2, 5}, {3, 2}, {7, 1}}

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:

Example for small prime:

ShowPrimeDecomposition[1337]

gives:

1337 = 7 191

Examples for large primes:

Table[AbsoluteTiming[ShowPrimeDecomposition[2^a-1]]//Print[#[[1]]," sec"]&,{a,50,150,10}];

gives back:

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

function [ outputPrimeDecomposition ] = primedecomposition ( inputValue )

outputPrimeDecomposition = factor ( inputValue ) ;

Using the built-in function:

(%i1) display2d: false$ /* disable rendering exponents as superscripts */

(%i2) factor(2016);

(%o2) 2^5*3^2*7



Using the underlying language:

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

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

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

Based on python solution:

import strutils, math, sequtils, times



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

result = 1 + n*4 - int64(n /% 2)*2



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

var res: seq[int64] = @[]

var maxq = int64(floor(sqrt(float(n))))

var d = 1

var q: int64 = (n %% 2) and 2 or 3 # either 2 or 3, alternating

while (q <= maxq) and ((n %% q) != 0):

q = getStep(d)

d += 1

if q <= maxq:

var q1: seq[int64] = primeFac(n /% q)

var q2: seq[int64] = primeFac(q)

res = concat(q2, q1, res)

else:

res.add(n)

result = res



var is_prime: seq[Bool] = @[]

is_prime.add(False)

is_prime.add(False)



iterator primes(limit: int): int =

for n in high(is_prime) .. limit+2: is_prime.add(True)

for n in 2 .. limit + 1:

if is_prime[n]:

yield n

for i in countup((n *% n), limit+1, n): # start at ``n`` squared

try:

is_prime[i] = False

except EInvalidIndex: break



# Example: calculate factors of Mersenne numbers to M59 #



for m in primes(59):

var p = int64(pow(2.0,float(m)) - 1)

write(stdout,"2**$1-1 = $2, with factors: " % [$m, $p] )

var start = cpuTime()

var f = primeFac(p)

for factor in f:

write(stdout, factor)

write(stdout, ", ")

FlushFile(stdout)

writeln(stdout, "=> $#ms" % $int(1000*(cpuTime()-start)) )

Output:

compiled with options -x:off -opt:speed

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

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

r = factor ( 120202039393 )

Oforth handles aribitrary precision integers.

: 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 ;

Output:

>2 128 pow 1 - dup println factors println 340282366920938463463374607431768211455 [3, 5, 17, 257, 641, 65537, 274177, 6700417, 67280421310721] ok

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:

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 )

} ;

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 .

Output:

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

Optimization:

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 .

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

Trivial trial division (very slow) [ edit ]

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

" ;

Better trial division [ edit ]

This is much faster than the trivial version above.

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 ;

}

Modules [ edit ]

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.

use ntheory qw/factor forprimes/ ;

use bigint ;



forprimes {

my $p = 2 ** $_ - 1 ;

print "2**$_-1: " , join ( " " , factor ( $p ) ) , "

" ;

} 100 , 150 ;

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

use Math :: Pari qw/:int factorint isprime/ ;



# 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 ) ) , "

" ;

}

With the same output.

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

include mpfr.e

atom t0 = time()

mpz z = mpz_init()

for i=1 to 17 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_prime_factors(z,20000))

if fs!=zs then zs &= " = "&fs end if

printf(1,"2^%d-1 = %s

",{pi,zs})

end for

string s = "600851475143"

mpz_set_str(z,s)

printf(1,"%s = %s

",{s,mpz_factorstring(mpz_prime_factors(z,500))})

?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 600851475143 = 71*839*1471*6857 "0.1s"

Note that mpz_prime_factors() needs to be told how far to push things before giving up, but if pushed to (say) 20,000,000 primes, performance can suffer quite dramatically.

t0 = time()

--for i=1 to 17 do

for i=18 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_prime_factors(z,20000000))

if fs!=zs then zs &= " = "&fs end if

printf(1,"2^%d-1 = %s

",{pi,zs})

end for

s = "100000000000000000037"

mpz_set_str(z,s)

printf(1,"%s = %s

",{s,mpz_factorstring(mpz_prime_factors(z,5000000))})

?elapsed(time()-t0)

Output:

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 100000000000000000037 = 31*821*59004541*66590107 "23.1s"

Whereas the code below produces this:

2^61-1 = 2305843009213693951 2^67-1 = 147573952589676412927 (not prime) 2^71-1 = 2361183241434822606847 (not prime) 2^73-1 = 9444732965739290427391 = 439*21514198099633918969 (last factor is not prime) 2^79-1 = 604462909807314587353087 (not prime) 2^83-1 = 9671406556917033397649407 = 167*57912614113275649087721 2^89-1 = 618970019642690137449562111 2^97-1 = 158456325028528675187087900671 (not prime) "0.0s"

The default of 100 (as in get_prime(100) yields 541) is quite low, but fast (as is that 20,000 above).

Obviously, were you not actually going to make any use of factors>541, then that's all you'd need.

Below I commented said changes to mpz_prime_factors() and added the "(as per docs)" part.

Note the latter was not needed above, via careful choice of 20000/500/20000000/5000000.

include mpfr.e

atom t0 = time()

mpz z = mpz_init()

randstate state = gmp_randinit_mt()

for i=18 to 25 do

integer pi = get_prime(i)

mpz_ui_pow_ui(z,2,pi)

mpz_sub_ui(z,z,1)

-- sequence f = mpz_prime_factors(z,20000000)

sequence f = mpz_prime_factors(z)

string zs = mpz_get_str(z),

fs = mpz_factorstring(f)

if fs!=zs then zs &= " = "&fs end if

if length(f[$])=1 then -- (as per docs)

mpz_set_str(z,f[$][1])

if not mpz_probable_prime_p(z, state) then

if length(f)=1 then

zs &= " (not prime)"

else

zs &= " (last factor is not prime)"

end if

end if

end if

printf(1,"2^%d-1 = %s

",{pi,zs})

end for

?elapsed(time()-t0)

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.

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

Output:

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



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 ;



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 ]



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



Output:

2 2 3 2 2 5 5

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 )

)

) .

Output:

?- 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 ] .

Simple version [ edit ]

Translation of: Erlang

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

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 )

) .

Expression Tree version [ edit ]

Uses a 2*3*5*7 factor wheel, but the main feature is that it returns the decomposition as a fully simplified expression tree.



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



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.

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;

Works with: PureBasic version 4.41



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 )

Output:

Factored 9007199254740991 in 0.27 seconds. 6361 69431 20394401

Python: Using Croft Spiral sieve [ edit ]

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

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

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 [ edit ]

Here a shorter and marginally faster algorithm:

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

Output:

576460752303423487 = [3203431780337, 179951] Needed 0.9240529537200928s

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

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

Recursive Approach (Less efficient for large numbers) [ edit ]



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

}



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

}

}



Alternate solution [ edit ]

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

(require math)

(define (factors n)

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



Or, an explicit (and less efficient) computation:



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



(formerly Perl 6)

Pure Raku [ edit ]

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⁴⁰.

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

}

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 [ edit ]

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.

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

}

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.

optimized slightly [ edit ]

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 (t