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

F decompose(BigInt number)
   [BigInt] result
   V n = number
   BigInt i = 2
   L n % i == 0
      result.append(i)
      n I/= i
   i = 3
   L n >= i * i
      L n % i == 0
         result.append(i)
         n I/= i
      i += 2
   I n != 1
      result.append(n)
   R result
L(i) 2..9
   print(decompose(i))
print(decompose(1023 * 1024))
print(decompose(2 * 3 * 5 * 7 * 11 * 11 * 13 * 17))
print(decompose(BigInt(16860167264933) * 179951))

Output:

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:

/* ARM assembly AARCH64 Raspberry PI 3B */
/*  program primeDecomp64.s   */
/*******************************************/
/* Constantes file                         */
/*******************************************/
/* for this file see task include a file in language AArch64 assembly*/
.include "../includeConstantesARM64.inc"
.equ NBFACT,            100
/*******************************************/
/* Structures                               */
/********************************************/
/* structurea area factors  */
    .struct  0
fac_value:                     // factor
    .struct  fac_value + 8
fac_number:                    // number of identical factors
    .struct  fac_number + 8
fac_end:
/*******************************************/
/* Initialized data */
/*******************************************/
.data
szMessStartPgm:            .asciz "Program start \n"
szMessEndPgm:              .asciz "Program normal end.\n"
szMessNotPrime:            .asciz "Not prime.\n"
szMessPrime:               .asciz "Prime\n"
szCarriageReturn:          .asciz "\n"
szSpaces:                  .asciz "  "
szMessNumber:              .asciz " The factors of @ are :\n"
/*******************************************/
/* UnInitialized data                      */
/*******************************************/
.bss
sZoneConv:           .skip 32
.align 4
tbZoneDecom:         .skip fac_end * NBFACT
/*******************************************/
/*  code section                           */
/*******************************************/
.text
.global main
main:                                    // program start
    ldr x0,qAdrszMessStartPgm            // display start message
    bl affichageMess
    ldr x20,qVal
    //mov x20,17
    mov x0,x20
    ldr x1,qAdrtbZoneDecom
    bl decompFact                        // decomposition
    cmp x0,#0
    beq 1f
    mov x2,x0
    mov x0,x20
    ldr x1,qAdrtbZoneDecom
    bl displayFactors                    // display factors
    b 2f
1:
    ldr x0,qAdrszMessPrime              // prime
    bl affichageMess
2:
    ldr x0,qAdrszMessEndPgm            // display end message
    bl affichageMess
100:                                   // standard end of the program
    mov x0,0                           // return code
    mov x8,EXIT                        // request to exit program
    svc 0                              // perform system call
qAdrszMessStartPgm:        .quad szMessStartPgm
qAdrszMessEndPgm:          .quad szMessEndPgm
qAdrszCarriageReturn:      .quad szCarriageReturn
qAdrszMessNotPrime:        .quad szMessNotPrime
qAdrszMessPrime:           .quad szMessPrime
qAdrtbZoneDecom:           .quad tbZoneDecom
//qVal:                      .quad 2 <<31
qVal:                      .quad 1047552       // test not prime
//qVal:                      .quad 1429671721    // test not prime (37811 * 37811)
/******************************************************************/
/*     prime decomposition                                               */
/******************************************************************/
/* x0 contains the number */
/* x1 contains address factors array */
/* REMARK no save register x9-x19   */
decompFact:
    stp x1,lr,[sp,-16]!       // save  registers
    mov x12,x0                // save number
    bl isPrime                // prime ?
    cbnz x0,12f               // yes -> no decomposition
    mov x19,fac_end           // element area size
    mov x18,0                 // raz indice
    mov x16,0                 // prev divisor
    mov x17,0                 // number of identical divisors
    mov x13,2                 // first divisor
2:
    cmp x12,1
    beq 10f
    udiv x14,x12,x13          // division
    msub x15,x14,x13,x12      // remainder = x12 -(x13*x14)
    cbnz x15,5f               // if remainder <> zero  x13 not divisor
    mov x12,x14               // quotient -> new dividende
    cmp x13,x16               // same divisor ?
    beq 4f                    // yes
    cbz x16,3f                // yes it is first divisor ?
    madd x11,x18,x19,x1       // no -> store prev divisor in the area
    str x16,[x11,fac_value]
    str x17,[x11,fac_number]  // and store number of same factor
    add x18,x18,1             // increment indice
    mov x17,0                 // raz number of same factor
3:
    mov x16,x13               // save new divisor
4:
    add x17,x17,1             // increment number of same factor
    mov x0,x12                // the new dividende is prime ?
    bl isPrime
    cbnz x0,10f               // yes
    b 2b                      // else loop
5:                            // divisor is not a factor
    cmp x13,2                 // begin ?
    cinc x13,x13,ne           // if divisor <> 2 add 1
    add x13,x13,1
    b 2b                      // and loop
10:                           // new dividende is prime
    cmp x16,x12               // divisor = dividende ?
    cinc x17,x17,eq           //add 1 if last dividende = diviseur
    madd x11,x18,x19,x1
    str x16,[x11,fac_value]   // store divisor in area
    str x17,[x11,fac_number]  // and store number
    add x18,x18,1             // increment indice
    cmp x16,x12               //store last dividende if <> diviseur
    beq 11f
    madd x11,x18,x19,x1
    str x12,[x11,fac_value]   // sinon stockage dans la table
    mov x17,1
    str x17,[x11,fac_number]  // store 1 in number
    add x18,x18,1
11:
    mov x0,x18                // return nb factors
    b 100f
12:
    mov x0,#0                 // number is prime
    b 100f
100:
    ldp x1,lr,[sp],16         // restaur des  2 registres
    ret                       // retour adresse lr x30
/******************************************************************/
/*     prime decomposition                                               */
/******************************************************************/
/* x0 contains the number */
/* x1 contains address factors array */
/* x2 number of factors */
displayFactors:
    stp x1,lr,[sp,-16]!        // save  registres
    mov x19,fac_end            // element area size
    mov x13,x1                 // save area address
    ldr x1,qAdrsZoneConv       // load zone conversion address
    bl conversion10
    ldr x0,qAdrszMessNumber
    bl strInsertAtCharInc                   // insert result at Second @ character
    bl affichageMess
    mov x9,0                   // indice
1:
    madd x10,x9,x19,x13        // compute address area element
    ldr x0,[x10,fac_value]
    ldr x12,[x10,fac_number]
    bl conversion10            // decimal conversion
2:
    mov x0,x1
    bl affichageMess
    ldr x0,qAdrszSpaces
    bl affichageMess
    subs x12,x12,#1
    bgt 2b
    add x9,x9,1
    cmp x9,x2
    blt 1b
    ldr x0,qAdrszCarriageReturn
    bl affichageMess
100:
    ldp x1,lr,[sp],16          // restaur des  2 registres
    ret                        // retour adresse lr x30
qAdrsZoneConv:          .quad sZoneConv
qAdrszSpaces:           .quad szSpaces
qAdrszMessNumber:       .quad szMessNumber
/******************************************************************/
/*     test if number is prime                                    */
/******************************************************************/
/* x0 contains the number  */
/* x0 return 1 if prime else return 0 */
isPrime:
    stp x1,lr,[sp,-16]!       // save  registers
    cmp x0,1                  // <= 1 ?
    ble 98f
    cmp x0,3                  // 2 and 3 prime
    ble 97f
    tst x0,1                  //  even ?
    beq 98f
    mov x9,3                  // first divisor
1:
    udiv x11,x0,x9
    msub x10,x11,x9,x0        // compute remainder
    cbz x10,98f               // end if zero
    add x9,x9,#2              // increment divisor
    cmp x9,x11                // divisors<=quotient ?
    ble 1b                    // loop
97:
    mov x0,1                  // return prime
    b 100f
98:
    mov x0,0                  // not prime
    b 100f
100:
    ldp x1,lr,[sp],16         // restaur  2 registers
    ret                       // return to address lr x30
/********************************************************/
/*        File Include fonctions                        */
/********************************************************/
/* for this file see task include a file in language AArch64 assembly */
.include "../includeARM64.inc"

Output:

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

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;

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;

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:

with Ada.Text_IO;                 use Ada.Text_IO;
with Unbounded_Unsigneds;         use Unbounded_Unsigneds;
with Unbounded_Unsigneds.Primes;  use Unbounded_Unsigneds.Primes;
with Strings_Edit.Unbounded_Unsigned_Edit;
use  Strings_Edit.Unbounded_Unsigned_Edit;
with Generic_Unbounded_Array;
procedure Prime_Decomposition is
   type Unbounded_Unsigned_Array is
      array (Positive range <>) of Unbounded_Unsigned;
   function Decompose (X : Unbounded_Unsigned)
      return Unbounded_Unsigned_Array is
      package Unbounded_Arrays is
         new Generic_Unbounded_Array
             (  Positive,
                Unbounded_Unsigned,
                Unbounded_Unsigned_Array,
                Zero
             );
      Result : Unbounded_Arrays.Unbounded_Array;
      Count  : Natural            := 0;
      Factor : Unbounded_Unsigned := Two;
      Value  : Unbounded_Unsigned := X;
      Limit  : Unbounded_Unsigned := Sqrt (X);
   begin
      loop
         if Is_Zero (Value mod Factor) then
            loop
               Count := Count + 1;
               Result.Put (Count, Factor);
               Div (Value, Factor);
               exit when not Is_Zero (Value mod Factor);
            end loop;
            if Is_Prime (Value, 10) = Prime then
               Count := Count + 1;
               Result.Put (Count, Value);
               exit;
            end if;
         end if;
         Next_Prime (Factor, 10);
         exit when Factor > Limit;
      end loop;
      if Count = 0 then
         return (1..0 => Zero);
      else
         return Result.Vector (1..Count);
      end if;
   end Decompose;
   procedure Print (A : Unbounded_Unsigned_Array) is
   begin
      for I in A'Range loop
         if I > A'First then
            Put (", ");
         end if;
         Put (Image (A (I)));
      end loop;
      New_Line;
   end Print;
begin
   Print (Decompose (Two * 2 * 3));
   Print (Decompose (Two * 3 * 5 * 7 * 11 * 11 * 13 * 17));
   Print (Decompose (From_Half_Word (233) * 1103 * 2089));
   Print (Decompose (From_Half_Word (431) * 9719 * 2099863));
   Print (Decompose (From_Half_Word (179951) * 16860167264933));
end Prime_Decomposition;

Output:

begin
integer procedure mod(a, b);
  value a, b; integer a, b;
begin
  mod := a - entier(a/b) * b;
end;
comment
  Decompose n into its prime factors and store
  in the array pf, returning the number found.
  If n is prime, it will be stored as the first
  and only factor;
integer procedure primefactors(n, pf);
  value n; integer n; integer array pf;
begin
  integer i, count;
  count := 1;
  i := 2;
  for i := i while (i * i) <= n do
    begin
      if mod(n, i) = 0 then
        begin
          pf[count] := i;
          count := count + 1;
          n := n / i;
        end
      else
        i := i + 1;
    end;
  pf[count] := n;
  primefactors := count;
end;
comment
  exercise the procedure by displaying the prime
  factors of the odd numbers from 77 to 99;
integer i, k, nfound;
integer array factors[1:32];
for i := 77 step 2 until 99 do
  begin
    nfound := primefactors(i, factors);
    outinteger(1,i);
    outstring(1,":  ");
    for k := 1 step 1 until nfound do
      outinteger(1,factors[k]);
    outstring(1,"\n");
  end;
end

Output:

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

BEGIN
INTEGER I, K, NFOUND;
INTEGER ARRAY FACTORS[1:16];
COMMENT - RETURN P MOD Q;
INTEGER FUNCTION MOD (P, Q);
INTEGER P, Q;
BEGIN
    MOD := P - Q * (P / Q);
END;
COMMENT
  FIND THE PRIME FACTORS OF N AND STORE IN THE EXTERNAL
  ARRAY "FACTORS", RETURNING THE NUMBER FOUND. IF N IS
  PRIME, IT WILL BE STORED AS THE FIRST AND ONLY FACTOR;
INTEGER FUNCTION PRIMEFACTORS(N);
INTEGER N;
BEGIN
  INTEGER P, COUNT;
  P := 2;
  COUNT := 1;
  WHILE N >= P * P DO
    BEGIN
      IF MOD(N, P) = 0 THEN
        BEGIN
          FACTORS[COUNT] := P;
          COUNT := COUNT + 1;
          N := N / P;
        END
      ELSE
        P := P + 1;
    END;
  FACTORS[COUNT] := N;
  PRIMEFACTORS := COUNT;
END;
COMMENT -- EXERCISE THE ROUTINE;
FOR I := 77 STEP 2 UNTIL 99 DO
  BEGIN
    WRITE(I,":");
    NFOUND := PRIMEFACTORS(I);
    COMMENT - PRINT OUT THE FACTORS THAT WERE FOUND;
    FOR K := 1 STEP 1 UNTIL NFOUND DO
      BEGIN
        WRITEON(FACTORS[K]);
      END;
  END;
END

Output:

begin % find the prime decompositionmtion of some integers                  %
    % increments n and returns the new value                                %
    integer procedure inc ( integer value result n ) ; begin n := n + 1; n end;
    % divides n by d and returns the result                                 %
    integer procedure over ( integer value result n
                           ; integer value d
                           ) ; begin n := n div d; n end;
    % sets the elements of f to the prime factors of n                      %
    %      the bounds of f should be 0 :: x where x is large enough to hold %
    %      all the factors, f( 0 ) will contain 6he number of factors       %
    procedure decompose ( integer value n; integer array f ( * ) ) ;
    begin
        integer d, v;
        f( 0 ) := 0;
        v      := abs n;
        if v > 0 and v rem 2 = 0 then begin
            f( inc( f( 0 ) ) ) := 2;
            while over( v, 2 ) > 0 and v rem 2 = 0 do f( inc( f( 0 ) ) ) := 2;
        end if_2_divides_v ;
        d := 3;
        while d * d <= v do begin
            if v rem d = 0 then begin
                f( inc( f( 0 ) ) ) := d;
                while over( v, d ) > 0 and v rem d = 0 do f( inc( f( 0 ) ) ) := d;
            end if_d_divides_v ;
            d := d + 2
        end while_d_squared_le_v ;
        if v > 1 then f( inc( f( 0 ) ) ) := v
    end factorise ;
    % some test cases                                                        %
    for n := 0, 1, 7, 31, 127, 2047, 8191, 131071, 524287, 2520, 32767, 8855, 441421750 do begin
       integer array f( 0 :: 20 );
       decompose( n, f );
       write( s_w := 0, n, ": " );
       for fPos := 1 until f( 0 ) do writeon( i_w := 1, s_w := 0, " ", f( fPos ) );
    end for_n ;
end.

Output:

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

Output:

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

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

num := 8388607, output := ""
for i, p in prime_numbers(num)
    output .= p " * "
MsgBox % num " = " Trim(output, " * ")
return

Output:

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

100 PROGRAM PrimeDecomposition
110 REM -(2^31) has most prime factors (31 twos) than other 32-bit signed integer.
120 DIM Facs(0 TO 30)
130 INPUT PROMPT "Enter a number: ": N 
140 CALL CalcFacs(N, Facs, FacsCnt)
150 REM There is at least one factor
160 FOR I = 0 TO FacsCnt - 1 
170    PRINT Facs(I);
180 NEXT I
190 PRINT
200 END
210 REM **
220 EXTERNAL SUB CalcFacs(N, Facs(), FacsCnt)
230 LET N = ABS(N)
240 LET FacsCnt = 0
250 IF N >= 2 THEN    
260    LET I = 2
270    DO WHILE I * I <= N      
280       IF MOD(N, I) = 0 THEN         
290          LET N = INT(N / I)
300          LET Facs(FacsCnt) = I 
310          LET FacsCnt = FacsCnt + 1
320          LET I = 2      
330       ELSE
340          LET I = I + 1
350       END IF
360    LOOP
370    LET Facs(FacsCnt) = N
380    LET FacsCnt = FacsCnt + 1
390 END IF
400 END SUB

Output:

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

REM Prime decomposition
DIM Facs(14)
REM -(2^15) has most prime factors (15 twos) than other 16-bit signed integer.
PRINT "Enter a number";
INPUT N
GOSUB CalcFacs:
FacsCntM1 = FacsCnt - 1
FOR I = 0 TO FacsCntM1
  PRINT Facs(I);
NEXT I
PRINT
END
CalcFacs:
  N = ABS(N)
  FacsCnt = 0
  IF N >= 2 THEN
    I = 2
    SqrI = I * I
    WHILE SqrI <= N
      NModI = N MOD I
      IF NModI = 0 THEN
        N = N / I
        Facs(FacsCnt) = I
        FacsCnt = FacsCnt + 1
        I = 2
      ELSE
        I = I + 1
      ENDIF
      SqrI = I * I
    WEND
    Facs(FacsCnt) = N
    FacsCnt = FacsCnt + 1
  ENDIF
RETURN

Output:

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

define limit = 20, loops = 0
dim list[limit]
input "loops?", loops
for x = 1 to loops
	let n = x
	print n, " : ",
	gosub collectprimefactors
	for y = 0 to c
		if list[y] then
			print list[y], " ",
			let list[y] = 0
		endif
	next y
	print ""
next x
end
sub collectprimefactors
	let c = 0
	do
		if int(n mod 2) = 0 then
			let n = int(n / 2)
			let list[c] = 2
			let c = c + 1
		endif
		wait
	loop int(n mod 2) = 0
	for i = 3 to sqrt(n) step 2
		do
			if int(n mod i) = 0 then
				let n = int(n / i)
				let list[c] = i
				let c = c + 1
			endif
			wait
		loop int(n mod i) = 0
	next i
	if n > 2 then
		let list[c] = n
		let c = c + 1
	endif
return

Output:

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

10 REM PRIME DECOMPOSITION
20 INPUT "ENTER A NUMBER"N
30 PRINT "--------------"
40 LET N=ABS(N)
50 IF N<2 STOP
60 LET I=2
70 IF I*I>N GOTO 150
80 LET M=N-(N/I)*I
90 IF M#0 GOTO 130
100 LET N=N/I
110 PRINT I
120 LET I=2
130 IF M#0 LET I=I+1
140 GOTO 70
150 PRINT N
160 STOP

Output:

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:

' Prime decomposition
DECLARE SUB Test (N AS INTEGER)
DECLARE FUNCTION FacsCnt% (BYVAL N AS INTEGER, Facs() AS INTEGER)
CONST MAXFACINDEX% = 30
Test 2
Test 2520
Test 13
END
FUNCTION FacsCnt% (BYVAL N AS INTEGER, Facs() AS INTEGER)
  DIM I AS INTEGER, Cnt AS INTEGER
  N = ABS(N)
  Cnt = 0
  IF N >= 2 THEN
    I = 2
    DO WHILE I * I <= N
      IF N MOD I = 0 THEN
        N = N \ I
        Facs(Cnt) = I
        Cnt = Cnt + 1
        I = 2
      ELSE
        I = I + 1
      END IF
    LOOP
    Facs(Cnt) = N
    Cnt = Cnt + 1
  END IF
  FacsCnt% = Cnt
END FUNCTION
SUB Test (N AS INTEGER)
  DIM Facs(0 TO MAXFACINDEX%) AS INTEGER
  DIM Cnt AS INTEGER, I AS INTEGER
  PRINT N; "=>";
  Cnt = FacsCnt%(N, Facs())
  FOR I = 0 TO Cnt - 2
    PRINT Facs(I);
  NEXT I
  PRINT Facs(Cnt - 1)
END SUB

Output:

rem - return p mod q
function mod(p, q = integer) = integer
end = p - q * (p/q)
dim integer factors(16)  rem log2(maxint) is sufficiently large
comment
  Find the prime factors of n and store in global array factors
  (arrays cannot be passed as parameters) and return the number
  found. If n is prime, it will be stored as the only factor.
end
function primefactors(n = integer) = integer
  var p, count = integer
  p = 2
  count = 1
  while n >= (p * p) do
     begin
       if mod(n, p) = 0 then
         begin
           factors(count) = p
           count = count + 1
           n = n / p
         end
       else
         p = p + 1
     end
   factors(count) = n
end = count
rem -- exercise the routine by checking odd numbers from 77 to 99
var i, k, nfound = integer
for i = 77 to 99 step 2
  nfound = primefactors(i)
  print i;"; ";
  for k = 1 to nfound
    print factors(k);
  next k
  print
next i
end

Output:

::prgmPREMIER
Disp "FACTEURS PREMIER"
Prompt N
If N<1:Stop
ClrList L1 ,L2
0→K
iPart(√(N))→L
N→M
For(I,2,L)
0→J
While fPart(M/I)=0
J+1→J
M/I→M
End
If J≠0
Then
K+1→K
I→L 1(K)
J→L2(K)
I→Z:prgmVSTR
"   "+Str0→Str1
If J≠1
Then
J→Z:prgmVSTR
Str1+"^"+Str0→Str1
End
Disp Str1
End
If M=1:Stop
End
If M≠1
Then
If M≠N
Then
M→Z:prgmVSTR
"   "+Str0→Str1
Disp Str1
Else
Disp "PREMIER"
End
End
::prgmVSTR
{Z,Z}→L5
{1,2}→L6
LinReg(ax+b)L6,L5,Y ₀
Equ►String(Y₀,Str0)
length(Str0)→O
sub(Str0,4,O-3)→Str0
ClrList L5,L6
DelVar Y

Output:

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

Output:

Function PrimeFactors(n)
    arrP = Split(ListPrimes(n)," ")
    divnum = n
    Do Until divnum = 1
        'The -1 is to account for the null element of arrP
        For i = 0 To UBound(arrP)-1
            If divnum = 1 Then
                Exit For
            ElseIf divnum Mod arrP(i) = 0 Then
                divnum = divnum/arrP(i)
                PrimeFactors = PrimeFactors & arrP(i) & " "
            End If
        Next
    Loop
End Function
Function IsPrime(n)
    If n = 2 Then
        IsPrime = True
    ElseIf n <= 1 Or n Mod 2 = 0 Then
        IsPrime = False
    Else
        IsPrime = True
        For i = 3 To Int(Sqr(n)) Step 2
            If n Mod i = 0 Then
                IsPrime = False
                Exit For
            End If
        Next
    End If
End Function
Function ListPrimes(n)
    ListPrimes = ""
    For i = 1 To n
        If IsPrime(i) Then
            ListPrimes = ListPrimes & i & " "
        End If
    Next
End Function
WScript.StdOut.Write PrimeFactors(CInt(WScript.Arguments(0)))
WScript.StdOut.WriteLine

Output:

@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

& 211p > : 1 - #v_ 25*, @ > 11g:. /    v
                > : 11g %!|
                          > 11g 1+ 11p v
       ^                               <

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

Output:

   > ⋈⟜Factor¨ 1232123+↕4  # Some factored numbers
┌─                           
╵ 1232123 ⟨ 29 42487 ⟩       
  1232124 ⟨ 2 2 3 102677 ⟩   
  1232125 ⟨ 5 5 5 9857 ⟩     
  1232126 ⟨ 2 7 17 31 167 ⟩  
                            ┘

:import std/Combinator .
:import std/List .
:import std/Math .
factors \divs primes
	divs y [[&[[&[[3 ⋅ 3 >? 4 case-1 (=?0 case-2 case-3)]] (quot-rem 2 1)]]]]
		case-1 4 >? (+1) {}4 empty
		case-2 3 : (5 1 (3 : 2))
		case-3 5 4 2
main [factors <$> ({ (+42) → (+50) })]

Output:

blsq ) 12fC
{2 2 3}

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

#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 *)&lambda; 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;
    }
}

Output:

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

typedef unsigned long long int ulong; // define a type that represent the limit (64-bit)
ulong mod_mul(ulong a, ulong b, const ulong mod) {
	ulong res = 0, c; // return (a * b) % mod, avoiding overflow errors while doing modular multiplication.
	for (b %= mod; a; a & 1 ? b >= mod - res ? res -= mod : 0, res += b : 0, a >>= 1, (c = b) >= mod - b ? c -= mod : 0, b += c);
	return res % mod;
}
ulong mod_pow(ulong n, ulong exp, const ulong mod) {
	ulong res = 1; // return (n ^ exp) % mod
	for (n %= mod; exp; exp & 1 ? res = mod_mul(res, n, mod) : 0, n = mod_mul(n, n, mod), exp >>= 1);
	return res;
}
ulong square_root(const ulong N) {
	ulong res = 0, rem = N, c, d;
	for (c = 1 << 62; c; c >>= 2) {
		d = res + c;
		res >>= 1;
		if (rem >= d)
			rem -= d, res += c;
	} // returns the square root of N.
	return res;
}
int is_prime(const ulong N) {
	ulong i = 1; // return a truthy value about the primality of N.
	if (N > 1) for (; i < 64 && mod_pow(i, N - 1, N) <= 1; ++i);
	return i == 64;
}
ulong pollard_rho(const ulong N) {
	// Require : N is a composite number, not a square.
	// Ensure : res is a non-trivial factor of N.
	// Option : change the timeout, change the rand function.
	static const int timeout = 18;
	static unsigned long long rand_val = 2994439072U;
	rand_val = (rand_val * 1025416097U + 286824428U) % 4294967291LLU;
	ulong res = 1, a, b, c, i = 0, j = 1, x = 1, y = 1 + rand_val % (N - 1);
	for (; res == 1; ++i) {
		if (i == j) {
			if (j >> timeout)
				break;
			j <<= 1;
			x = y;
		}
		a = y, b = y; // performs y = (y * y) % N
		for (y = 0; a; a & 1 ? b >= N - y ? y -= N : 0, y += b : 0, a >>= 1, (c = b) >= N - b ? c -= N : 0, b += c);
		y = (1 + y) % N;
		for (a = y > x ? y - x : x - y, b = N; (a %= b) && (b %= a);); // compute the gcd(abs(y - x), N);
		res = a | b;
	}
	return res;
}
void factor(const ulong N, ulong *array) {
	// very basic manager that fill the given array (the size of the result is the first array element)	
	// it does not perform initial trial divisions, which is generally highly recommended.
	if (N < 4 || is_prime(N)) {
		if (N > 1 || !*array) array[++*array] = N;
		return;
	}
	ulong x = square_root(N);
	if (x * x != N) x = pollard_rho(N);
	factor(x, array);
	factor(N / x, array);
}
#include <stdio.h>
int main(void) {
	// simple test.
	unsigned long long n = 18446744073709551615U;
	ulong fac[65] = {0};
	factor(n, fac);
	for (ulong i = 1; i <= *fac; ++i)
		printf("* %llu\n", fac[i]);
}