A basic forest/count-based implementation.

Although forest provides keys to evaluate stuff and has basic if keys, I used a macro that evaluates everything. As it relies on counts, I don’t think it works with forest’s keys as they use pgfmath. (Which suffers from the typical problem that it uses TeX’s dimen registers which cannot be greater than 18-thousand-something; with the fp package and TikZ’s fixedpointarithmetic library this could be extended.)

The greatest representable number is 2 147 483 644 (2^31-4), the next integer 2 147 483 645 (2^31-3) fails.

The algorithm is not very efficient.
For the input number p, the factors 2, 3, 5, 7, 9, …, 2n + 1 until p/2 are checked whether they divide p without remainder.

The code includes some comments on the algorithm.

The \num macro from the siunitx package was used to typeset the numbers (even the exponents even though they cannot be greater than 30).

Code

\documentclass[tikz]{standalone}
\usepackage{forest,mathtools,siunitx}
\makeatletter
\def\ifNum#1{\ifnum#1\relax
  \expandafter\pgfutil@firstoftwo\else
  \expandafter\pgfutil@secondoftwo\fi}
\forestset{
  num content/.style={
    delay={
      content/.expanded={\noexpand\num{\forestoption{content}}}}},
  pt@prime/.style={draw, circle},
  pt@start/.style={},
  pt@normal/.style={},
  start primeTree/.style={%
    /utils/exec=%
      % \pt@start holds the current minimum factor, we'll start with 2
      \def\pt@start{2}%
      % \pt@result will hold the to-be-typeset factorization, we'll start with
      % \pgfutil@gobble since we don't want a initial \times
      \let\pt@result\pgfutil@gobble
      % \pt@start@cnt holds the number of ^factors for the current factor
      \def\pt@start@cnt{0}%
      % \pt@lStart will later hold "l"ast factor used
      \let\pt@lStart\pgfutil@empty,
    alias=pt-start,
    pt@start/.try,
    delay={content/.expanded={$\noexpand\num{\forestove{content}}
                            \noexpand\mathrlap{{}= \noexpand\pt@result}$}},
    primeTree},
  primeTree/.code=%
    % take the content of the node and save it in the count
    \c@pgf@counta\forestove{content}\relax
    % if it's 2 we're already finished with the factorization
    \ifNum{\c@pgf@counta=2}{%
      % add the factor
      \pt@addfactor{2}%
      % finalize the factorization of the result
      \pt@addfactor{}%
      % and set the style to the prime style
      \forestset{pt@prime/.try}%
    }{%
      % this simply calculates content/2 and saves it in \pt@end
      % this is later used for an early break of the recursion since no factor
      % can be greater then content/2 (for integers of course)
      \edef\pt@content{\the\c@pgf@counta}%
      \divide\c@pgf@counta2\relax
      \advance\c@pgf@counta1\relax % to be on the safe side
      \edef\pt@end{\the\c@pgf@counta}%
      \pt@do}}

%%% our main "function"
\def\pt@do{%
  % let's test if the current factor is already greather then the max factor
  \ifNum{\pt@end<\pt@start}{%
    % great, we're finished, the same as above
    \expandafter\pt@addfactor\expandafter{\pt@content}%
    \pt@addfactor{}%
    \def\pt@next{\forestset{pt@prime/.try}}%
  }{%
    % this calculates int(content/factor)*factor
    % if factor is a factor of content (without remainder), the result will
    % equal content. The int(content/factor) is saved in \pgf@temp.
    \c@pgf@counta\pt@content\relax
    \divide\c@pgf@counta\pt@start\relax
    \edef\pgf@temp{\the\c@pgf@counta}%
    \multiply\c@pgf@counta\pt@start\relax
    \ifNum{\the\c@pgf@counta=\pt@content}{%
      % yeah, we found a factor, add it to the result and ...
      \expandafter\pt@addfactor\expandafter{\pt@start}%
      % ... add the factor as the first child with style pt@prime
      % and the result of int(content/factor) as another child.
      \edef\pt@next{\noexpand\forestset{%
        append={[\pt@start, pt@prime/.try]},
        append={[\pgf@temp, pt@normal/.try]},
        % forest is complex, this makes sure that for the second child, the
        % primeTree style is not executed too early (there must be a better way).
        delay={
          for descendants={
            delay={if n'=1{primeTree, num content}{}}}}}}%
    }{%
      % Alright this is not a factor, let's get the next factor
      \ifNum{\pt@start=2}{%
        % if the previous factor was 2, the next one will be 3
        \def\pt@start{3}%
      }{%
        % hmm, the previos factor was not 2,
        % let's add 2, maybe we'll hit the next prime number
        % and maybe a factor
        \c@pgf@counta\pt@start
        \advance\c@pgf@counta2\relax
        \edef\pt@start{\the\c@pgf@counta}%
      }%
      % let's do that again
      \let\pt@next\pt@do
    }%
  }%
  \pt@next
}

%%% this builds the \pt@result macro with the factors
\def\pt@addfactor#1{%
  \def\pgf@tempa{#1}%
  % is it the same factor as the previous one
  \ifx\pgf@tempa\pt@lStart
    % add 1 to the counter
    \c@pgf@counta\pt@start@cnt\relax
    \advance\c@pgf@counta1\relax
    \edef\pt@start@cnt{\the\c@pgf@counta}%
  \else
    % a new factor! Add the previous one to the product of factors
    \ifx\pt@lStart\pgfutil@empty\else
      % as long as there actually is one, the \ifnum makes sure we do not add ^1
      \edef\pgf@tempa{\noexpand\num{\pt@lStart}\ifnum\pt@start@cnt>1 
                                           ^{\noexpand\num{\pt@start@cnt}}\fi}%
      \expandafter\pt@addfactor@\expandafter{\pgf@tempa}%
    \fi
    % setup the macros for the next round
    \def\pt@lStart{#1}% <- current (new) factor
    \def\pt@start@cnt{1}% <- first time
  \fi
}
%%% This simply appends "\times #1" to \pt@result, with etoolbox this would be
%%% \appto\pt@result{\times#1}
\def\pt@addfactor@#1{%
  \expandafter\def\expandafter\pt@result\expandafter{\pt@result \times #1}}

%%% Our main macro:
%%% #1 = possible optional argument for forest (can be tikz too)
%%% #2 = the number to factorize
\newcommand*{\PrimeTree}[2][]{%
  \begin{forest}%
    % as the result is set via \mathrlap it doesn't update the bounding box
    % let's fix this:
    tikz={execute at end scope={\pgfmathparse{width("${}=\pt@result$")}%
                         \path ([xshift=\pgfmathresult pt]pt-start.east);}},
    % other optional arguments
    #1
    % And go!
    [#2, start primeTree]
  \end{forest}}
\makeatother
\begin{document}
\PrimeTree{36}
\PrimeTree{90}
\PrimeTree{112}
\PrimeTree{612}
\PrimeTree{7875}
\PrimeTree{22230}
\PrimeTree{1073741824}
\PrimeTree{2147483644}
\end{document}

Output

edit (2017/09/01 & 02):

Release 1.2o deprecates some macros which were used here, so the answer was updated. On this occasion I replaced uses of \if...\else...\fi with the xint provided \xintiiifOne etc... which safely select a branch à la firstoftwo/secondoftwo way and avoid things like \expandafter\expandafter\expandafter at user level.

Note: the xint manual also contains code implementing expandably Rabin-Miller strong pseudo-primality tests...

This answer has evolved in stages. New contributions were added at the bottom (apart from the images which I put on top).

1. macro \factorize to produce programmatically the factorization of input N. The output is put in macro \factors. For example for N=36, \factors expands to {{}{}{36}}{{2}{2}{9}}{{3}{2}{1}}. Each successive triplet is {p}{k}{m} where p^k is the p-factor of N, and m the result of dividing N by it and all previous ones. The first triple {}{}{N} is a bit of a nuisance, and this explains \expandafter\@gobble\factors in some the displaying codes. Initially the displays used tabulars.

2. then I added code using pst-tree to produce the trees from the result available in \factors. I also added code for just printing the factorization inline.

3. I now add code using TikZ+forest, which I have learned a bit since seeing it used beautifully in Qrrbrbirlbel's answer. The code for the forest tree also incoporates the inline version of the factorization.

I use pst-tree and forest only to display, not to compute the factorization. Both pst-tree and forest syntax for the trees allowed a simple approach to work from \factors to the trees: two token lists are filled-up step by step; if using the native TikZ syntax with child, that would be much more difficult because of braces { and }.

Images of trees produced with forest:

Images of trees using pst-tree (as were done manually by PSTikZ ^[1].) Here are some samples (the code also has the variant to not display the 1's when they are exponents). For some other ways to do trees, e.g with TikZ's childs and nodes, the \factors macro prepared by the \factorize command should probably be done in a slightly different manner.

I too propose half of an answer using package xint ^[2] for the arithmetic computations (on arbitrarily big numbers). The command \factorize{N} results in macro \factors containing a list of triplets {p}{k}{m}, where p is a prime number showing up in the factorization, k is its exponent, m is the result of division of N by p^k as well as all previous factors corresponding to smaller primes. So the last triplet has m=1, and the first is {}{}{N}.

The tree could then be constructed from this macro \factors; here I just have a \displayfactorization to display the result using tabulars (and I use some macros of xint ^[3] to transform \factors into what goes into the tabular). There is an optional argument to set up the width of the second column.

I use the numprint ^[4] package to print very long numbers without going beyond the page physical limits.

The algorithm is very lame: first divisions by 2 are tested, then 3, then 5, and all successive odd integers until N has been reduced to 1.

There is no impact on TeX memory, apart from \factors being filled up.

The algorithm would be of course faster if done using only TeX \count registers (and eTeX \numexpr for convenience), but this would restrict it to numbers < 2^{31}.

update: I have incorporated the usual halting test to not go beyond square root of n, makes the thing a bit more efficient when the factorization ends in a 'big' prime (two such examples added).

\documentclass{article}
\usepackage{xint}
\usepackage{xintexpr}
\usepackage[T1]{fontenc}
\usepackage{array}
\usepackage[np]{numprint}
\AtBeginDocument{\npthousandsep{,\hskip .1pt plus .02pt minus .02pt}}

\catcode`\_ 11

\def\factorize#1{%
    \edef\factorize_N{#1}%
    \def\factorize_exp{0}%
    \edef\factors{{{}{}{\factorize_N}}}%
    \factorize_i
}

\def\factorize_i{%
    \xintiiifOdd{\factorize_N}%
      {\factorize_ii}%
      {\edef\factorize_exp{\xintInc{\factorize_exp}}%
       \edef\factorize_N  {\xintHalf{\factorize_N}}%
       \factorize_i}%
}

\def\factorize_ii{%
    \xintiiifZero{\factorize_exp}%
      {}%
      {\edef\factors{\factors{{2}{\factorize_exp}{\factorize_N}}}}%
    \xintiiifOne{\factorize_N}%
      {}%
      {\def\factorize_M{3}%
       \def\factorize_exp{0}%
       \factorize_iii}%
}

\def\factorize_iii{%
    \xintAssign\xintiiDivision\factorize_N\factorize_M\to\factorize_Q\factorize_R
    \xintiiifSgn{\factorize_R}%
      {% never happens: remainder can not be negative
      }%
      {% case of vanishing remainder
       \edef\factorize_exp{\xintInc{\factorize_exp}}%
       \let\factorize_N\factorize_Q 
       \factorize_iii
      }%
      {\factorize_iv}% 
}

\def\factorize_iv{%
    \xintiiifZero{\factorize_exp}%
      {}%
      {\edef\factors{\factors{{\factorize_M}{\factorize_exp}{\factorize_N}}}}%
    \xintiiifOne{\factorize_N}%
      {}%
      {% here N>1, N=QM+R (0<R<Q) is < M(Q+1) and N has no prime factors
       % at most equal to M. If a prime P>M divides N, the
       % quotient N/P will be < Q+1, hence at most Q. If Q<=M, then
       % N/P must be 1 else there would be some prime <=M dividing N.
       % no \xintiiifGeq ...
       % \xintiifCmp will have branches for each of <, =, >, less convenient
       % So we use \xintiiifLt which exists, and permute the branches
       % compared to original code
       \xintiiifLt\factorize_M\factorize_Q
         {% we go on testing with bigger factors
          % or \edef\factorize_M{\xintInc{\xintInc{\factorize_M}}} perhaps
          \edef\factorize_M{\xintiiAdd \factorize_M 2}%
          \def\factorize_exp{0}%
          \factorize_iii
         }%
         {% implies that N is prime
          \edef\factors{\factors{{\factorize_N}{1}{1}}}% we stop here
         }%
      }%
}

\catcode`\_ 8

\def\auxiliarymacroa #1{\auxiliarymacrob #1}
\def\auxiliarymacrob #1#2#3{${#1}^{#2}$&\np{#3}\tabularnewline\hline}

\newcommand*\displayfactorization[1][.2\linewidth]{%
   \begin{tabular}{>{\rule{0pt}{11pt}}c|>{\raggedright}p{#1}}%
   \xintApplyUnbraced\auxiliarymacroa{\factors}
   \end{tabular}\hskip.5em plus .125em minus .125em }

%for testing:
% \def\displayfactorization{\meaning\factors}

\pagestyle{empty}
\begin{document}\thispagestyle{empty}\ttfamily

\noindent\factorize{36}\displayfactorization
\factorize{90}\displayfactorization
\factorize{112}\displayfactorization
\factorize{612}\displayfactorization
\factorize{7875}\displayfactorization
\factorize{22230}\displayfactorization
\factorize{1073741824}\displayfactorization
\factorize{2147483644}\displayfactorization
\factorize{\xintiiPow 2{40}}\displayfactorization
\factorize{\xinttheiiexpr 2^37 * 3^23 * 17^13\relax}%
\displayfactorization
% two examples ending with a (somewhat) `big' prime,
\factorize{10968058712}\displayfactorization
\factorize{1689242184972}\displayfactorization

\factorize{\xinttheiiexpr 111^13 * 371^7 * 1271^35\relax}
\displayfactorization[.75\linewidth]

% \factorize{2147483647}% does not exhaust memory takes about 2s

% \clearpage

% about 6s total last time I tried :

% \def\factorizeanddisplay #1{%
%     \pdfresettimer
%     \factorize{#1}%
%     \edef\z{\the\pdfelapsedtime}%
%     (needed \xintRound{2}{\z/65536} seconds)
%     \displayfactorization
%     \par\medskip
% }

% \factorizeanddisplay{2147483642}
% \factorizeanddisplay{2147483643}
% \factorizeanddisplay{2147483644}
% \factorizeanddisplay{2147483645}
% \factorizeanddisplay{2147483646}
% \factorizeanddisplay{2147483647}
% \factorizeanddisplay{2147483648}
% \factorizeanddisplay{2147483649}
% \factorizeanddisplay{2147483650}
\end{document}

and an example with big integers:

And now the code repeated, together with code to produce trees in the format of the OP question, using pst-tree:

This resuires latex+dvips or can also use xelatex.

% COMPILE WITH LATEX+DVIPS

\documentclass[border=3pt,preview,varwidth,multi]{standalone}

\usepackage{pst-tree}
\psset{levelsep=1,treesep=1,nodesep=2pt}

\usepackage{xint}
\usepackage{xintexpr}

\catcode`\_ 11

% This code (non-expandable) produces
% {{}{}{N}} followed with successive braced triplets
% {{p}{k}{m}} where p is a prime factor of N, 
% k its exponent in N, and m is the result of dividing
% N by p^k and all previous powers of smaller primes
% So the last triplet has m=1

% The code uses package xint to be able to deal
% with numbers larger than the TeX limit of 2^{31}-1
% on count registers. 

\def\factorize#1{%
    \edef\factorize_N{#1}%
    \def\factorize_exp{0}%
    \edef\factors{{{}{}{\factorize_N}}}%
    \factorize_i
}

\def\factorize_i{%
    \xintiiifOdd{\factorize_N}%
      {\factorize_ii}%
      {\edef\factorize_exp{\xintInc{\factorize_exp}}%
       \edef\factorize_N  {\xintHalf{\factorize_N}}%
       \factorize_i}%
}

\def\factorize_ii{%
    \xintiiifZero{\factorize_exp}%
      {}%
      {\edef\factors{\factors{{2}{\factorize_exp}{\factorize_N}}}}%
    \xintiiifOne{\factorize_N}%
      {}%
      {\def\factorize_M{3}%
       \def\factorize_exp{0}%
       \factorize_iii}%
}

\def\factorize_iii{%
    \xintAssign\xintiiDivision\factorize_N\factorize_M\to\factorize_Q\factorize_R
    \xintiiifSgn{\factorize_R}%
      {% never happens: remainder can not be negative
      }%
      {% case of vanishing remainder
       \edef\factorize_exp{\xintInc{\factorize_exp}}%
       \let\factorize_N\factorize_Q 
       \factorize_iii
      }%
      {\factorize_iv}% 
}

\def\factorize_iv{%
    \xintiiifZero{\factorize_exp}%
      {}%
      {\edef\factors{\factors{{\factorize_M}{\factorize_exp}{\factorize_N}}}}%
    \xintiiifOne{\factorize_N}%
      {}%
      {% here N>1, N=QM+R (0<R<Q) is < M(Q+1) and N has no prime factors
       % at most equal to M. If a prime P>M divides N, the
       % quotient N/P will be < Q+1, hence at most Q. If Q<=M, then
       % N/P must be 1 else there would be some prime <=M dividing N.
       % no \xintiiifGeq ...
       % \xintiifCmp will have branches for each of <, =, >, less convenient
       % So we use \xintiiifLt which exists, and permute the branches
       % compared to original code
       \xintiiifLt\factorize_M\factorize_Q
         {% we go on testing with bigger factors
          % or \edef\factorize_M{\xintInc{\xintInc{\factorize_M}}} perhaps
          \edef\factorize_M{\xintiiAdd \factorize_M 2}%
          \def\factorize_exp{0}%
          \factorize_iii
         }%
         {% implies that N is prime
          \edef\factors{\factors{{\factorize_N}{1}{1}}}% we stop here
         }%
      }%
}

\catcode`\_ 8

%% We now define the macro \FactorTree which will produce a tree
%% displaying the factorization, using PSTricks


\newtoks\FactorTreeA
\newtoks\FactorTreeB

\makeatletter

\newcommand*{\FactorsToTree}[1]{%
    \FactorsToTree@ #1%
}

% macro which was used to produce the images; variant follows which skips the
% exponents equal to 1.

% \newcommand*{\FactorsToTree@}[3]{%
%     \xintSgnFork{\xintCmp{#3}{1}}% check to see if end has been reached
%     {}%
%     {\FactorTreeA\expandafter{\the\FactorTreeA
%                               \Tcircle{${#1}^{#2}$}%
%                               \TR{1}%
%                               }}%
%     {\FactorTreeA\expandafter{\the\FactorTreeA 
%                              \Tcircle{${#1}^{#2}$}%
%                              \psTree{\TR{#3}}}%
%      \FactorTreeB\expandafter{\the\FactorTreeB \endpsTree}}%
% }


% This variant will not print the exponents equal to 1:

\newcommand*{\FactorsToTree@}[3]{%
    \ifnum 0#2=1 % first triplet has an empty #2, hence the trick with 0
       \expandafter\@firstoftwo
    \else
       \expandafter\@secondoftwo
    \fi
    % exponent #2 is 1, so don't print it
    {\xintSgnFork{\xintCmp{#3}{1}}% check to see if end has been reached
    {}%
    {\FactorTreeA\expandafter{\the\FactorTreeA
                              \Tcircle{${#1}$}%
                              \TR{1}%
                              }}%
    {\FactorTreeA\expandafter{\the\FactorTreeA 
                             \Tcircle{${#1}$}%
                             \psTree{\TR{#3}}}%
     \FactorTreeB\expandafter{\the\FactorTreeB \endpsTree}}}
    % exponent #2 is > 1 (or absent in the {}{}{N} triplet)
    {\xintSgnFork{\xintCmp{#3}{1}}% check to see if end has been reached
    {}%
    {\FactorTreeA\expandafter{\the\FactorTreeA
                              \Tcircle{${#1}^{#2}$}%
                              \TR{1}%
                              }}%
    {\FactorTreeA\expandafter{\the\FactorTreeA 
                             \Tcircle{${#1}^{#2}$}%
                             \psTree{\TR{#3}}}%
     \FactorTreeB\expandafter{\the\FactorTreeB \endpsTree}}}%
}


\newcommand*{\FactorTree}[1]{%
    \factorize{#1}%
    \FactorTreeA{\@gobbletwo}%
    \FactorTreeB{}%
    \xintApplyUnbraced\FactorsToTree{\factors}%
    \the\FactorTreeA\the\FactorTreeB
}

\makeatother


\begin{document}

%% \preview\FactorTree{1}\endpreview  %% (ok)

\preview\FactorTree{13}\endpreview

\preview\FactorTree{36}\endpreview

\preview\FactorTree{90}\endpreview

\preview\FactorTree{112}\endpreview

\preview\FactorTree{612}\endpreview

\preview\FactorTree{7875}\endpreview

\preview\FactorTree{22230}\endpreview

\preview\FactorTree{1073741824}\endpreview

\preview\FactorTree{2147483644}\endpreview

\preview\FactorTree{\xintiiPow 2{40}}\endpreview

\preview\FactorTree{\xinttheiiexpr 2^37 * 3^23 * 17^13\relax}%
\endpreview
% two examples ending with a (somewhat) `big' prime,

\preview\FactorTree{10968058712}\endpreview

\preview\FactorTree{1689242184972}\endpreview

\end{document}

Code for inline product decomposition:

% The command \FactorizeInline produces (for math mode, at it uses \times) the
% product decomposition of its argument into prime powers, the exponents equal
% to 1 are not printed. The argument is supposed to be an integer > 1
% (arbitrarily big, but computation times will be large if it is not a product
% of reasonably small primes).
% $N = \FactorizeInline{N}$

\makeatletter
\def\@factorinliner  #1{\@factorinliner@ #1}
\def\@factorinliner@ #1#2#3{\ifnum #2>1 \expandafter\@firstoftwo\else
                                        \expandafter\@secondoftwo\fi
                            {{#1}^{#2}}{#1}}
\newcommand*{\FactorizeInline}[1]{%
    \factorize{#1}% 
    \xintListWithSep\times
            {\xintApply\@factorinliner{\expandafter\@gobble\factors}}%
}%
\makeatother

$13=\FactorizeInline{13}$

$36=\FactorizeInline{36}$

$90=\FactorizeInline{90}$

$112=\FactorizeInline{112}$

$612=\FactorizeInline{612}$

$7875=\FactorizeInline{7875}$

$22230=\FactorizeInline{22230}$

$1073741824=\FactorizeInline{1073741824}$

$2147483642=\FactorizeInline{2147483642}$

% $2147483643=\FactorizeInline{2147483643}$ % 4.5 seconds on my laptop

$2147483644=\FactorizeInline{2147483644}$

$2147483645=\FactorizeInline{2147483645}$

% $2147483646=\FactorizeInline{2147483646}$

% $2147483647=\FactorizeInline{2147483647}$ % 8 seconds on my 2012 laptop

$2147483648=\FactorizeInline{2147483648}$

% $2147483649=\FactorizeInline{2147483649}$ % 4.5 seconds on my laptop

$2147483650=\FactorizeInline{2147483650}$

% $19928968819=\FactorizeInline{19928968819}$ % 25 seconds on my 2012 laptop

$\xintiiPow 2{40} = \FactorizeInline{\xintiiPow 2{40}}$

% two examples ending with a (somewhat) `big' prime,

$10968058712=\FactorizeInline{10968058712}$

$1689242184972=\FactorizeInline{1689242184972}$

%\xinttheiiexpr 2^37 * 3^23 * 17^13\relax

$128154740640622513993964746937443811328=\FactorizeInline{128154740640622513993964746937443811328}$

Code for doing the trees with forest:

\documentclass[border=3pt,varwidth,multi={forest}]{standalone}

\usepackage{calc} % for \widthof

\usepackage{tikz}
\usetikzlibrary{calc}
\usepackage{forest}

\usepackage{xint}
\usepackage{xintexpr}


%  macro \factorize as above
\catcode`\_ 11

% This code (non-expandable) produces
% {{}{}{N}} followed with successive braced triplets
% {{p}{k}{m}} where p is a prime factor of N, 
% k its exponent in N, and m is the result of dividing
% N by p^k and all previous powers of smaller primes
% So the last triplet has m=1

% The code uses package xint to be able to deal
% with numbers larger than the TeX limit of 2^{31}-1
% on count registers. 


\def\factorize#1{%
    \edef\factorize_N{#1}%
    \def\factorize_exp{0}%
    \edef\factors{{{}{}{\factorize_N}}}%
    \factorize_i
}

\def\factorize_i{%
    \xintiiifOdd{\factorize_N}%
      {\factorize_ii}%
      {\edef\factorize_exp{\xintInc{\factorize_exp}}%
       \edef\factorize_N  {\xintHalf{\factorize_N}}%
       \factorize_i}%
}

\def\factorize_ii{%
    \xintiiifZero{\factorize_exp}%
      {}%
      {\edef\factors{\factors{{2}{\factorize_exp}{\factorize_N}}}}%
    \xintiiifOne{\factorize_N}%
      {}%
      {\def\factorize_M{3}%
       \def\factorize_exp{0}%
       \factorize_iii}%
}

\def\factorize_iii{%
    \xintAssign\xintiiDivision\factorize_N\factorize_M\to\factorize_Q\factorize_R
    \xintiiifSgn{\factorize_R}%
      {% never happens: remainder can not be negative
      }%
      {% case of vanishing remainder
       \edef\factorize_exp{\xintInc{\factorize_exp}}%
       \let\factorize_N\factorize_Q 
       \factorize_iii
      }%
      {\factorize_iv}% 
}

\def\factorize_iv{%
    \xintiiifZero{\factorize_exp}%
      {}%
      {\edef\factors{\factors{{\factorize_M}{\factorize_exp}{\factorize_N}}}}%
    \xintiiifOne{\factorize_N}%
      {}%
      {% here N>1, N=QM+R (0<R<Q) is < M(Q+1) and N has no prime factors
       % at most equal to M. If a prime P>M divides N, the
       % quotient N/P will be < Q+1, hence at most Q. If Q<=M, then
       % N/P must be 1 else there would be some prime <=M dividing N.
       % no \xintiiifGeq ...
       % \xintiifCmp will have branches for each of <, =, >, less convenient
       % So we use \xintiiifLt which exists, and permute the branches
       % compared to original code
       \xintiiifLt\factorize_M\factorize_Q
         {% we go on testing with bigger factors
          % or \edef\factorize_M{\xintInc{\xintInc{\factorize_M}}} perhaps
          \edef\factorize_M{\xintiiAdd \factorize_M 2}%
          \def\factorize_exp{0}%
          \factorize_iii
         }%
         {% implies that N is prime
          \edef\factors{\factors{{\factorize_N}{1}{1}}}% we stop here
         }%
      }%
}

\catcode`\_ 8

%% We now define the macro \FactorTree which will produce a tree
%% displaying the factorization, 
%% using TikZ+forest (for the bracket syntax, much easier to deal with compared
%% to the braces-based child-node native TikZ tree syntax)


\makeatletter

\newtoks\FactorTreeA
\newtoks\FactorTreeB

\newcommand*{\FactorsToTree@}[3]{%
    \ifnum #2=1 % 
       \expandafter\@firstoftwo
    \else
       \expandafter\@secondoftwo
    \fi
    % exponent #2 is 1, so don't print it
    {\xintSgnFork{\xintCmp{#3}{1}}% check to see if end has been reached
    {}%
        % here, this is the last triplet and it has the shape {P}{1}{1}
        % and P was already inserted as tree node in the previous step
        % and the forest syntax allows to insert options here
    {\FactorTreeA\expandafter{\the\FactorTreeA ,draw,circle}}%
    {\FactorTreeA\expandafter{\the\FactorTreeA 
                             [{${#1}$}]
                             [{${#3}$}%
                             }%
     \FactorTreeB\expandafter{\the\FactorTreeB ]}%
    }}%
    % exponent #2 is > 1
    {\xintSgnFork{\xintCmp{#3}{1}}% check to see if end has been reached
    {}%
    {\FactorTreeA\expandafter{\the\FactorTreeA
                              [{${#1}^{#2}$}]
                              %[1]%
                              }%
    }%
    {\FactorTreeA\expandafter{\the\FactorTreeA 
                             [{${#1}^{#2}$}]
                             [{$#3$}%
                             }%
     \FactorTreeB\expandafter{\the\FactorTreeB ]}%
    }}%
}

\newcommand*{\FactorsToTree}[1]{\FactorsToTree@ #1}

% for factors displayed inline 

\def\@factorinliner  #1{\@factorinliner@ #1}
\def\@factorinliner@ #1#2#3{\ifnum #2>1 \expandafter\@firstoftwo\else
                                        \expandafter\@secondoftwo\fi
                            {{#1}^{#2}}{#1}}
\def\FactorsInline{%
    \xintListWithSep\times
             {\xintApply\@factorinliner{\expandafter\@gobble\factors}}%
}

\newlength{\horizontalshift}  % for positioning of the edges from the tree top
\newsavebox{\NandFactors}

\newcommand*{\FactorTree}[1]{%
    \factorize{#1}%
    \sbox{\NandFactors}{$#1=\FactorsInline$}%
    \setlength{\horizontalshift}{(\wd\NandFactors-\widthof{$#1$})/2}%
    \FactorTreeA{}%
    \FactorTreeB{}%
    \bracketset{action character=@}%
    \xintApplyUnbraced\FactorsToTree{\expandafter\@gobble\factors}%
    \begin{forest}
      for tree={edge path={\noexpand\path [\forestoption{edge}]
                                (!u.south)--(.child anchor);}},
      where={level()==1}
        {edge path= {\noexpand\path [\forestoption{edge}]
            ($(!u.south)-(\the\horizontalshift,0cm)$)--(.child anchor);}}{},
      where={n()==1}{draw,circle}{},
      [{\box\NandFactors}, right=\horizontalshift,
      @\the\FactorTreeA
      @\the\FactorTreeB
      ]
    \end{forest}
}

\makeatother



\begin{document}

\FactorTree{13}

\FactorTree{36}

\FactorTree{90}

\FactorTree{112}

\FactorTree{612}

\FactorTree{7875}

\FactorTree{22230}

\FactorTree{1073741824}

\FactorTree{2147483644}

\FactorTree{\xintiiPow 2{40}}

\FactorTree{\xinttheiiexpr 2^37 * 3^23 * 17^13\relax}%

% two examples ending with a (somewhat) `big' prime,

\FactorTree{10968058712}

\FactorTree{1689242184972}

\end{document}

This is half the answer, it just factorizes numbers. TeX is Turing-complete, no extra packages are required. :-) I'm not really in the mood for pst-tree, so I'll skip the second half. Besides, @Qrrbrbirlbel's answer (which I upvoted) covers that handsomely. If, for some reason, one prefers my code, it should be easy to patch the generic factor processing macros that I have with whatever is needed to generate the trees.

The algorithm I used is a quite simple one: it checks 2 and 3 and then checks all numbers that are 6k-1 or 6k+1, up to the square root of n. It is supposed to work with arbitrarily long numbers (up to 2^31, TeX's limit for counters), provided it does not exhaust TeX's memory when recursing. This happens, e.g., when you give it 2^31-1, which happens to be a Mersenne prime.

\documentclass[12pt]{article}

\makeatletter

\newcount\factorize@n  % the number to be factorized
\newcount\factorize@t  % temporary
\newcount\factorize@p  % a candidate factor
\newcount\factorize@c  % counter of factors

% machinery for factorization
%
\def\factorize#1{%
  \factorize@begin{#1}%
  \factorize@n=#1
  \factorize@c=0
  \factorize@p=2
  \factorize@try%
  \factorize@p=3
  \factorize@try%
  \factorize@p=5
  \factorize@loop%
  \ifnum\factorize@c>0%
    \ifnum\factorize@n>1%
      \factorize@process{\the\factorize@n}%
    \fi%
  \else%
    \factorize@process{\the\factorize@n}%
  \fi%
  \factorize@end{#1}%
}
\def\factorize@loop{%
  \factorize@t=\factorize@p
  \multiply\factorize@t by \factorize@p
  \ifnum\factorize@t>\factorize@n\else%
    \factorize@try%
    \advance\factorize@p by 2
    \factorize@t=\factorize@p
    \multiply\factorize@t by \factorize@p
    \ifnum\factorize@t>\factorize@n\else%
      \factorize@try%
      \advance\factorize@p by 4
      \factorize@loop%
    \fi%
  \fi%
}

\def\factorize@try{%
  \factorize@t=\factorize@n
  \divide\factorize@t by \factorize@p
  \multiply\factorize@t by \factorize@p
  \ifnum\factorize@n=\factorize@t%
    \factorize@process{\the\factorize@p}%
    \divide\factorize@n by \factorize@p
    \advance\factorize@c by 1
    \factorize@try%
  \fi%
}

% Stubs to be called at start, end, and when a factor is found
%
\def\factorize@begin#1{%
  \noindent%
  $#1 =%$
}
\def\factorize@end#1{%
  $%$
  \par%
}
\def\factorize@process#1{%
  \ifnum\factorize@c>0%
    \times%
  \fi%
  #1%
}

\makeatother

% Demo
%
\begin{document}
  \factorize{42}
  \factorize{5040}
  \factorize{1234567890}
  %\factorize{2147483647} exhausts TeX's memory
\end{document}

A small variation groups factors of multiplicity greater than one.

\documentclass[12pt]{article}

\makeatletter

\newcount\factorize@n  % the number to be factorized
\newcount\factorize@t  % temporary
\newcount\factorize@p  % a candidate factor
\newcount\factorize@c  % counter of factors (total)
\newcount\factorize@w  % counter of factors (power of given candidate)

% machinery for factorization
%
\def\factorize#1{%
  \factorize@begin{#1}%
  \factorize@n=#1
  \factorize@c=0
  \factorize@p=2
  \factorize@try%
  \factorize@p=3
  \factorize@try%
  \factorize@p=5
  \factorize@loop%
  \ifnum\factorize@c>0
    \ifnum\factorize@n>1
      \factorize@process{\the\factorize@n}{1}%
    \fi%
  \else%
    \factorize@process{\the\factorize@n}{1}%
  \fi%
  \factorize@end{#1}%
}
\def\factorize@loop{%
  \factorize@t=\factorize@p
  \multiply\factorize@t by \factorize@p
  \ifnum\factorize@t>\factorize@n\else%
    \factorize@try%
    \advance\factorize@p by 2
    \factorize@t=\factorize@p
    \multiply\factorize@t by \factorize@p
    \ifnum\factorize@t>\factorize@n\else%
      \factorize@try%
      \advance\factorize@p by 4
      \factorize@loop%
    \fi%
  \fi%
}
\def\factorize@try{%
  \factorize@w=0
  \factorize@try@aux%
  \ifnum\factorize@w>0
    \factorize@process{\the\factorize@p}{\the\factorize@w}%
    \advance\factorize@c by \factorize@w
  \fi%
}
\def\factorize@try@aux{%
  \factorize@t=\factorize@n
  \divide\factorize@t by \factorize@p
  \multiply\factorize@t by \factorize@p
  \ifnum\factorize@n=\factorize@t
    \divide\factorize@n by \factorize@p
    \advance\factorize@w by 1
    \factorize@try@aux%
  \fi%
}

% Stubs to be called at start, end, and when a factor is found
%
\def\factorize@begin#1{%
  \noindent%
  $#1 =%$
}
\def\factorize@end#1{%
  $%$
  \par%
}
\def\factorize@process#1#2{%
  \ifnum\factorize@c>0
    \times%
  \fi%
  \ifnum#2>1
    #1^{#2}%
  \else%
    #1%
  \fi%
}

\makeatother

% Demo
%
\begin{document}
  \factorize{42}
  \factorize{5040}
  \factorize{1234567890}
  %\factorize{2147483647} exhausts TeX's memory
\end{document}

MWE with Asymptote solution (not optimized). struct PrimeTree that draws a tree is defined in the preamble and is used in asy pictures. It consists of a function fillPrimeList(), which fills the list of prime numbers, a function factors(), which collects prime factors of the number, and a function draw(), which draws a tree.

% pfactors.tex :
\documentclass{article}
\usepackage{subcaption}
\usepackage{lmodern}
\usepackage[inline]{asymptote}
\usepackage[left=2cm,right=2cm]{geometry}
\begin{asydef}
struct PrimeTree{
  int num;
  int maxPrimeInd;
  int[] primelist;
  pen defpen,linkPen;

  void fillPrimeList(){
    int numd=num;
    int nmax;
    bool primetest;
    if(numd%2==0)numd=(int)(numd/2);
    if(numd%3==0)numd=(int)(numd/3);
    nmax=(int)(floor(sqrt(numd)));
    primelist=new int[nmax];
    primelist[0]=2;
    primelist[1]=3;

    maxPrimeInd=2;

    for(int i=5;i<=nmax;i+=2){
      primetest=true;
      for(int j=0;primetest && (j<maxPrimeInd);++j){
        if((i%primelist[j]==0)){
          primetest=false;
        }
      }
      if(primetest){
        if(numd%i==0){
          numd=(int)(numd/i);
          nmax=(int)(floor(sqrt(numd)));
        }
        primelist[maxPrimeInd]=i;
        ++maxPrimeInd;
      }    
    }
  }

  int[] factors(){
    fillPrimeList();
    int[] fs;
    int a,b;
    a=num;
    for(int i=0;i<maxPrimeInd;++i){
      b=primelist[i];
      while(a%b==0){
        fs.push(b);
        a=(int)(a/b);
      }
    }
    if(a>1)fs.push(a);
    return fs;
  }

  void draw(){  
    int[] fs=factors();
    int a,b;

    real row=0, col=0;
    real dr=1.618, dc=1.382;

    a=num; 
    for(int i=0; a>1 && i<fs.length;++i){
      b=fs[i];
      if(a!=b){
        label(string(a),(0,row),W);
        draw(Label(string(b),(col+dc,row),E),roundbox);
        draw((0,row)--(col+dc,row),linkPen);
        if(i>0) draw((0,row)--(0,row+dr),linkPen);
      }else{
        draw(Label(string(b),(col+dc,row),E),roundbox);  
        if(i>0) draw((0,row+dr)--(col+dc,row),linkPen);
      }
      row-=dr;
      a=(int)(a/b);
    }
  }
  void operator init(int num, pen defpen=darkblue+fontsize(10pt), pen linkPen=orange){
      assert(num>0);
      this.num=num;
      this.defpen=defpen;
      this.linkPen=linkPen; 
      defaultpen(defpen);
      draw();
  }
}

\end{asydef}

\begin{document}
\begin{figure}
\captionsetup[subfigure]{justification=centering}
\centering
\begin{subfigure}{0.24\textwidth}
\begin{asy}
unitsize(10pt);
PrimeTree(22230);
\end{asy}
\caption{}
\label{fig:1a}
\end{subfigure}
%
\begin{subfigure}{0.24\textwidth}
\begin{asy}
unitsize(10pt);
PrimeTree(1073741824);
\end{asy}
\caption{}
\label{fig:1b}
\end{subfigure}
%
\begin{subfigure}{0.24\textwidth}
\begin{asy}
unitsize(10pt);
PrimeTree(2147483644);
\end{asy}
\caption{}
\label{fig:1c}
\end{subfigure}
%
\begin{subfigure}{0.24\textwidth}
\begin{asy}
unitsize(10pt);
PrimeTree(19928968819);
\end{asy}
\caption{}
\label{fig:1d}
\end{subfigure}
\end{figure}
\end{document}    
%
% To process it with `latexmk`, create file `latexmkrc`:
% 
%     sub asy {return system("asy '$_[0]'");}
%     add_cus_dep("asy","eps",0,"asy");
%     add_cus_dep("asy","pdf",0,"asy");
%     add_cus_dep("asy","tex",0,"asy");
% 
% and run `latexmk -pdf pfactors.tex`.

A slightly different approach is to combine the typesetting of LaTeX with the abilities of other programming languages. Compared to the other solutions the code executes almost instantaneously, is relatively readable to a layman and much shorter.

Here python was used based on the great pythontex package and the forest package to draw the trees. Several sophisticated packages for number theory are additionally available for python, so a much more advanced algorithm can also be used. The code is a little messy, as I squished the creation of the string for the forest package into the factorization algorithm.

\documentclass{article}

\usepackage{tikz}
\usepackage{forest}
\usepackage{pythontex}


\begin{pycode}
from math import sqrt

def generate_tree(number):
        global tree_text
        tree_text = ""
        def prime_factors_recursive(n, level): 
                """Finds the prime factors of 'n' and generates text representation
                     according to tikz forest package.
                """ 
                limit = int(sqrt(n)) + 1
                divisor = 2
                num = n
                level += 1
                global tree_text
                tree_text = tree_text + "["
                for divisor in range(2, limit): 
                        if (num % divisor == 0): 
                                num /= divisor 
                                tree_text = tree_text + "%d [%d,circle,draw] " %(n, divisor)
                                return (n, (divisor, prime_factors_recursive(num, level)))
                tree_text = tree_text + "%d,circle,draw" %(n)
                for i in range(level):
                        tree_text = tree_text + "]"
                return (n,)

        prime_factors_recursive(number, 0)
        output = r'\begin{forest}' + tree_text + r'\end{forest}'
        return output

\end{pycode}

\newcommand{\PrimeTree}[1]{%
\py{generate_tree(#1)}
}

\begin{document}
\PrimeTree{36}
\PrimeTree{90}
\PrimeTree{112}
\PrimeTree{612}
\PrimeTree{7875}
\PrimeTree{22230}
\PrimeTree{19928968819}

\end{document}

I only discovered this topic today. So here is a very late attempt with MetaPost ^[1], probably very clumsy, but it seems to work nicely. I've not read the other contributions carefully yet, I may come back to catch some nice ideas and include them :-) Any remark is welcome.

Basically, the main program prime_tree calls a recursive macro, prime_factorization on the integer n, which first deals with 2 as divisor and then with the others (odd) prime divisors by calling another macro, odd_ prime_factorization.

Edit The code has been much simplified. In particular, no more odd_prime_factorization macro: prime_factorization now takes care of the main job alone. More improvements are probably to come.

Edit 2 The code has been thoroughly revised and completed in order to make any overlapping between the different nodes impossible. Also, it has been included in a LuaLaTeX program, which seems oddly enough a little quicker than standalone MetaPost.

\documentclass[12pt, border=5mm]{standalone}
\usepackage{luatex85, luamplib}
    \mplibtextextlabel{enable}
    \mplibnumbersystem{double}
\begin{document}
\begin{mplibcode}

vardef xradius(expr pic) = %horizontal radius of a picture
    save bbox_pic; path bbox_pic; bbox_pic = bbox pic;
    xpart(urcorner bbox_pic - center bbox_pic)
enddef;

vardef yradius(expr pic) = % vertical radius of a picture
    save bbox_pic; path bbox_pic; bbox_pic = bbox pic;
    ypart(urcorner bbox_pic - center bbox_pic)
enddef;

vardef radius(expr pic) = % circling a picture
    save bbox_pic; path bbox_pic; bbox_pic = bbox pic;
    max(xradius(bbox_pic), yradius(bbox_pic))
enddef;

vardef clearedlabel(expr str, pos) =    
    save newbox; path newbox; newbox = bbox thelabel(str, pos);     
    unfill newbox; label(str, pos);    
enddef;   

vardef circledlabel(expr str, pos) =   
    save circle; path circle;           
    circle = fullcircle scaled (2radius(textext(str))) shifted pos;       
    unfill circle; draw circle; label(str, pos);    
enddef;     

def prime_factorization(expr p, n, pos) =   
    if p**2 > n: % end test
        if pos = origin: circledlabel(decimal n, pos); % First n given was prime number! 
        else: % avoiding vertical overlapping with box above
            pair endpos; endpos = pos - (0, voffset - v + max(r1,r2));
            draw pos--endpos;
            circledlabel(decimal n, endpos);
        fi  
    % Main recursion step: branches below n
    elseif n mod p = 0: 
        picture pic_divisor, pic_dividend; numeric r[], dividend, hoffset, voffset;
        dividend = n div p;
        pic_divisor = textext(decimal p); pic_dividend = textext(decimal dividend);
        r1 = radius(pic_divisor); r2 = xradius(pic_dividend);  
        % Avoiding horizontal and vertical overlapping between circle and boxes
        hoffset = .5(h+r1+r2);
        voffset = v - r1;
        pair pos_dividend, pos_divisor; 
        pos_dividend =  pos + (hoffset, voffset); pos_divisor = pos + (-hoffset, voffset);
        draw pos -- pos_dividend; draw pos -- pos_divisor;
        clearedlabel(decimal n, pos); circledlabel(decimal p, pos_divisor);    
        prime_factorization(p, dividend, pos_dividend);
    elseif p = 2: prime_factorization (3, n, pos);      
    else: prime_factorization(p+2, n, pos); fi 
enddef;    

def prime_tree(expr n, pos) =
    numeric v, h; h = 0.75cm; v = -.65cm;
    prime_factorization(2, n, pos);
enddef;    

beginfig(1);
    L := 4cm; 
    prime_tree(36, origin);
    prime_tree(90, (L, 0));
    prime_tree(112, (2L, 0));
    prime_tree(612, (3L, 0));
    L := 5cm; l := -5cm;
    prime_tree(7875, (0, l));
    prime_tree(22230, (L, l)); 
    prime_tree(765434567654565, (2L, l));    
endfig;

\end{mplibcode}
\end{document}

Result:

The factorization of one of the examples provided by jfbu above took much longer time than the rest: 19928968819. The result is however as expected, after more than a minute of computations on my old laptop (MacBook Pro from 2008).

The next step is (I hope) the introduction of the exponents…

Edit 3 Done at last: a version taking the exponents into account. It has needed some time, the necessary changes having been difficult for me to sort out, mostly the appropriate way to deal with the last factors.

Like the former one, it has been incorporated in a LuaLaTeX program since it seems to run (slightly) faster than standalone MetaPost.

\documentclass[12pt, border=5mm]{standalone}
\usepackage{luatex85, luamplib}
    \mplibtextextlabel{enable}
    \mplibnumbersystem{double}
\begin{document}
\begin{mplibcode}

def str_number expr n = "$" & decimal n & "$" enddef;

def str_prime_exponent(expr p, cnt) =
    if cnt = 1: "$" & decimal p & "$"
    else: "$" & decimal p & "^{" & decimal cnt & "}$" fi
enddef;

vardef xradius(expr pic) = %horizontal radius of a picture
    save bbox_pic; path bbox_pic; bbox_pic = bbox pic;
    xpart(urcorner bbox_pic - center bbox_pic)
enddef;

vardef yradius(expr pic) = % vertical radius of a picture
    save bbox_pic; path bbox_pic; bbox_pic = bbox pic;
    ypart(urcorner bbox_pic - center bbox_pic)
enddef;

vardef radius(expr pic) = % circling a picture
    save bbox_pic; path bbox_pic; bbox_pic = bbox pic;
    max(xradius(bbox_pic), yradius(bbox_pic))
enddef;

vardef clearedlabel(expr str, pos) =    
    save newbox; path newbox; newbox = bbox thelabel(str, pos);     
    unfill newbox; label(str, pos);    
enddef;   

vardef circledlabel(expr str, pos) =   
    save circle; path circle;           
    circle = fullcircle scaled (2radius(textext(str))) shifted pos;       
    unfill circle; draw circle; label(str, pos);    
enddef; 

def show_divisor(expr p, cnt) =
    pair pos_divisor;
    pos_divisor = pos + (-hoffset, voffset);
    draw pos -- pos_divisor;
    circledlabel(str_prime_exponent(p, cnt), pos_divisor);
enddef;

def prepare_place_dividend_of(expr N) =
    pair pos_dividend;
    pos_dividend =  pos + (hoffset, voffset); 
    draw pos -- pos_dividend; 
    clearedlabel(str_number N, pos); 
enddef;   

def prime_factorization(expr p, n) =

    if p**2 <= n:

        if n mod p = 0:
            cnt := incr cnt;
            prime_factorization(p, n div p);
        else:
            if cnt > 0:
                picture pic_divisor, pic_dividend; numeric r[], hoffset, voffset;
                pic_divisor = textext(str_prime_exponent(p, cnt)); 
                pic_dividend = textext(str_number n);
                r1 = radius(pic_divisor); r2 = xradius(pic_dividend);
                % Avoiding horizontal and vertical overlapping between circle and boxes
                hoffset = .5(h+r1+r2); 
                voffset = v - r1;
                show_divisor(p, cnt);
                prepare_place_dividend_of(N);
                cnt := 0; N := n; pos := pos_dividend; 
            fi
            prime_factorization(if p = 2: 3 else: p+2 fi, n); 
        fi

    elseif cnt > 0:
        if p = n:
                cnt := incr cnt;
                picture pic_divisor; numeric r, hoffset, voffset;
                pic_divisor = textext(str_prime_exponent(p, cnt)); 
                r = radius(pic_divisor);
                hoffset = .5(h+r); 
                voffset = v - r;
                show_divisor(p, cnt);
                clearedlabel(str_number N, pos);            
        else:
                picture pic_divisor, pic_dividend; numeric r[], hoffset, voffset;
                pic_divisor = textext(str_prime_exponent(p, cnt)); 
                pic_dividend = textext(str_number(n));
                r1 = radius(pic_divisor); r2 = radius(pic_dividend);
                % Avoiding horizontal and vertical overlapping between circles and boxes
                hoffset = .5(h+r1+r2); 
                voffset = v - max(r1, r2);
                show_divisor(p, cnt);
                prepare_place_dividend_of(N);               
                circledlabel(str_number n, pos_dividend);
        fi

    elseif pos = origin: 
            circledlabel("$" & decimal n & "$", pos);

    else: % avoiding vertical overlapping with box above
          pair endpos; 
          numeric r; r = radius(textext(str_number(n)));
          endpos = pos - (0, voffset - v + max(r1,r));
          draw pos--endpos;
          circledlabel(str_number n, endpos); 
    fi
enddef;

def prime_tree(expr n) =
    numeric v, h, cnt, N; 
    h = 0.75cm; v = -.65cm; cnt = 0; N = n;
    prime_factorization(2, n);
enddef;    

beginfig(1);
    numeric L; L = 4cm; numeric l; l := -5cm;
    pair pos;
    pos := origin; prime_tree(36);
    pos := (L, 0); prime_tree(90);
    pos := (2L, 0); prime_tree(112);
    pos := (3L, 0); prime_tree(612); 
    L := 4.5cm;
    pos := (0, l); prime_tree(7875);
    pos := (L, l); prime_tree(22230); 
    pos := (2L, l); prime_tree(765434567654565);  
endfig;

\end{mplibcode}
\end{document}

Output:

This is an answer to a topic raised by the OP in a comment, not an answer to the initial question of trees (this is why I make a second answer).

PSTikZ ^[1] asked for macros computing the greatest common divisor and the least common multiple of a set of integers.

The initial version of this answer provided such macros based on \xintGCDand \xintLCMfrom package xintgcd ^[2].

Since 1.09a (2013/09/24) the package actually provides \xintGCDofand \xintLCMof and xintexpr ^[3] has functions gcd and lcm. Thus the answer here got (belatedly) updated to simply use them.

\documentclass{article}
\usepackage{amsmath}
\usepackage{xintgcd}% \xintGCDof and \xintLCMof of macros
\usepackage{xintexpr}% gcd and lcm functions

\newcommand\test[1]{$\texttt{\detokenize{\xintGCDof{#1}}}\to\xintGCDof{#1}$\par
  $\texttt{\detokenize{\xintLCMof{#1}}}\to\xintLCMof{#1}$\par}

\begin{document}

\test {{12}}

\test {{12}{8}}

\test {{144}{216}{360}}

\test {{125}{625}{1000}}

\begin{align*}
  \mathrm{GCD}
  \begin{pmatrix}
    64*9*125*49*121*13\\ 16*81*25*343*1331*169\\ 128*27*5*7*121*169\\
    8*27*25*49*121*13
  \end{pmatrix}
  &=
  \xinttheiiexpr gcd(64*9*125*49*121*13,
  16*81*25*343*1331*169, 128*27*5*7*121*169, 8*27*25*49*121*13)\relax\\
8*9*5*7*121*13
 &= \xinttheiiexpr 8*9*5*7*121*13\relax\\
\mathrm{LCM}
  \begin{pmatrix}
    64*9*125*49*121*13\\ 16*81*25*343*1331*169\\ 128*27*5*7*121*169\\
    8*27*25*49*121*13
  \end{pmatrix}
  &=
  \xinttheiiexpr lcm(64*9*125*49*121*13,
  16*81*25*343*1331*169, 128*27*5*7*121*169, 8*27*25*49*121*13)\relax\\
128*81*125*343*1331*169
&= \xinttheiiexpr 128*81*125*343*1331*169\relax
\end{align*}

\end{document}