Prime decomposition
You are encouraged to solve this task according to the task description, using any language you may know.
The prime decomposition of a number is defined as a list of prime numbers which when all multiplied together, are equal to that number. Example: 12 = 2 × 2 × 3, so its prime decomposition is {2, 2, 3}
Write a function which returns an array or collection which contains the prime decomposition of a given number, n, greater than 1. If your language does not have an isPrime-like function available, you may assume that you have a function which determines whether a number is prime (note its name before your code).
If you would like to test code from this task, you may use code from trial division or the Sieve of Eratosthenes.
Note: The program must not be limited by the word size of your computer or some other artificial limit; it should work for any number regardless of size (ignoring the physical limits of RAM etc).
ABAP
<lang ABAP>class ZMLA_ROSETTA definition
public create public .
public section.
types:
enumber TYPE N LENGTH 60,
listof_enumber TYPE TABLE OF enumber .
class-methods FACTORS
importing
value(N) type ENUMBER
exporting
value(ORET) type LISTOF_ENUMBER .
protected section.
private section.
ENDCLASS.
CLASS ZMLA_ROSETTA IMPLEMENTATION.
- <SIGNATURE>---------------------------------------------------------------------------------------+
- | Static Public Method ZMLA_ROSETTA=>FACTORS
- +-------------------------------------------------------------------------------------------------+
- | [--->] N TYPE ENUMBER
- | [<---] ORET TYPE LISTOF_ENUMBER
- +--------------------------------------------------------------------------------------</SIGNATURE>
method FACTORS.
CLEAR oret.
WHILE n mod 2 = 0.
n = n / 2.
APPEND 2 to oret.
ENDWHILE.
DATA: lim type enumber,
i type enumber.
lim = sqrt( n ).
i = 3.
WHILE i <= lim.
WHILE n mod i = 0.
APPEND i to oret.
n = n / i.
lim = sqrt( n ).
ENDWHILE.
i = i + 2.
ENDWHILE.
IF n > 1.
APPEND n to oret.
ENDIF.
endmethod.
ENDCLASS.</lang>
ACL2
<lang Lisp>(include-book "arithmetic-3/top" :dir :system)
(defun prime-factors-r (n i)
(declare (xargs :mode :program))
(cond ((or (zp n) (zp (- n i)) (zp i) (< i 2) (< n 2))
(list n))
((= (mod n i) 0)
(cons i (prime-factors-r (floor n i) 2)))
(t (prime-factors-r n (1+ i)))))
(defun prime-factors (n)
(declare (xargs :mode :program)) (prime-factors-r n 2))</lang>
Ada
<lang ada>with Ada.Text_IO; use Ada.Text_IO;
procedure Test_Prime is
generic
type Number is private;
Zero : Number;
One : Number;
Two : Number;
with function Image (X : Number) return String is <>;
with function "+" (X, Y : Number) return Number is <>;
with function "/" (X, Y : Number) return Number is <>;
with function "mod" (X, Y : Number) return Number is <>;
with function ">=" (X, Y : Number) return Boolean is <>;
package Prime_Numbers is
type Number_List is array (Positive range <>) of Number;
function Decompose (N : Number) return Number_List;
procedure Put (List : Number_List);
end Prime_Numbers;
package body Prime_Numbers is
function Decompose (N : Number) return Number_List is
Size : Natural := 0;
M : Number := N;
K : Number := Two;
begin
-- Estimation of the result length from above
while M >= Two loop
M := (M + One) / Two;
Size := Size + 1;
end loop;
M := N;
-- Filling the result with prime numbers
declare
Result : Number_List (1..Size);
Index : Positive := 1;
begin
while N >= K loop -- Divisors loop
while Zero = (M mod K) loop -- While divides
Result (Index) := K;
Index := Index + 1;
M := M / K;
end loop;
K := K + One;
end loop;
return Result (1..Index - 1);
end;
end Decompose;
procedure Put (List : Number_List) is
begin
for Index in List'Range loop
Put (Image (List (Index)));
end loop;
end Put;
end Prime_Numbers;
package Integer_Numbers is new Prime_Numbers (Natural, 0, 1, 2, Positive'Image);
use Integer_Numbers;
begin
Put (Decompose (12));
end Test_Prime;</lang> The solution is generic. The package is instantiated by a type that supports necessary operations +, /, mod, >=. The constants 0, 1, 2 are parameters too, because the type might have no literals. The package also provides a procedure to output an array of prime numbers and a function to convert a number to string (as a parameter). The function Decompose first estimates the maximal result length as log2 of the argument. Then it allocates the result and starts to enumerate divisors. It does not care to check if the divisors are prime, because non-prime divisors will be automatically excluded. In the example provided, the package is instantiated with plain integer type. Sample output:
2 2 3
ALGOL 68
- note: This specimen retains the original Python coding style.
<lang algol68>#IF long int possible THEN #
MODE LINT = LONG INT; LINT lmax int = long max int; OP LLENG = (INT i)LINT: LENG i,
LSHORTEN = (LINT i)INT: SHORTEN i;
- 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
)</lang> Output:
2**2-1 = 3, with factors: 3 2**3-1 = 7, with factors: 7 2**5-1 = 31, with factors: 31 2**7-1 = 127, with factors: 127 2**11-1 = 2047, with factors: 23 89 2**13-1 = 8191, with factors: 8191 2**17-1 = 131071, with factors: 131071 2**19-1 = 524287, with factors: 524287 2**23-1 = 8388607, with factors: 47 178481 2**29-1 = 536870911, with factors: 233 1103 2089 2**31-1 = 2147483647, with factors: 2147483647 2**37-1 = 137438953471, with factors: 223 616318177 2**41-1 = 2199023255551, with factors: 13367 164511353 2**43-1 = 8796093022207, with factors: 431 9719 2099863 2**47-1 = 140737488355327, with factors: 2351 4513 13264529 2**53-1 = 9007199254740991, with factors: 6361 69431 20394401 2**59-1 = 576460752303423487, with factors: 179951 3203431780337
Note: ALGOL 68G took 49,109,599 BogoMI and ELLA ALGOL 68RS took 1,127,634 BogoMI to complete the example.
AutoHotkey
<lang AutoHotkey>MsgBox % factor(8388607) ; 47 * 178481
factor(n) {
If (n = 1)
Return
f = 2
While (f <= n)
{
If (Mod(n, f) = 0)
{
next := factor(n / f)
factors = %f%`n%next%
Return factors
}
f++
}
}</lang>
AWK
As the examples show, pretty large numbers can be factored in tolerable time:
<lang awk>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) }</lang>
$ awk -f primefac.awk 36 2 2 3 3 77 7 11 536870911 233 1103 2089 8796093022207 431 9719 2099863
Commodore BASIC
It's not easily possible to have arbitrary precision integers in PET basic, so here is at least a version using built-in data types (reals). On return from the subroutine starting at 9000 the global array pf contains the number of factors followed by the factors themselves: <lang zxbasic>9000 REM ----- function generate 9010 REM in ... i ... number 9020 REM out ... pf() ... factors 9030 REM mod ... ca ... pf candidate 9040 pf(0)=0 : ca=2 : REM special case 9050 IF i=1 THEN RETURN 9060 IF INT(i/ca)*ca=i THEN GOSUB 9200 : GOTO 9050 9070 FOR ca=3 TO INT( SQR(i)) STEP 2 9080 IF i=1 THEN RETURN 9090 IF INT(i/ca)*ca=i THEN GOSUB 9200 : GOTO 9080 9100 NEXT 9110 IF i>1 THEN ca=i : GOSUB 9200 9120 RETURN 9200 pf(0)=pf(0)+1 9210 pf(pf(0))=ca 9220 i=i/ca 9230 RETURN</lang>
Befunge
Handles safely integers only up to 250 (or ones which don't have prime divisors greater than 250). <lang Befunge>& 211p > : 1 - #v_ 25*, @ > 11g:. / v
> : 11g %!|
> 11g 1+ 11p v
^ <</lang>
Burlesque
<lang burlesque> blsq ) 12fC {2 2 3} </lang>
C
Relatively sophiscated sieve method based on size 30 prime wheel. The code does not pretend to handle prime factors larger than 64 bits. All 32-bit primes are cached with 137MB data. Cache data takes about a minute to compute the first time the program is run, which is also saved to the current directory, and will be loaded in a second if needed again. <lang c>#include <inttypes.h>
- 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"\n", s); } sieve(s); } fp = fopen("primebits", "w"); fwrite(pbits, 1, PBITS, fp); fclose(fp); }
int is_prime(xint x) { pint p; if (x > 5) { if (x < MAX_PRIME) return pbits[x/30] & bit_pos[x % 30];
for (p = 2; p && (xint)p * p <= x; p = next_prime(p)) if (x % p == 0) return 0;
return 1; } return x == 2 || x == 3 || x == 5; }
void sieve(pint p) { unsigned char b[8]; off_t ofs[8]; int i, q;
for (i = 0; i < 8; i++) { q = rem_num[i] * p; b[i] = ~bit_pos[q % 30]; ofs[i] = q / 30; }
for (q = ofs[1], i = 7; i; i--) ofs[i] -= ofs[i-1];
for (ofs[0] = p, i = 1; i < 8; i++) ofs[0] -= ofs[i];
for (i = 1; q < PBITS; q += ofs[i = (i + 1) & 7]) pbits[q] &= b[i]; }
pint next_prime(pint p) { off_t addr; uint8_t bits, rem;
if (p > 5) { addr = p / 30; bits = bit_pos[ p % 30 ] << 1; for (rem = 0; (1 << rem) < bits; rem++); while (pbits[addr] < bits || !bits) { if (++addr >= PBITS) return 0; bits = 1; rem = 0; } if (addr >= PBITS) return 0; while (!(pbits[addr] & bits)) { rem++; bits <<= 1; } return p = addr * 30 + rem_num[rem]; }
switch(p) { case 2: return 3; case 3: return 5; case 5: return 7; } return 2; }
int decompose(xint n, xint *f) { pint p = 0; int i = 0;
/* check small primes: not strictly necessary */ if (n <= MAX_PRIME && is_prime(n)) { f[0] = n; return 1; }
while (n >= (xint)p * p) { if (!(p = next_prime(p))) break; while (n % p == 0) { n /= p; f[i++] = p; } } if (n > 1) f[i++] = n; return i; }
int main() { int i, len; pint p = 0; xint f[MAX_FACTORS], po;
init_primes();
for (p = 1; p < 64; p++) { po = (1LLU << p) - 1; printf("2^%"PRIuPINT" - 1 = %"PRIuXINT, p, po); fflush(stdout); if ((len = decompose(po, f)) > 1) for (i = 0; i < len; i++) printf(" %c %"PRIuXINT, i?'x':'=', f[i]); putchar('\n'); }
return 0; }</lang>
Using GNU Compiler Collection gcc extensions
Note: The following code sample is experimental as it implements python style iterators for (potentially) infinite sequences. C is not normally written this way, and in the case of this sample it requires the GCC "nested procedure" extension to the C language. <lang c>#include <limits.h>
- 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("\n",m);
if( m >= 59 )BREAK;
CONTINUE;
}
}</lang> Output:
2**2-1 = 3, with factors: 3 2**3-1 = 7, with factors: 7 2**5-1 = 31, with factors: 31 2**7-1 = 127, with factors: 127 2**11-1 = 2047, with factors: 23 89 2**13-1 = 8191, with factors: 8191 2**17-1 = 131071, with factors: 131071 2**19-1 = 524287, with factors: 524287 2**23-1 = 8388607, with factors: 47 178481 2**29-1 = 536870911, with factors: 233 1103 2089 2**31-1 = 2147483647, with factors: 2147483647 2**37-1 = 137438953471, with factors: 223 616318177 2**41-1 = 2199023255551, with factors: 13367 164511353 2**43-1 = 8796093022207, with factors: 431 9719 2099863 2**47-1 = 140737488355327, with factors: 2351 4513 13264529 2**53-1 = 9007199254740991, with factors: 6361 69431 20394401 2**59-1 = 576460752303423487, with factors: 179951 3203431780337
Note: gcc took 487,719 BogoMI to complete the example.
To understand what was going on with the above code, pass it through cpp and read the outcome. Translated into normal C code sans the function call overhead, it's really this (the following uses a adjustable cache, although setting it beyond a few thousands doesn't gain further benefit):<lang c>#include <stdio.h>
- 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 " : "\n"); }
return 0; }</lang>
C++
<lang cpp>#include <iostream>
- include <gmpxx.h>
// This function template works for any type representing integers or // nonnegative integers, and has the standard operator overloads for // arithmetic and comparison operators, as well as explicit conversion // from int. // // OutputIterator must be an output iterator with value_type Integer. // It receives the prime factors. template<typename Integer, typename OutputIterator>
void decompose(Integer n, OutputIterator out)
{
Integer i(2);
while (n != 1)
{
while (n % i == Integer(0))
{
*out++ = i;
n /= i;
}
++i;
}
}
// this is an output iterator similar to std::ostream_iterator, except // that it outputs the separation string *before* the value, but not // before the first value (i.e. it produces an infix notation). template<typename T> class infix_ostream_iterator:
public std::iterator<T, std::output_iterator_tag>
{
class Proxy;
friend class Proxy;
class Proxy
{
public:
Proxy(infix_ostream_iterator& iter): iterator(iter) {}
Proxy& operator=(T const& value)
{
if (!iterator.first)
{
iterator.stream << iterator.infix;
}
iterator.stream << value;
}
private:
infix_ostream_iterator& iterator;
};
public:
infix_ostream_iterator(std::ostream& os, char const* inf):
stream(os),
first(true),
infix(inf)
{
}
infix_ostream_iterator& operator++() { first = false; return *this; }
infix_ostream_iterator operator++(int)
{
infix_ostream_iterator prev(*this);
++*this;
return prev;
}
Proxy operator*() { return Proxy(*this); }
private:
std::ostream& stream; bool first; char const* infix;
};
int main() {
std::cout << "please enter a positive number: ";
mpz_class number;
std::cin >> number;
if (number <= 0)
std::cout << "this number is not positive!\n;";
else
{
std::cout << "decomposition: ";
decompose(number, infix_ostream_iterator<mpz_class>(std::cout, " * "));
std::cout << "\n";
}
}</lang>
C#
<lang csharp>using System; using System.Collections.Generic;
namespace PrimeDecomposition {
class Program
{
static void Main(string[] args)
{
getPrimes(12);
}
static List<int> getPrimes(decimal n)
{
List<int> storage = new List<int>();
while (n > 1)
{
int i = 1;
while (true)
{
if (isPrime(i))
{
if (((decimal)n / i) == Math.Round((decimal) n / i))
{
n /= i;
storage.Add(i);
break;
}
}
i++;
}
}
return storage;
}
static bool isPrime(int n)
{
if (n <= 1) return false;
for (int i = 2; i <= Math.Sqrt(n); i++)
if (n % i == 0) return false;
return true;
}
}
}</lang>
Simple trial division
This version a translation from Java of the sample presented by Robert C. Martin during a TDD talk at NDC 2011.
Although this three-line algorithm does not mention anything about primes, the fact that factors are taken out of the number n in ascending order garantees the list will only contain primes.
<lang csharp>using System.Collections.Generic;
namespace PrimeDecomposition { public class Primes {
public List<int> FactorsOf(int n)
{ var factors = new List<int>();
for (var divisor = 2; n > 1; divisor++) for (; n % divisor == 0; n /= divisor) factors.Add(divisor);
return factors; } }</lang>
Clojure
<lang clojure>;;; No stack consuming algorithm (defn factors
"Return a list of factors of N."
([n]
(factors n 2 ()))
([n k acc]
(if (= 1 n)
acc
(if (= 0 (rem n k))
(recur (quot n k) k (cons k acc))
(recur n (inc k) acc)))))</lang>
Common Lisp
<lang Lisp>;;; Recursive algorithm (defun factor (n)
"Return a list of factors of N." (when (> n 1) (loop with max-d = (isqrt n)
for d = 2 then (if (evenp d) (+ d 1) (+ d 2)) do (cond ((> d max-d) (return (list n))) ; n is prime ((zerop (rem n d)) (return (cons d (factor (truncate n d)))))))))</lang>
D
<lang d>import std.traits: Unqual;
Unqual!T[] decompose(T)(T number) /*pure nothrow*/ in {
assert(number > 1);
} body {
alias UT = Unqual!T; typeof(return) result; UT n = number;
for (UT i = 2; n % i == 0;) {
result ~= i;
n /= i;
}
for (UT i = 3; n >= i * i; i += 2) {
while (n % i == 0) {
result ~= i;
n /= i;
}
}
if (n != 1)
result ~= n;
return result;
}
void main() {
import std.stdio, std.bigint, std.algorithm;
foreach (immutable n; 2 .. 10)
writeln(decompose(n));
writeln(decompose(1023 * 1024)); writeln(decompose(BigInt(2 * 3 * 5 * 7 * 11 * 11 * 13 * 17))); writeln(decompose(BigInt(16860167264933UL) * 179951)); writeln(group(decompose(BigInt(2) ^^ 100_000)));
}</lang>
- Output:
[2] [3] [2, 2] [5] [2, 3] [7] [2, 2, 2] [3, 3] [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 11, 31] [2, 3, 5, 7, 11, 11, 13, 17] [179951, 16860167264933] [Tuple!(BigInt, uint)(2, 100000)]
E
This example assumes a function isPrime and was tested with this one. It could use a self-referential implementation such as the Python task, but the original author of this example did not like the ordering dependency involved.
<lang e>def primes := {
var primesCache := [2]
/** A collection of all prime numbers. */
def primes {
to iterate(f) {
primesCache.iterate(f)
for x in (int > primesCache.last()) {
if (isPrime(x)) {
f(primesCache.size(), x)
primesCache with= x
}
}
}
}
}
def primeDecomposition(var x :(int > 0)) {
var factors := []
for p in primes {
while (x % p <=> 0) {
factors with= p
x //= p
}
if (x <=> 1) {
break
}
}
return factors
}</lang>
Erlang
<lang erlang>% no stack consuming version
factors(N) ->
factors(N,2,[]).
factors(1,_,Acc) -> Acc; factors(N,K,Acc) when N rem K == 0 ->
factors(N div K,K, [K|Acc]);
factors(N,K,Acc) ->
factors(N,K+1,Acc).</lang>
F#
<lang Fsharp>let decompose_prime n =
let rec loop c p = if c < (p * p) then [c] elif c % p = 0I then p :: (loop (c/p) p) else loop c (p + 1I) loop n 2I
decompose_prime 600851475143I</lang>
FALSE
<lang false>[2[\$@$$*@>~][\$@$@$@$@\/*=$[%$." "$@\/\0~]?~[1+1|]?]#%.]d: 27720d;! {2 2 2 3 3 5 7 11}</lang>
Factor
Word factors from dictionary math.primes.factors converts a number into a sequence of its prime divisors; the rest of the code prints this sequence. <lang factor>USING: io kernel math math.parser math.primes.factors sequences ;
27720 factors [ number>string ] map " " join print ;</lang>
Forth
<lang forth>: decomp ( n -- )
2
begin 2dup dup * >=
while 2dup /mod swap
if drop 1+ 1 or \ next odd number
else -rot nip dup .
then
repeat
drop . ;</lang>
Fortran
<lang fortran>module PrimeDecompose
implicit none
integer, parameter :: huge = selected_int_kind(18) ! => integer(8) ... more fails on my 32 bit machine with gfortran(gcc) 4.3.2
contains
subroutine find_factors(n, d) integer(huge), intent(in) :: n integer, dimension(:), intent(out) :: d
integer(huge) :: div, next, rest integer :: i
i = 1
div = 2; next = 3; rest = n
do while ( rest /= 1 )
do while ( mod(rest, div) == 0 )
d(i) = div
i = i + 1
rest = rest / div
end do
div = next
next = next + 2
end do
end subroutine find_factors
end module PrimeDecompose</lang>
<lang fortran>program Primes
use PrimeDecompose implicit none
integer, dimension(100) :: outprimes integer i
outprimes = 0
call find_factors(12345649494449_huge, outprimes)
do i = 1, 100
if ( outprimes(i) == 0 ) exit
print *, outprimes(i)
end do
end program Primes</lang>
Frink
Frink has a built-in factoring function which uses wheel factoring, trial division, Pollard p-1 factoring, and Pollard rho factoring. It also recognizes some special forms (e.g. Mersenne numbers) and handles them efficiently. <lang frink>println[factor[2^508-1]]</lang>
Output (total process time including JVM startup = 1.515 s):
[[3, 1], [5, 1], [509, 1], [18797, 1], [26417, 1], [72118729, 1], [140385293, 1], [2792688414613, 1], [8988357880501, 1], [90133566917913517709497, 1], [56713727820156410577229101238628035243, 1], [170141183460469231731687303715884105727, 1]]
Note that this means 31 * 51 * ...
GAP
Built-in function : <lang gap>FactorsInt(2^67-1);
- [ 193707721, 761838257287 ]</lang>
Or using the FactInt package : <lang gap>FactInt(2^67-1);
- [ [ 193707721, 761838257287 ], [ ] ]</lang>
Go
<lang go>package main
import ( "fmt" "math/big" )
var ( ZERO = big.NewInt(0) ONE = big.NewInt(1) )
func Primes(n *big.Int) []*big.Int { res := []*big.Int{} mod, div := new(big.Int), new(big.Int) for i := big.NewInt(2); i.Cmp(n) != 1; { div.DivMod(n, i, mod) for mod.Cmp(ZERO) == 0 { res = append(res, new(big.Int).Set(i)) n.Set(div) div.DivMod(n, i, mod) } i.Add(i, ONE) } return res }
func main() { vals := []int64{ 1 << 31, 1234567, 333333, 987653, 2 * 3 * 5 * 7 * 11 * 13 * 17, } for _, v := range vals { fmt.Println(v, "->", Primes(big.NewInt(v))) } }</lang> Output:
2147483648 -> [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2] 1234567 -> [127 9721] 333333 -> [3 3 7 11 13 37] 987653 -> [29 34057] 510510 -> [2 3 5 7 11 13 17]
Groovy
This solution uses the fact that a given factor must be prime if no smaller factor divides it evenly, so it does not require an "isPrime-like function", assumed or otherwise. <lang groovy>def factorize = { long target ->
if (target == 1) return [1L]
if (target < 4) return [1L, target]
def targetSqrt = Math.sqrt(target)
def lowfactors = (2L..targetSqrt).findAll { (target % it) == 0 }
if (lowfactors == []) return [1L, target]
def nhalf = lowfactors.size() - ((lowfactors[-1]**2 == target) ? 1 : 0)
[1] + lowfactors + (0..<nhalf).collect { target.intdiv(lowfactors[it]) }.reverse() + [target]
}
def decomposePrimes = { target ->
def factors = factorize(target) - [1]
def primeFactors = []
factors.eachWithIndex { f, i ->
if (i==0 || factors[0..<i].every {f % it != 0}) {
primeFactors << f
def pfPower = f*f
while (target % pfPower == 0) {
primeFactors << f
pfPower *= f
}
}
}
primeFactors
}</lang>
Test #1: <lang groovy>((1..30) + [97*4, 1000, 1024, 333333]).each { println ([number:it, primes:decomposePrimes(it)]) }</lang>
Output #1:
[number:1, primes:[]] [number:2, primes:[2]] [number:3, primes:[3]] [number:4, primes:[2, 2]] [number:5, primes:[5]] [number:6, primes:[2, 3]] [number:7, primes:[7]] [number:8, primes:[2, 2, 2]] [number:9, primes:[3, 3]] [number:10, primes:[2, 5]] [number:11, primes:[11]] [number:12, primes:[2, 2, 3]] [number:13, primes:[13]] [number:14, primes:[2, 7]] [number:15, primes:[3, 5]] [number:16, primes:[2, 2, 2, 2]] [number:17, primes:[17]] [number:18, primes:[2, 3, 3]] [number:19, primes:[19]] [number:20, primes:[2, 2, 5]] [number:21, primes:[3, 7]] [number:22, primes:[2, 11]] [number:23, primes:[23]] [number:24, primes:[2, 2, 2, 3]] [number:25, primes:[5, 5]] [number:26, primes:[2, 13]] [number:27, primes:[3, 3, 3]] [number:28, primes:[2, 2, 7]] [number:29, primes:[29]] [number:30, primes:[2, 3, 5]] [number:388, primes:[2, 2, 97]] [number:1000, primes:[2, 2, 2, 5, 5, 5]] [number:1024, primes:[2, 2, 2, 2, 2, 2, 2, 2, 2, 2]] [number:333333, primes:[3, 3, 7, 11, 13, 37]]
Test #2: <lang groovy>def isPrime = {factorize(it).size() == 2} (1..60).step(2).findAll(isPrime).each { println ([number:"2**${it}-1", value:2**it-1, primes:decomposePrimes(2**it-1)]) }</lang>
Output #2:
[number:2**3-1, value:7, primes:[7]] [number:2**5-1, value:31, primes:[31]] [number:2**7-1, value:127, primes:[127]] [number:2**11-1, value:2047, primes:[23, 89]] [number:2**13-1, value:8191, primes:[8191]] [number:2**17-1, value:131071, primes:[131071]] [number:2**19-1, value:524287, primes:[524287]] [number:2**23-1, value:8388607, primes:[47, 178481]] [number:2**29-1, value:536870911, primes:[233, 1103, 2089]] [number:2**31-1, value:2147483647, primes:[2147483647]] [number:2**37-1, value:137438953471, primes:[223, 616318177]] [number:2**41-1, value:2199023255551, primes:[13367, 164511353]] [number:2**43-1, value:8796093022207, primes:[431, 9719, 2099863]] [number:2**47-1, value:140737488355327, primes:[2351, 4513, 13264529]] [number:2**53-1, value:9007199254740991, primes:[6361, 69431, 20394401]] [number:2**59-1, value:576460752303423487, primes:[179951, 3203431780337]]
Perhaps a more sophisticated algorithm is in order. It took well over 1 hour to calculate the last three decompositions using this solution.
Haskell
The task description hints at using isPrime function:
<lang haskell>factorize_ n | n > 1 = concat [divs n p | p <- [2..n], isPrime p]
where divs n p = if rem n p==0 then p:divs (quot n p) p else []</lang>
but it is not very efficient, if at all. Inlining and optimizing gets:
<lang haskell>factorize n | n > 1 = go n primesList
where
go n ds@(d:t)
| d*d > n = [n]
| r == 0 = d : go q ds
| otherwise = go n t
where
(q,r) = quotRem n d</lang>
See Sieve of Eratosthenes or Primality by trial division for a source of primes to use with this function. Actually as some other entries notice, for any ascending order list containing all primes, used in place of primesList, the factors found by this function are guaranteed to be prime, so no separate testing for primality is needed; however using just primes is more efficient.
Yet another way using filters and test driven development:
<lang haskell>import Test.HUnit import Data.List
factor::Int->[Int] factor 1 = [] factor n
| product primeDivisor == n = primeDivisor | otherwise = primeDivisor ++ factor (div n $ product primeDivisor) where primeDivisor = filter isPrime $ filter (isDivisor n) [2 .. n] isDivisor n d = (==) 0 $ mod n d isPrime n = not $ any (isDivisor n) [2 .. n-1]
tests = TestList[TestCase $ assertEqual "1 has no prime factors" [] $ factor 1
,TestCase $ assertEqual "2 is 2" [2] $ factor 2
,TestCase $ assertEqual "3 is 3" [3] $ factor 3
,TestCase $ assertEqual "4 is 2*2" [2, 2] $ factor 4
,TestCase $ assertEqual "5 is 5" [5] $ factor 5
,TestCase $ assertEqual "6 is 2*3" [2, 3] $ factor 6
,TestCase $ assertEqual "7 is 7" [7] $ factor 7
,TestCase $ assertEqual "8 is 2*2*2" [2, 2, 2] $ factor 8
,TestCase $ assertEqual "large number" [2, 3, 3, 5, 7, 11, 13] $ sort $ factor (2*9*5*7*11*13)]</lang>
Icon and Unicon
<lang Icon>procedure main() factors := primedecomp(2^43-1) # a big int end
procedure primedecomp(n) #: return a list of factors local F,o,x F := []
every writes(o,n|(x := genfactors(n))) do {
\o := "*" /o := "=" put(F,x) # build a list of factors to satisfy the task }
write() return F end
link factors</lang>
Uses genfactors and prime from factors
Sample Output showing factors of a large integer:
8796093022207=431*9719*2099863
J
<lang j>q:</lang>
Example use: <lang j>q: 3684 2 2 3 307
_1+2^128x
340282366920938463463374607431768211455
q: _1+2^128x
3 5 17 257 641 65537 274177 6700417 67280421310721
*/ q: _1+2^128x
340282366920938463463374607431768211455</lang>
Java
This is a version for arbitrary-precision integers which assumes the existence of a function with the signature: <lang java>public boolean prime(BigInteger i);</lang> You will need to import java.util.List, java.util.LinkedList, and java.math.BigInteger. <lang java>public static List<BigInteger> primeFactorBig(BigInteger a){
List<BigInteger> ans = new LinkedList<BigInteger>();
//loop until we test the number itself or the number is 1
for (BigInteger i = BigInteger.valueOf(2); i.compareTo(a) <= 0 && !a.equals(BigInteger.ONE);
i = i.add(BigInteger.ONE)){
while (a.remainder(i).equals(BigInteger.ZERO) && prime(i)) { //if we have a prime factor
ans.add(i); //put it in the list
a = a.divide(i); //factor it out of the number
}
}
return ans;
}</lang>
Alternate version, optimised to be faster. <lang java>private static final BigInteger two = BigInteger.valueOf(2);
public List<BigInteger> primeDecomp(BigInteger a) {
// impossible for values lower than 2
if (a.compareTo(two) < 0) {
return null;
}
//quickly handle even values
List<BigInteger> result = new ArrayList<BigInteger>();
while (a.and(BigInteger.ONE).equals(BigInteger.ZERO)) {
a = a.shiftRight(1);
result.add(two);
}
//left with odd values
if (!a.equals(BigInteger.ONE)) {
BigInteger b = BigInteger.valueOf(3);
while (b.compareTo(a) < 0) {
if (b.isProbablePrime(10)) {
BigInteger[] dr = a.divideAndRemainder(b);
if (dr[1].equals(BigInteger.ZERO)) {
result.add(b);
a = dr[0];
}
}
b = b.add(two);
}
result.add(b); //b will always be prime here...
}
return result;
}</lang>
Simple but very inefficient method, because it will test divisibility of all numbers from 2 to max prime factor. When decomposing a large prime number this will take O(n) trial divisions instead of more common O(log n). <lang java>public static List<BigInteger> primeFactorBig(BigInteger a){
List<BigInteger> ans = new LinkedList<BigInteger>();
for(BigInteger divisor = BigInteger.valueOf(2); a.compareTo(ONE) > 0; divisor = divisor.add(ONE))
while(a.mod(divisor).equals(ZERO)){ ans.add(divisor); a = a.divide(divisor); }
return ans;
}</lang>
JavaScript
This code uses the BigInteger Library jsbn and jsbn2 <lang javascript>function run_factorize(input, output) {
var n = new BigInteger(input.value, 10);
var TWO = new BigInteger("2", 10);
var divisor = new BigInteger("3", 10);
var prod = false;
if (n.compareTo(TWO) < 0)
return;
output.value = "";
while (true) {
var qr = n.divideAndRemainder(TWO);
if (qr[1].equals(BigInteger.ZERO)) {
if (prod)
output.value += "*";
else
prod = true;
output.value += "2";
n = qr[0];
}
else
break;
}
while (!n.equals(BigInteger.ONE)) {
var qr = n.divideAndRemainder(divisor);
if (qr[1].equals(BigInteger.ZERO)) {
if (prod)
output.value += "*";
else
prod = true;
output.value += divisor;
n = qr[0];
}
else
divisor = divisor.add(TWO);
}
}</lang>
Without any library. <lang javascript>function run_factorize(n) { if (n <= 3) return [n];
var ans = []; var done = false; while (!done) { if (n%2 === 0){ ans.push(2); n /= 2; continue; } if (n%3 === 0){ ans.push(3); n /= 3; continue; } if ( n === 1) return ans; var sr = Math.sqrt(n); done = true; // try to divide the checked number by all numbers till its square root. for (var i=6; i<=sr; i+=6){ if (n%(i-1) === 0){ // is n divisible by i-1? ans.push( (i-1) ); n /= (i-1); done = false; break; } if (n%(i+1) === 0){ // is n divisible by i+1? ans.push( (i+1) ); n /= (i+1); done = false; break; } } } ans.push( n ); return ans; }</lang>
Logo
<lang logo>to decompose :n [:p 2]
if :p*:p > :n [output (list :n)] if less? 0 modulo :n :p [output (decompose :n bitor 1 :p+1)] output fput :p (decompose :n/:p :p)
end</lang>
Lua
The code of the used auxiliary function "IsPrime(n)" is located at Primality by trial division#Lua
<lang lua>function PrimeDecomposition( n )
local f = {}
if IsPrime( n ) then
f[1] = n
return f
end
local i = 2
repeat
while n % i == 0 do
f[#f+1] = i
n = n / i
end
repeat
i = i + 1
until IsPrime( i )
until n == 1
return f
end</lang>
Mathematica
Bare built-in function does: <lang Mathematica> FactorInteger[2016] => {{2, 5}, {3, 2}, {7, 1}}</lang>
Read as: 2 to the power 5 times 3 squared times 7 (to the power 1). To show them nicely we could use the following functions: <lang Mathematica>supscript[x_,y_]:=If[y==1,x,Superscript[x,y]] ShowPrimeDecomposition[input_Integer]:=Print@@{input," = ",Sequence@@Riffle[supscript@@@FactorInteger[input]," "]}</lang>
Example for small prime: <lang Mathematica> ShowPrimeDecomposition[1337]</lang> gives: <lang Mathematica> 1337 = 7 191</lang> Examples for large primes: <lang Mathematica> Table[AbsoluteTiming[ShowPrimeDecomposition[2^a-1]]//Print[#1," sec"]&,{a,50,150,10}];</lang> gives back: <lang Mathematica>1125899906842623 = 3 11 31 251 601 1801 4051 0.000231 sec 1152921504606846975 = 3^2 5^2 7 11 13 31 41 61 151 331 1321 0.000146 sec 1180591620717411303423 = 3 11 31 43 71 127 281 86171 122921 0.001008 sec 1208925819614629174706175 = 3 5^2 11 17 31 41 257 61681 4278255361 0.000340 sec 1237940039285380274899124223 = 3^3 7 11 19 31 73 151 331 631 23311 18837001 0.000192 sec 1267650600228229401496703205375 = 3 5^3 11 31 41 101 251 601 1801 4051 8101 268501 0.000156 sec 1298074214633706907132624082305023 = 3 11^2 23 31 89 683 881 2971 3191 201961 48912491 0.001389 sec 1329227995784915872903807060280344575 = 3^2 5^2 7 11 13 17 31 41 61 151 241 331 1321 61681 4562284561 0.000374 sec 1361129467683753853853498429727072845823 = 3 11 31 131 2731 8191 409891 7623851 145295143558111 0.024249 sec 1393796574908163946345982392040522594123775 = 3 5^2 11 29 31 41 43 71 113 127 281 86171 122921 7416361 47392381 0.009419 sec 1427247692705959881058285969449495136382746623 = 3^2 7 11 31 151 251 331 601 1801 4051 100801 10567201 1133836730401 0.007705 sec</lang>
MATLAB
<lang Matlab>function [outputPrimeDecomposition] = primedecomposition(inputValue)
outputPrimeDecomposition = factor(inputValue);</lang>
Maxima
Using the built-in function: <lang maxima>(%i1) display2d: false$ /* disable rendering exponents as superscripts */ (%i2) factor(2016); (%o2) 2^5*3^2*7 </lang> Using the underlying language: <lang maxima>prime_dec(n) := apply(append, create_list(makelist(a[1], a[2]), a, ifactors(n)))$
/* or, slighlty more "functional" */ prime_dec(n) := apply(append, map(lambda([a], apply(makelist, a)), ifactors(n)))$
prime_dec(2^4*3^5*5*7^2); /* [2, 2, 2, 2, 3, 3, 3, 3, 3, 5, 7, 7] */</lang>
MUMPS
<lang MUMPS>ERATO1(HI)
SET HI=HI\1 KILL ERATO1 ;Don't make it new - we want it to remain after the quit NEW I,J,P FOR I=2:1:(HI**.5)\1 DO .FOR J=I*I:I:HI DO ..SET P(J)=1 ;$SELECT($DATA(P(J))#10:P(J)+1,1:1) ;WRITE !,"Prime numbers between 2 and ",HI,": " FOR I=2:1:HI DO .S:'$DATA(P(I)) ERATO1(I)=I ;WRITE $SELECT((I<3):"",1:", "),I KILL I,J,P QUIT
PRIMDECO(N)
;Returns its results in the string PRIMDECO ;Kill that before the first call to this recursive function QUIT:N<=1 IF $D(PRIMDECO)=1 SET PRIMDECO="" D ERATO1(N) SET N=N\1,I=0 FOR SET I=$O(ERATO1(I)) Q:+I<1 Q:'(N#I) IF I>1 SET PRIMDECO=$S($L(PRIMDECO)>0:PRIMDECO_"^",1:"")_I D PRIMDECO(N/I) ;that is, if I is a factor of N, add it to the string QUIT</lang>
Usage:
USER>K ERATO1,PRIMDECO D PRIMDECO^ROSETTA(31415) W PRIMDECO 5^61^103 USER>K ERATO,PRIMDECO D PRIMDECO^ROSETTA(31318) W PRIMDECO 2^7^2237 USER>K ERATO,PRIMDECO D PRIMDECO^ROSETTA(34) W PRIMDECO 2^17 USER>K ERATO,PRIMDECO D PRIMDECO^ROSETTA(68) W PRIMDECO 2^2^17 USER>K ERATO,PRIMDECO D PRIMDECO^ROSETTA(7) W PRIMDECO 7 USER>K ERATO,PRIMDECO D PRIMDECO^ROSETTA(777) W PRIMDECO 3^7^37
Nimrod
Based on python solution: <lang nimrod>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)) )</lang>
- 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
OCaml
<lang ocaml>open Big_int;;
let prime_decomposition x =
let rec inner c p =
if lt_big_int p (square_big_int c) then
[p]
else if eq_big_int (mod_big_int p c) zero_big_int then
c :: inner c (div_big_int p c)
else
inner (succ_big_int c) p
in
inner (succ_big_int (succ_big_int zero_big_int)) x;;</lang>
Octave
<lang octave>r = factor(120202039393)</lang>
PARI/GP
GP normally returns factored integers as a matrix with the first column representing the primes and the second their exponents. Thus factor(12)==[2,2;3,1] is true. But it's simple enough to convert this to a vector with repetition:
<lang parigp>pd(n)={
my(f=factor(n),v=f[,1]~);
for(i=1,#v,
while(f[i,2]--,
v=concat(v,f[i,1])
)
);
vecsort(v)
};</lang>
Pascal
<lang pascal>Program PrimeDecomposition(output);
type
DynArray = array of integer;
procedure findFactors(n: Int64; var d: DynArray);
var
divisor, next, rest: Int64;
i: integer;
begin
i := 0;
divisor := 2;
next := 3;
rest := n;
while (rest <> 1) do
begin
while (rest mod divisor = 0) do
begin
setlength(d, i+1);
d[i] := divisor;
inc(i);
rest := rest div divisor;
end;
divisor := next;
next := next + 2;
end;
end;
var
factors: DynArray; j: integer;
begin
setlength(factors, 1); findFactors(1023*1024, factors); for j := low(factors) to high(factors) do writeln (factors[j]);
end.</lang> Output:
% ./PrimeDecomposition 2 2 2 2 2 2 2 2 2 2 3 11 31
Perl
Simple-minded trial division: <lang perl>sub prime_factors { my ($n, $d, @out) = (shift, 1); while ($n > 1 && $d++) { $n /= $d, push @out, $d until $n % $d; } @out }
print "@{[prime_factors(1001)]}\n";</lang>
Perl 6
<lang perl6>constant @primes = 2, 3, 5, -> $n is copy {
repeat { $n += 2 } until $n %% none @primes ... { $_ * $_ >= $n }
$n;
} ... *;
sub factors(Int $remainder is copy) {
return 1 if $remainder <= 1;
gather for @primes -> $factor {
if $factor * $factor > $remainder {
take $remainder if $remainder > 1;
last;
}
# How many times can we divide by this prime?
while $remainder %% $factor {
take $factor;
last if ($remainder div= $factor) === 1;
}
}
}
say factors 536870911;</lang>
Output:
233 1103 2089
PicoLisp
The following solution generates a sequence of "trial divisors" (2 3 5 7 11 13 17 19 23 29 31 37 ..), as described by Donald E. Knuth, "The Art of Computer Programming", Vol.2, p.365. <lang PicoLisp>(de factor (N)
(make
(let (D 2 L (1 2 2 . (4 2 4 2 4 6 2 6 .)) M (sqrt N))
(while (>= M D)
(if (=0 (% N D))
(setq M (sqrt (setq N (/ N (link D)))))
(inc 'D (pop 'L)) ) )
(link N) ) ) )
(factor 1361129467683753853853498429727072845823)</lang> Output:
-> (3 11 31 131 2731 8191 409891 7623851 145295143558111)
PL/I
<lang PL/I> test: procedure options (main, reorder);
declare (n, i) fixed binary (31);
get list (n);
put edit ( n, '[' ) (x(1), a);
restart:
if is_prime(n) then
do;
put edit (trim(n), ']' ) (x(1), a);
stop;
end;
do i = n/2 to 2 by -1;
if is_prime(i) then
if (mod(n, i) = 0) then
do;
put edit ( trim(i) ) (x(1), a);
n = n / i;
go to restart;
end;
end;
put edit ( ' ]' ) (a);
is_prime: procedure (n) options (reorder) returns (bit(1));
declare n fixed binary (31); declare i fixed binary (31);
if n < 2 then return ('0'b);
if n = 2 then return ('1'b);
if mod(n, 2) = 0 then return ('0'b);
do i = 3 to sqrt(n) by 2;
if mod(n, i) = 0 then return ('0'b);
end;
return ('1'b);
end is_prime;
end test; </lang> Results from various runs:
1234567 [ 9721 127 ]
32768 [ 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ]
99 [ 11 3 3 ]
9876543 [ 14503 227 3 ]
100 [ 5 5 2 2 ]
9999999 [ 4649 239 3 3 ]
5040 [ 7 5 3 3 2 2 2 2 ]
Prolog
<lang Prolog>prime_decomp(N, L) :- SN is sqrt(N), prime_decomp_1(N, SN, 2, [], L).
prime_decomp_1(1, _, _, L, L) :- !.
% Special case for 2, increment 1 prime_decomp_1(N, SN, D, L, LF) :- ( 0 is N mod D -> Q is N / D, SQ is sqrt(Q), prime_decomp_1(Q, SQ, D, [D |L], LF) ; D1 is D+1, ( D1 > SN -> LF = [N |L] ; prime_decomp_2(N, SN, D1, L, LF) ) ).
% General case, increment 2 prime_decomp_2(1, _, _, L, L) :- !.
prime_decomp_2(N, SN, D, L, LF) :- ( 0 is N mod D -> Q is N / D, SQ is sqrt(Q), prime_decomp_2(Q, SQ, D, [D |L], LF); D1 is D+2, ( D1 > SN -> LF = [N |L] ; prime_decomp_2(N, SN, D1, L, LF) ) ).</lang>
Output : <lang Prolog> ?- time(prime_decomp(9007199254740991, L)). % 138,882 inferences, 0.344 CPU in 0.357 seconds (96% CPU, 404020 Lips) L = [20394401,69431,6361].
?- time(prime_decomp(576460752303423487, L)).
% 2,684,734 inferences, 0.672 CPU in 0.671 seconds (100% CPU, 3995883 Lips) L = [3203431780337,179951].
?- time(prime_decomp(1361129467683753853853498429727072845823, L)).
% 18,080,807 inferences, 7.953 CPU in 7.973 seconds (100% CPU, 2273422 Lips) L = [145295143558111,7623851,409891,8191,2731,131,31,11,3].
</lang>
Pure
<lang pure>factor n = factor 2 n with
factor k n = k : factor k (n div k) if n mod k == 0;
= if n>1 then [n] else [] if k*k>n; = factor (k+1) n if k==2; = factor (k+2) n otherwise; end;</lang>
PureBasic
<lang PureBasic> CompilerIf #PB_Compiler_Debugger
CompilerError "Turn off the debugger if you want reasonable speed in this example."
CompilerEndIf
Define.q
Procedure Factor(Number, List Factors())
Protected I = 3
While Number % 2 = 0
AddElement(Factors())
Factors() = 2
Number / 2
Wend
Protected Max = Number
While I <= Max And Number > 1
While Number % I = 0
AddElement(Factors())
Factors() = I
Number/I
Wend
I + 2
Wend
EndProcedure
Number = 9007199254740991 NewList Factors() time = ElapsedMilliseconds() Factor(Number, Factors()) time = ElapsedMilliseconds()-time S.s = "Factored " + Str(Number) + " in " + StrD(time/1000, 2) + " seconds." ForEach Factors()
S + #CRLF$ + Str(Factors())
Next MessageRequester("", S)</lang>
Factored 9007199254740991 in 0.27 seconds. 6361 69431 20394401
Python
Note: the program below is imported as a library here.
<lang python> import sys
def is_prime(n):
return zip((True, False), decompose(n))[-1][0]
class IsPrimeCached(dict):
def __missing__(self, n):
r = is_prime(n)
self[n] = r
return r
is_prime_cached = IsPrimeCached()
def primes():
yield 2
n = 3
while n < sys.maxint - 2:
yield n
n += 2
while n < sys.maxint - 2 and not is_prime_cached[n]:
n += 2
def decompose(n):
for p in primes():
if p*p > n: break
while n % p == 0:
yield p
n /=p
if n > 1:
yield n
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 = {0:d}, with factors:".format(m, p) )
start = time.time()
for factor in decompose(p):
print factor,
sys.stdout.flush()
print( "=> {0:.2f}s".format( time.time()-start ) )
if m >= 59:
break
</lang> Output:
2**2-1 = 3, with factors: 3 => 0.00s 2**3-1 = 7, with factors: 7 => 0.00s 2**5-1 = 31, with factors: 31 => 0.00s 2**7-1 = 127, with factors: 127 => 0.00s 2**11-1 = 2047, with factors: 23 89 => 0.00s 2**13-1 = 8191, with factors: 8191 => 0.00s 2**17-1 = 131071, with factors: 131071 => 0.00s 2**19-1 = 524287, with factors: 524287 => 0.01s 2**23-1 = 8388607, with factors: 47 178481 => 0.00s 2**29-1 = 536870911, with factors: 233 1103 2089 => 0.01s 2**31-1 = 2147483647, with factors: 2147483647 => 1.67s 2**37-1 = 137438953471, with factors: 223 616318177 => 0.02s 2**41-1 = 2199023255551, with factors: 13367 164511353 => 0.01s 2**43-1 = 8796093022207, with factors: 431 9719 2099863 => 0.01s 2**47-1 = 140737488355327, with factors: 2351 4513 13264529 => 0.00s 2**53-1 = 9007199254740991, with factors: 6361 69431 20394401 => 1.17s 2**59-1 = 576460752303423487, with factors: 179951 3203431780337 => 211.07s
Note: Python took 740,238 BogoMI to complete the example.
Modifying the primes() and is_prime() functions as below increases performance. <lang python> primelist = [2, 3] def is_prime(n):
if n in primelist: return True
if n < primelist[-1]: return False
for y in primes():
if not n % y: return False
if n < y * y: return True
def primes():
for n in primelist: yield n
n = primelist[-1]
while True:
n += 2
for x in primelist:
if not n % x: break
if x * x > n:
primelist.append(n)
yield n
break
</lang>
Here a shorter and generally way faster algorithm: <lang python>def fac(n):
step = lambda x: 1 + x*4 - (x/2)*2
maxq = long(math.floor(math.sqrt(n)))
d = 1
q = n % 2 == 0 and 2 or 3
while q <= maxq and n % q != 0:
q = step(d)
d += 1
res = []
if q <= maxq:
res.extend(fac(n//q))
res.extend(fac(q))
else: res=[n]
return res
if __name__ == '__main__':
import time start = time.time() tocalc = 2**59-1 print "%s = %s" % (tocalc, fac(tocalc)) print "Needed %ss" % (time.time() - start)</lang>
Output:
576460752303423487 = [3203431780337L, 179951] Needed 0.621000051498s
R
<lang R>findfactors <- function(n) {
d <- c()
div <- 2; nxt <- 3; rest <- n
while( rest != 1 ) {
while( rest%%div == 0 ) {
d <- c(d, div)
rest <- floor(rest / div)
}
div <- nxt
nxt <- nxt + 2
}
d
}
print(findfactors(1005025))</lang>
Racket
<lang Racket>
- lang racket
(require math) (define (factors n)
(append-map (λ (x) (make-list (cadr x) (car x))) (factorize n)))
</lang>
Or, an explicit (and less efficient) computation: <lang Racket>
- lang racket
(define (factors number)
(let loop ([n number] [i 2])
(if (= n 1)
'()
(let-values ([(q r) (quotient/remainder n i)])
(if (zero? r) (cons i (loop q i)) (loop n (add1 i)))))))
</lang>
REXX
No (error) checking was done for the input arguments to test their validity. <lang rexx>/*REXX program fins the prime factors of a (or some) positive integer(s)*/ numeric digits 100 /*bump up precision of the nums. */ parse arg low high . /*get the argument(s). */ if low== then low=1 /*no LOW? Then make one up. */ if high== then high=low /*no HIGH? Then make one up. */ w=length(high) /*get max width for pretty tell. */
do n=low to high /*process single number | a range*/
say right(n,w) 'prime factors =' factr(n)
end /*n*/
exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────FACTR subroutine────────────────────*/ factr: procedure; parse arg x 1 z,,list /*sets X&Z to arg1, LIST to null*/ if x <1 then return /*Too small? Then return null.*/ if x==1 then return 1 /*special case for unity. */
do j=2 to 5; call buildF; end /*fast builds for list: 2 ──► 5. */
j=5 /*start were we left off (J=5). */
do y=0 by 2; j=j+2+y//4 /*insure it's not divisible by 3.*/
if right(j,1)==5 then iterate /*fast check for divisible by 5.*/
if j>z then leave /*num. reduced to a small number?*/
if j*j>x then leave /*are we higher than the √ of X ?*/
call buildF /*add a prime factor to list (J).*/
end /*y*/
if z==1 then return strip(list) /*if residual=unity, don't append*/
return strip(list z) /*return list, append residual. */
/*──────────────────────────────────BUILDF subroutine───────────────────*/ buildF: do forever /*keep dividing until it hurts. */
if z//j\==0 then return /*can't divide any more? */
list=list j /*add number to the list (J). */
z=z%j /*do an integer divide. */
end /*forever*/</lang>
output when the arguments specified are: 1 150
1 prime factors = 1 2 prime factors = 2 3 prime factors = 3 4 prime factors = 2 2 5 prime factors = 5 6 prime factors = 2 3 7 prime factors = 7 8 prime factors = 2 2 2 9 prime factors = 3 3 10 prime factors = 2 5 11 prime factors = 11 12 prime factors = 2 2 3 13 prime factors = 13 14 prime factors = 2 7 15 prime factors = 3 5 16 prime factors = 2 2 2 2 17 prime factors = 17 18 prime factors = 2 3 3 19 prime factors = 19 20 prime factors = 2 2 5 21 prime factors = 3 7 22 prime factors = 2 11 23 prime factors = 23 24 prime factors = 2 2 2 3 25 prime factors = 5 5 26 prime factors = 2 13 27 prime factors = 3 3 3 28 prime factors = 2 2 7 29 prime factors = 29 30 prime factors = 2 3 5 31 prime factors = 31 32 prime factors = 2 2 2 2 2 33 prime factors = 3 11 34 prime factors = 2 17 35 prime factors = 5 7 36 prime factors = 2 2 3 3 37 prime factors = 37 38 prime factors = 2 19 39 prime factors = 3 13 40 prime factors = 2 2 2 5 41 prime factors = 41 42 prime factors = 2 3 7 43 prime factors = 43 44 prime factors = 2 2 11 45 prime factors = 3 3 5 46 prime factors = 2 23 47 prime factors = 47 48 prime factors = 2 2 2 2 3 49 prime factors = 7 7 50 prime factors = 2 5 5 51 prime factors = 3 17 52 prime factors = 2 2 13 53 prime factors = 53 54 prime factors = 2 3 3 3 55 prime factors = 5 11 56 prime factors = 2 2 2 7 57 prime factors = 3 19 58 prime factors = 2 29 59 prime factors = 59 60 prime factors = 2 2 3 5 61 prime factors = 61 62 prime factors = 2 31 63 prime factors = 3 3 7 64 prime factors = 2 2 2 2 2 2 65 prime factors = 5 13 66 prime factors = 2 3 11 67 prime factors = 67 68 prime factors = 2 2 17 69 prime factors = 3 23 70 prime factors = 2 5 7 71 prime factors = 71 72 prime factors = 2 2 2 3 3 73 prime factors = 73 74 prime factors = 2 37 75 prime factors = 3 5 5 76 prime factors = 2 2 19 77 prime factors = 7 11 78 prime factors = 2 3 13 79 prime factors = 79 80 prime factors = 2 2 2 2 5 81 prime factors = 3 3 3 3 82 prime factors = 2 41 83 prime factors = 83 84 prime factors = 2 2 3 7 85 prime factors = 5 17 86 prime factors = 2 43 87 prime factors = 3 29 88 prime factors = 2 2 2 11 89 prime factors = 89 90 prime factors = 2 3 3 5 91 prime factors = 7 13 92 prime factors = 2 2 23 93 prime factors = 3 31 94 prime factors = 2 47 95 prime factors = 5 19 96 prime factors = 2 2 2 2 2 3 97 prime factors = 97 98 prime factors = 2 7 7 99 prime factors = 3 3 11 100 prime factors = 2 2 5 5 101 prime factors = 101 102 prime factors = 2 3 17 103 prime factors = 103 104 prime factors = 2 2 2 13 105 prime factors = 3 5 7 106 prime factors = 2 53 107 prime factors = 107 108 prime factors = 2 2 3 3 3 109 prime factors = 109 110 prime factors = 2 5 11 111 prime factors = 3 37 112 prime factors = 2 2 2 2 7 113 prime factors = 113 114 prime factors = 2 3 19 115 prime factors = 5 23 116 prime factors = 2 2 29 117 prime factors = 3 3 13 118 prime factors = 2 59 119 prime factors = 7 17 120 prime factors = 2 2 2 3 5 121 prime factors = 11 11 122 prime factors = 2 61 123 prime factors = 3 41 124 prime factors = 2 2 31 125 prime factors = 5 5 5 126 prime factors = 2 3 3 7 127 prime factors = 127 128 prime factors = 2 2 2 2 2 2 2 129 prime factors = 3 43 130 prime factors = 2 5 13 131 prime factors = 131 132 prime factors = 2 2 3 11 133 prime factors = 7 19 134 prime factors = 2 67 135 prime factors = 3 3 3 5 136 prime factors = 2 2 2 17 137 prime factors = 137 138 prime factors = 2 3 23 139 prime factors = 139 140 prime factors = 2 2 5 7 141 prime factors = 3 47 142 prime factors = 2 71 143 prime factors = 11 13 144 prime factors = 2 2 2 2 3 3 145 prime factors = 5 29 146 prime factors = 2 73 147 prime factors = 3 7 7 148 prime factors = 2 2 37 149 prime factors = 149 150 prime factors = 2 3 5 5
output when the argument specified is: 9007199254740991
9007199254740991 prime factors = 6361 69431 20394401
Ruby
Ruby's standard library can do prime division. Its 'mathn' package adds Integer#prime_division.
<lang ruby>irb(main):001:0> require 'mathn' => true irb(main):002:0> 2131447995319.prime_division => [[701, 1], [1123, 2], [2411, 1]]</lang>
Ruby 1.9 moves Integer#prime_division to its 'prime' package, though 'mathn' still works. MRI 1.9 is faster than MRI 1.8 with large integers; MRI 1.9 can decompose 2543821448263974486045199 in a few seconds, but MRI 1.8 might need hours.
<lang ruby>irb(main):001:0> require 'prime' => true irb(main):003:0> 2543821448263974486045199.prime_division => [[701, 1], [1123, 2], [2411, 1], [1092461, 2]]</lang>
Simple algorithm
<lang ruby># Get prime decomposition of integer _i_.
- This routine is terribly inefficient, but elegance rules.
def prime_factors(i)
v = (2..i-1).detect{|j| i % j == 0}
v ? ([v] + prime_factors(i/v)) : [i]
end
- Example: Decompose all possible Mersenne primes up to 2**31-1.
- This may take several minutes to show that 2**31-1 is prime.
(2..31).each do |i|
factors = prime_factors(2**i-1)
puts "2**#{i}-1 = #{2**i-1} = #{factors.join(' * ')}"
end</lang>
... 2**28-1 = 268435455 = 3 * 5 * 29 * 43 * 113 * 127 2**29-1 = 536870911 = 233 * 1103 * 2089 2**30-1 = 1073741823 = 3 * 3 * 7 * 11 * 31 * 151 * 331 2**31-1 = 2147483647 = 2147483647
Faster algorithm
<lang ruby># Get prime decomposition of integer _i_.
- This routine is more efficient than prime_factors,
- and quite similar to Integer#prime_division of MRI 1.9.
def prime_factors_faster(i)
factors = []
check = proc do |p|
while(q, r = i.divmod(p)
r.zero?)
factors << p
i = q
end
end
check[2]
check[3]
p = 5
while p * p <= i
check[p]
p += 2
check[p]
p += 4 # skip multiples of 2 and 3
end
factors << i if i > 1
factors
end
- Example: Decompose all possible Mersenne primes up to 2**70-1.
- This may take several minutes to show that 2**61-1 is prime,
- but 2**62-1 and 2**67-1 are not prime.
(2..70).each do |i|
factors = prime_factors_faster(2**i-1)
puts "2**#{i}-1 = #{2**i-1} = #{factors.join(' * ')}"
end</lang>
... 2**67-1 = 147573952589676412927 = 193707721 * 761838257287 2**68-1 = 295147905179352825855 = 3 * 5 * 137 * 953 * 26317 * 43691 * 131071 2**69-1 = 590295810358705651711 = 7 * 47 * 178481 * 10052678938039 2**70-1 = 1180591620717411303423 = 3 * 11 * 31 * 43 * 71 * 127 * 281 * 86171 * 122921
This benchmark compares the different implementations.
<lang ruby>require 'benchmark' require 'mathn' Benchmark.bm(24) do |x|
[2**25 - 6, 2**35 - 7].each do |i|
puts "#{i} = #{prime_factors_faster(i).join(' * ')}"
x.report(" prime_factors") { prime_factors(i) }
x.report(" prime_factors_faster") { prime_factors_faster(i) }
x.report(" Integer#prime_division") { i.prime_division }
end
end</lang>
With MRI 1.8, prime_factors is slow, Integer#prime_division is fast, and prime_factors_faster is very fast. With MRI 1.9, Integer#prime_division is also very fast.
Scala
<lang Scala>import annotation.tailrec import collection.parallel.mutable.ParSeq
object PrimeFactors extends App {
def factorize(n: Long): List[Long] = {
@tailrec
def factors(tuple: (Long, Long, List[Long], Int)): List[Long] = {
tuple match {
case (1, _, acc, _) => acc
case (n, k, acc, _) if (n % k == 0) => factors((n / k, k, acc ++ ParSeq(k), Math.sqrt(n / k).toInt))
case (n, k, acc, sqr) if (k < sqr) => factors(n, k + 1, acc, sqr)
case (n, k, acc, sqr) if (k >= sqr) => factors((1, k, acc ++ ParSeq(n), 0))
}
}
factors((n, 2, List[Long](), Math.sqrt(n).toInt))
}
def mersenne(p: Int): BigInt = (BigInt(2) pow p) - 1
def sieve(nums: Stream[Int]): Stream[Int] = Stream.cons(nums.head, sieve((nums.tail) filter (_ % nums.head != 0))) // An infinite stream of primes, lazy evaluation and memo-ized val oddPrimes = sieve(Stream.from(3, 2)) def primes = sieve(2 #:: oddPrimes)
oddPrimes takeWhile (_ <= 59) foreach { p =>
{ // Needs some intermediate results for nice formatting
val numM = s"M${p}"
val nMersenne = mersenne(p).toLong
val lit = f"${nMersenne}%30d"
val datum = System.nanoTime
val result = factorize(nMersenne)
val mSec = ((System.nanoTime - datum) / 1.e+6).round
def decStr = { if (lit.length > 30) f"(M has ${lit.length}%3d dec)" else "" }
def sPrime = { if (result.isEmpty) " is a Mersenne prime number." else "" }
println(
f"$numM%4s = 2^$p%03d - 1 = ${lit}%s${sPrime} ($mSec%,4d msec) composed of ${result.mkString(" × ")}")
}
}
}</lang>
- Output:
M3 = 2^003 - 1 = 7 ( 23 msec) composed of 7 M5 = 2^005 - 1 = 31 ( 0 msec) composed of 31 M7 = 2^007 - 1 = 127 ( 0 msec) composed of 127 M11 = 2^011 - 1 = 2047 ( 0 msec) composed of 23 × 89 M13 = 2^013 - 1 = 8191 ( 0 msec) composed of 8191 M17 = 2^017 - 1 = 131071 ( 1 msec) composed of 131071 M19 = 2^019 - 1 = 524287 ( 1 msec) composed of 524287 M23 = 2^023 - 1 = 8388607 ( 1 msec) composed of 47 × 178481 M29 = 2^029 - 1 = 536870911 ( 2 msec) composed of 233 × 1103 × 2089 M31 = 2^031 - 1 = 2147483647 ( 39 msec) composed of 2147483647 M37 = 2^037 - 1 = 137438953471 ( 8 msec) composed of 223 × 616318177 M41 = 2^041 - 1 = 2199023255551 ( 2 msec) composed of 13367 × 164511353 M43 = 2^043 - 1 = 8796093022207 ( 2 msec) composed of 431 × 9719 × 2099863 M47 = 2^047 - 1 = 140737488355327 ( 2 msec) composed of 2351 × 4513 × 13264529 M53 = 2^053 - 1 = 9007199254740991 ( 7 msec) composed of 6361 × 69431 × 20394401 M59 = 2^059 - 1 = 576460752303423487 ( 152 msec) composed of 179951 × 3203431780337
Getting the prime factors does not require identifying prime numbers. Since the problems seems to ask for it, here is one version that does it:
<lang Scala>class PrimeFactors(n: BigInt) extends Iterator[BigInt] {
val zero = BigInt(0) val one = BigInt(1) val two = BigInt(2) def isPrime(n: BigInt) = n.isProbablePrime(10) var currentN = n var prime = two
def nextPrime =
if (prime == two) {
prime += one
} else {
prime += two
while (!isPrime(prime)) {
prime += two
if (prime * prime > currentN)
prime = currentN
}
}
def next = {
if (!hasNext)
throw new NoSuchElementException("next on empty iterator")
while(currentN % prime != zero) {
nextPrime
}
currentN /= prime
prime
}
def hasNext = currentN != one && currentN > zero
}</lang>
The method isProbablePrime(n) has a chance of 1 - 1/(2^n) of correctly identifying a prime. Next is a version that does not depend on identifying primes, and works with arbitrary integral numbers:<lang Scala>class PrimeFactors[N](n: N)(implicit num: Integral[N]) extends Iterator[N] {
import num._ val two = one + one var currentN = n var divisor = two
def next = {
if (!hasNext)
throw new NoSuchElementException("next on empty iterator")
while(currentN % divisor != zero) {
if (divisor == two)
divisor += one
else
divisor += two
if (divisor * divisor > currentN)
divisor = currentN
}
currentN /= divisor
divisor
}
def hasNext = currentN != one && currentN > zero
}</lang>
- Output:
Both versions can be rather slow, as they accept arbitrarily big numbers, as requested. Test:
scala> BigInt(2) to BigInt(30) filter (_ isProbablePrime 10) map (p => (p, BigInt(2).pow(p.toInt) - 1)) foreach {
| case (prime, n) => println("2**"+prime+"-1 = "+n+", with factors: "+new PrimeFactors(n).mkString(", "))
| }
2**2-1 = 3, with factors: 3
2**3-1 = 7, with factors: 7
2**5-1 = 31, with factors: 31
2**7-1 = 127, with factors: 127
2**11-1 = 2047, with factors: 23, 89
2**13-1 = 8191, with factors: 8191
2**17-1 = 131071, with factors: 131071
2**19-1 = 524287, with factors: 524287
2**23-1 = 8388607, with factors: 47, 178481
2**29-1 = 536870911, with factors: 233, 1103, 2089
2**31-1 = 2147483647, with factors: 2147483647
2**37-1 = 137438953471, with factors: 223, 616318177
2**41-1 = 2199023255551, with factors: 13367, 164511353
2**43-1 = 8796093022207, with factors: 431, 9719, 2099863
2**47-1 = 140737488355327, with factors: 2351, 4513, 13264529
2**53-1 = 9007199254740991, with factors: 6361, 69431, 20394401
2**59-1 = 576460752303423487, with factors: 179951, 3203431780337
Scheme
<lang scheme>(define (factor number)
(define (*factor divisor number)
(if (> (* divisor divisor) number)
(list number)
(if (= (modulo number divisor) 0)
(cons divisor (*factor divisor (/ number divisor)))
(*factor (+ divisor 1) number))))
(*factor 2 number))
(display (factor 111111111111)) (newline)</lang> Output:
(3 7 11 13 37 101 9901)
Seed7
<lang seed7>const func array integer: factorise (in var integer: number) is func
result
var array integer: result is 0 times 0;
local
var integer: checker is 2;
begin
while checker * checker <= number do
if number rem checker = 0 then
result &:= [](checker);
number := number div checker;
else
incr(checker);
end if;
end while;
if number <> 1 then
result &:= [](number);
end if;
end func;</lang>
Original source: [1]
Slate
Admittedly, this is just based on the Smalltalk entry below: <lang slate>n@(Integer traits) primesDo: block "Decomposes the Integer into primes, applying the block to each (in increasing order)." [| div next remaining |
div: 2.
next: 3.
remaining: n.
[[(remaining \\ div) isZero]
whileTrue:
[block applyTo: {div}.
remaining: remaining // div].
remaining = 1] whileFalse:
[div: next.
next: next + 2] "Just look at the next odd integer."
].</lang>
Smalltalk
<lang smalltalk>Integer extend [
primesDo: aBlock [
| div next rest |
div := 2. next := 3.
rest := self.
[ [ rest \\ div == 0 ]
whileTrue: [
aBlock value: div.
rest := rest // div ].
rest = 1] whileFalse: [
div := next. next := next + 2 ]
]
] 123456 primesDo: [ :each | each printNl ]</lang>
Tcl
<lang tcl>namespace eval primes {}
proc primes::reset {} {
variable list [list] variable current_index end
}
namespace eval primes {reset}
proc primes::restart {} {
variable list
variable current_index
if {[llength $list] > 0} {
set current_index 0
}
}
proc primes::is_prime {candidate} {
variable list
if {$candidate in $list} {return true}
foreach prime $list {
if {$candidate % $prime == 0} {
return false
}
if {$prime * $prime > $candidate} {
return true
}
}
while true {
set largest [get_next_prime]
if {$largest * $largest >= $candidate} {
return [is_prime $candidate]
}
}
}
proc primes::get_next_prime {} {
variable list
variable current_index
if {$current_index ne "end"} {
set p [lindex $list $current_index]
if {[incr current_index] == [llength $list]} {
set current_index end
}
return $p
}
switch -exact -- [llength $list] {
0 {set candidate 2}
1 {set candidate 3}
default {
set candidate [lindex $list end]
while true {
incr candidate 2
if {[is_prime $candidate]} break
}
}
}
lappend list $candidate
return $candidate
}
- return the prime factors of a number in a dictionary.
- The keys will be the factors, the value will be the number
- of times the factor divides the given number
- example: 120 = 2**3 * 3 * 5, so
- [primes::factors 120] returns 2 3 3 1 5 1
- so: set prod 1
- dict for {p e} [primes::factors 120] {
- set prod [expr {$prod * $p**$e}]
- }
- expr {$prod == 120} ;# ==> true
proc primes::factors {num} {
restart
set factors [dict create]
for {set i [get_next_prime]} {$i <= $num} {} {
if {$num % $i == 0} {
dict incr factors $i
set num [expr {$num / $i}]
continue
} elseif {$i*$i > $num} {
dict incr factors $num
break
} else {
set i [get_next_prime]
}
}
return $factors
}</lang> Testing <lang tcl>primes::reset foreach m {2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59} {
set n [expr {2**$m - 1}]
catch {time {set f [dict create {*}[primes::factors $n]]} 1} tm
set primes [list]
dict for {p e} $f {lappend primes {*}[lrepeat $e $p]}
puts [format "2**%02d-1 = %-18s = %-22s => %s" $m $n [join $primes *] $tm]
}</lang> Outputs
2**02-1 = 3 = 3 => 20 microseconds per iteration 2**03-1 = 7 = 7 => 16 microseconds per iteration 2**05-1 = 31 = 31 => 27 microseconds per iteration 2**07-1 = 127 = 127 => 33 microseconds per iteration 2**11-1 = 2047 = 23*89 => 43 microseconds per iteration 2**13-1 = 8191 = 8191 => 159 microseconds per iteration 2**17-1 = 131071 = 131071 => 535 microseconds per iteration 2**19-1 = 524287 = 524287 => 911 microseconds per iteration 2**23-1 = 8388607 = 47*178481 => 162 microseconds per iteration 2**29-1 = 536870911 = 233*1103*2089 => 982 microseconds per iteration 2**31-1 = 2147483647 = 2147483647 => 138831 microseconds per iteration 2**37-1 = 137438953471 = 223*616318177 => 5154 microseconds per iteration 2**41-1 = 2199023255551 = 13367*164511353 => 2901 microseconds per iteration 2**43-1 = 8796093022207 = 431*9719*2099863 => 2141 microseconds per iteration 2**47-1 = 140737488355327 = 2351*4513*13264529 => 1102 microseconds per iteration 2**53-1 = 9007199254740991 = 6361*69431*20394401 => 97472 microseconds per iteration 2**59-1 = 576460752303423487 = 179951*3203431780337 => 12664437 microseconds per iteration
TXR
<lang txr>@(next :args) @(do
(defun factor (n)
(if (> n 1)
(for ((max-d (sqrt n))
(d 2))
(t)
((set d (if (evenp d) (+ d 1) (+ d 2))))
(cond ((> d max-d) (return (list n)))
((zerop (mod n d))
(return (cons d (factor (trunc n d))))))))))
@{num /[0-9]+/} @(bind factors @(factor (int-str num 10))) @(output) @num -> {@(rep)@factors, @(last)@factors@(end)} @(end)</lang>
$ txr factor.txr 1139423842450982345
1139423842450982345 -> {5, 19, 37, 12782467, 25359769}
$ txr factor.txr 1
1 -> {}
$ txr factor.txr 2
2 -> {2}
$ txr factor.txr 3
3 -> {3}
$ txr factor.txr 2
2 -> {2}
$ txr factor.txr 3
3 -> {3}
$ txr factor.txr 4
4 -> {2, 2}
$ txr factor.txr 5
5 -> {5}
$ txr factor.txr 6
6 -> {2, 3}
V
like in scheme (using variables) <lang v>[prime-decomposition
[inner [c p] let
[c c * p >]
[p unit]
[ [p c % zero?]
[c c p c / inner cons]
[c 1 + p inner]
ifte]
ifte].
2 swap inner].</lang>
(mostly) the same thing using stack (with out variables) <lang v>[prime-decomposition
[inner
[dup * <]
[pop unit]
[ [% zero?]
[ [p c : [c p c / c]] view i inner cons]
[succ inner]
ifte]
ifte].
2 inner].</lang>
Using it <lang v>|1221 prime-decomposition puts</lang>
=[3 11 37]
XSLT
Let's assume that in XSLT the application of a template is similar to the invocation of a function. So when the following template <lang xml><xsl:stylesheet xmlns:xsl="http://www.w3.org/1999/XSL/Transform" version="1.0">
<xsl:template match="/numbers">
<html>
<body>
-
<xsl:apply-templates />
</body>
</html>
</xsl:template>
<xsl:template match="number">
</xsl:template> <xsl:template match="value" mode="value"> <xsl:apply-templates /> </xsl:template> <xsl:template match="value" mode="factors"> <xsl:call-template name="generate"> <xsl:with-param name="number" select="number(current())" /> <xsl:with-param name="candidate" select="number(2)" /> </xsl:call-template> </xsl:template> <xsl:template name="generate"> <xsl:param name="number" /> <xsl:param name="candidate" /> <xsl:choose> <xsl:when test="$number = 1"></xsl:when> <xsl:when test="$candidate * $candidate > $number"> <xsl:value-of select="$number" /> </xsl:when> <xsl:when test="$number mod $candidate = 0"> <xsl:value-of select="$candidate" /> <xsl:text> </xsl:text> <xsl:call-template name="generate"> <xsl:with-param name="number" select="$number div $candidate" /> <xsl:with-param name="candidate" select="$candidate" /> </xsl:call-template> </xsl:when> <xsl:otherwise> <xsl:choose> <xsl:when test="$candidate = 2"> <xsl:call-template name="generate"> <xsl:with-param name="number" select="$number" /> <xsl:with-param name="candidate" select="$candidate + 1" /> </xsl:call-template> </xsl:when> <xsl:otherwise> <xsl:call-template name="generate"> <xsl:with-param name="number" select="$number" /> <xsl:with-param name="candidate" select="$candidate + 2" /> </xsl:call-template> </xsl:otherwise> </xsl:choose> </xsl:otherwise> </xsl:choose> </xsl:template> </xsl:stylesheet></lang> is applied against the document <lang xml><numbers> <number><value>1</value></number> <number><value>2</value></number> <number><value>4</value></number> <number><value>8</value></number> <number><value>9</value></number> <number><value>255</value></number> </numbers></lang> then the output contains the prime decomposition of each number: <lang html><html> <body>
- Number: 1 Factors:
- Number: 2 Factors: 2
- Number: 4 Factors: 2 2
- Number: 8 Factors: 2 2 2
- Number: 9 Factors: 3 3
- Number: 255 Factors: 3 5 17
</body> </html></lang>
- Programming Tasks
- Prime Numbers
- Arbitrary precision
- ABAP
- ACL2
- Ada
- ALGOL 68
- AutoHotkey
- AWK
- Commodore BASIC
- Befunge
- Burlesque
- C
- C++
- GMP
- C sharp
- Clojure
- Common Lisp
- D
- E
- Erlang
- F Sharp
- FALSE
- Factor
- Forth
- Fortran
- Frink
- GAP
- Go
- Groovy
- Haskell
- Icon
- Unicon
- Icon Programming Library
- J
- Java
- JavaScript
- Logo
- Lua
- Mathematica
- MATLAB
- Maxima
- MUMPS
- Nimrod
- OCaml
- Octave
- PARI/GP
- Pascal
- Perl
- Perl 6
- PicoLisp
- PL/I
- Prolog
- Pure
- PureBasic
- Python
- R
- Racket
- REXX
- Ruby
- Scala
- Scala Implementations
- Scheme
- Seed7
- Slate
- Smalltalk
- Tcl
- TXR
- V
- GUISS/Omit
- XSLT