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

   count in factors
   factors of an integer
   Sieve of Eratosthenes
   primality by trial division
   factors of a Mersenne number
   trial factoring of a Mersenne number
   partition an integer X into N primes
   sequence of primes by Trial Division


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