Simplifying square roots

Here's a LuaLaTeX-based solution. The code sets up a LaTeX macro named \rsqrt, which invokes a Lua function named rsqrt. The latter implements the simplification algorithm you've proposed -- with the following refinements:

• For n=0 or n=1, the code simply returns n (without a square root symbol); and

• Care is taken to omit the \sqrt{n/i²} term if it's equal to 1, i.e., if n is a "square number" (4, 9, 16, etc) or a product of smaller square numbers. E.g, if n=36, \rsqrt{36} shows 6 since 36 = 6^2 = 2^2*3^2.

No input sanity checking is performed, i.e., the user is responsible for providing an argument to \rsqrt that is either a non-negative whole number or evaluates to a non-negative whole number. Thus, it's ok to write \rsqrt{1e6} and \rsqrt{3.6e7}: the macro will return 1000 and 6000, respectively.

Observe that the macro \rsqrt must be used in math mode as it may output \sqrt directives.

% !TEX TS-program = lualatex  
%% Note: Code updated 2019/10/26 to work with LaTeX 2019-10-01 

%% Create an external file to contain the Lua code
\begin{filecontents*}[overwrite]{rsqrt.lua}

function rsqrt ( n )
  -- n : a non-negative whole number (or something
  --     that evaluates to a non-neg. whole number)
  if n == 0 or n == 1 then   -- Nothing to do
    return ( n )
  else
    i = math.floor ( math.sqrt ( n ) )
    while i > 1 do
      if ( n % i^2 == 0 ) then     -- n is divisible by i^2
        k = math.floor ( n / i^2 ) -- 'math.floor' makes k an explicit integer
        if k == 1 then  -- n is a "square" number (or a product of square numbers)
          return ( i )
        else
          return ( i .. "\\sqrt{" .. k .. "}" )
        end
      end
      i = i-1
    end
    -- No simplification possible:
    return ( "\\sqrt{" .. n .. "}" )
  end
end
-- Define a vector (in form of a Lua table) of whole numbers
nvec = {0,1,2,3,4,5,7,12,16,18,27,32}

-- Lua function to print 3-column array:
function PrintArray()
  for i=1,#nvec do
    u = nvec[i]
    tex.sprint ( math.floor(u) .. 
                "& \\sqrt{" .. math.floor(u) .. 
                "}&" .. rsqrt(u) .. "\\\\" )
  end
end

\end{filecontents*}

\documentclass{article}

%% Load Lua code from external file:
\directlua{dofile("rsqrt.lua")}

%% TeX-side code: "wrapper" macro that invokes the Lua function:
\newcommand\rsqrt[1]{\directlua{tex.sprint(rsqrt(#1))}}

\begin{document}
\[
\renewcommand\arraystretch{1.25}
\begin{array}{@{} rcc @{}}
  \verb+n+ & \verb+\sqrt+ & \verb+\rsqrt+ \\ % print header row
  \directlua{PrintArray()}    % create and print body of 'array'
\end{array}
\]
\end{document} 

The algorithm in the question is very inefficient: except of course if the original integer is a perfect square.

This answer (in chronological order):

1. an approach with macros which mimics the simplest factorization algorithm,

2. an expandable approach using the algorithm as in OP. update very embarrassingly, the author had not understood the OP's algorithm and after having found a simplifying I such as I^2 divided N, he went one recursively with N<-N/I^2 not understanding that the algorithm could stop there. (as a pale excuse, he had implemented first the "bottom-up" way, which needs recursivity, contrarily to the "top-down" (less efficient) way). The answer is thus updated, apologies to all generous trustful early up-voters.

   I am updating again (sorry) because I have now read more of xintexpr.sty's documentation, and I switched for efficiency from i=sqrt(N)..1 to i=-sqrt(N)++ (there is no -- available, hence a trick with minus sign). The former generates beforehand the whole list of floor(\sqrt{N}) numbers (sqrt means truncated square root in \xintiiexpr), the latter is an iterator which doesn't generate things beforehand. Furthermore the former syntax can only generate about 5000 values (sqrt(N)..[-1]..1 would have no such limitation, but still would generate everything beforehand).

3. an expandable implementation of faster (high-school factorization type) algorithm as in approach 1.

To be honest 2., 3. and even 1. would probably be better written entirely using \numexpr as their pretence at handling numbers bigger than 2^31 is a bit far-stretched, it takes time with a prime number of 10 digits to do tens of thousands of divisions to conclude it was square-free ... Implementation 2. has an intrinsic limit to 2^62 because the square root must be a TeX number (due to some construct inside).

In 2. and 3., we push a bit beyond reasonable range the possibilities of \xintexpr syntax with recursive sequences. The notation is a bit cumbersome. Besides xintexpr.sty 1.2g is needed because it changed relevant syntax.

4. finally (2017) I also add the no-package numexpr only expandable approach.

First approach (we change the algorithm)

Not to say that this is an easy problem. A bit of googling reveals that apparently mathematicians currently believe that finding the square free radical of an integer may be as hard as complete factorization: https://math.stackexchange.com/questions/171568/finding-the-radical-of-an-integer and https://math.stackexchange.com/questions/14667/square-free-integers-factorization.

Here is an approach (using macros) which mimicks simplest form of factorization algorithm.

Package xintexpr is used only in order to allow inputs such as 1e7 or even expressions. It also loads xinttools which is used in the syntax.

Apart from that all operations are done with macros available from xint. As we deal in the example almost only with numbers <2^31 we could employ a variant where all operations would be done using uniquely \numexpr, naturally that would be quite faster.

The code uses \xintiiDivision which computes simultaneously quotient and remainder. This is why \xintAssign is used to store them into two macros \A and \B. The code examines if \B vanishes to detect the divisibility by a Q=P^2.

\documentclass[a4paper]{article}

\usepackage{geometry}

\usepackage{xintexpr}

\makeatletter
\def\Rsqrt@ {%
    \let\Nrad\N
    \def\Nroot {1}% 
% we will always have original N = \Nrad times \Nroot^2
% first we check powers of 2
    \def\P{2}%
    \def\Q{4}% \Q is always square of \P
    \xintloop
% try to divide \Nrad by 4. If possible, multiply \Nroot by 2
        \xintAssign\xintiiDivision{\Nrad}{\Q}\to \A\B
        \xintiiifZero{\B}
           {\let\Nrad\A
            \edef\Nroot{\xintiiMul{\Nroot}{\P}}% 
            \iftrue}
           {\iffalse}%
    \repeat
% try to divide \Nrad by 9=3^2, then by 25=5^2, etc...
% unfortunately we divide by all odd integers, but only odd prime
% integers would be really needed
    \def\P{3}%
    \xintloop
      \edef\Q{\xintiiSqr{\P}}%
      \xintiiifGt{\Q}{\Nrad}
         {\iffalse}%
         {\xintloop
             \xintAssign\xintiiDivision{\Nrad}{\Q}\to \A\B
             \xintiiifZero{\B}
                {\let\Nrad\A
                 \edef\Nroot{\xintiiMul{\P}{\Nroot}}% 
                 \iftrue}
                {\iffalse}%
          \repeat 
          \edef\P{\xintiiAdd{2}{\P}}%
          \iftrue
         }%
    \repeat
% at this stage \N = \Nrad times \Nroot^2
% and \Nrad is square-free.
    \xintiiifOne{\Nroot}{}{\Nroot}%
    \xintiiifOne{\Nrad} {}{\sqrt{\Nrad}}%
}%

\newcommand* \Rsqrt[1]{%
   \begingroup
     \edef\N{\xinttheiexpr #1\relax}%
     \xintiiifSgn \N
        {\pm\edef\N{\xintiiAbs{\N}}\xintiiifOne\N{}{\Rsqrt@}i}
        {0}
        {\xintiiifOne \N{1}{\Rsqrt@}}
   \endgroup
}

\makeatother
\usepackage{multicol}

\begin{document}

\parindent0pt\def\columnseprule{.4pt}%

% testing
% \begin{multicols}{4}
% \xintFor* #1 in {\xintSeq {10000}{10100}}\do
%        {$\sqrt{#1}=\Rsqrt{#1}$\par}
% \end{multicols}

% $\Rsqrt{-10}, \Rsqrt{-1}, \Rsqrt{-16}$

$\Rsqrt {1e6}, \Rsqrt {1e7}, \Rsqrt{1e21}$

\pdfsetrandomseed 123456789

\begin{multicols}{3}
\xintiloop [1+1]
    \edef\Z {\xinttheiiexpr 
             (\pdfuniformdeviate 1000)^2
             *\pdfuniformdeviate 1000\relax }%
    $\sqrt{\Z}=\Rsqrt{\Z}$\par
\ifnum\xintiloopindex<50
\repeat
\end{multicols}
    
\end{document}

---

Second approach (same algorithm as in OP, expandable implementation)

With the original algorithm. Here we define \ExtractRadical which expandably returns A,B with N=A^2 B. A non-expandable wrapper re-cycles \Rsqrt from faster approach above to produce A\sqrt{B}, distinguishing cases of negative N or N=0, 1.

I have added code comments to explain implementation. An earlier version was very lame (see top of answer) and additionally required xintexpr 1.2g which is not the case anymore.

\documentclass[a4paper]{article}

\usepackage{geometry}

\usepackage{xintexpr}

% Aim: given N find largest I such as I^2 divides N.
% Then compute J=N/I^2 and print I\sqrt{J}.

% Algorithm: compute first truncated square root Imax of N.
% With I = Imax try to divide N by I^2, if it does not work
% repeat with I replaced by I-1 and so on. 
% As soon as it works the seeked-for I has been found.

% **** Notice: embarrassingly the author of this answer initially continued
% **** the algorithm recursively with N<-N/I^2, which was very stupid, but
% **** explainable from the fact that he had implemented first another (much
% **** faster) algorithm which divided not from the top down, but from the
% **** bottom up.

% The code is far simpler now. And it does not require xintexpr 1.2g anymore,
% earlier versions of xintexpr.sty work, too.

% The iteration over i used Imax..1 syntax which requires Imax 
% to be <2^31. Else we could use Imax..[-1]..1, but we don't
% really consider realistic to iterate over 2^31 or more values !

% After an update we use (-Imax)++ syntax; this also requires Imax<2^31.

\def\ExtractRadical #1{%
  \xinttheiiexpr 
  subs(
  % we return I, #1//I^2 where I is biggest integer such as I^2 divides #1.
     (I, #1//I^2), 
  % The I is computed via the "seq" here. Theoretically this "seq" 
  % evaluates as many values as the last list indicates.
  % But here we omit all i's such that i^2 does not divide #1
  % and as soon as we have found one, we stop here and now by 
  % "break". We work topdown, at the worst I=1.
% The i=A..B syntax pre-generates all values, which is wasteful
% and limited to about at most 5000 values.
%       I=seq((#1/:i^2)?{omit}{break(i)}, i=sqrt(#1)..1)
% On the contrary the N++ syntax does not pre-generate anything.
       I=seq((#1/:i^2)?{omit}{break(-i)}, i=-sqrt(#1)++)
% There is currently no "n--" only "n++", thus we tricked with a minus sign.
      )
  \relax
}

\makeatletter

\def\Rsqrt@ {\expandafter\Rsqrt@@\romannumeral-`0\ExtractRadical\N,}

% The #2#3 trick is to get rid of a space after the comma
% because \ExtractRadical does \xinttheiiexpr which in case
% of comma separated values on output always inserts such a space.
% Naturally as the typesetting is in math mode the space is
% not a real problem (it is not a problem either in \xintiiifOne 
% as here its argument is already expanded anyhow).
\def\Rsqrt@@ #1,#2#3,{\xintiiifOne{#1}{}{#1}\xintiiifOne{#2#3}{}{\sqrt{#2#3}}}

\newcommand* \Rsqrt[1]{%
   \begingroup
     \edef\N{\xinttheiexpr #1\relax}%
     \xintiiifSgn \N
        {\pm\edef\N{\xintiiAbs{\N}}\xintiiifOne\N{}{\Rsqrt@}i}
        {0}
        {\xintiiifOne \N{1}{\Rsqrt@}}
   \endgroup
 }

\makeatother
\usepackage{multicol}

\begin{document}

\parindent0pt\def\columnseprule{.4pt}%

% testing

\begin{multicols}{4}
\xintFor* #1 in {\xintSeq {10000}{10099}}\do
       {$\sqrt{#1}=\Rsqrt{#1}$\par}
\end{multicols}

% $\Rsqrt{-10}, \Rsqrt{-1}, \Rsqrt{-16}$

$\Rsqrt {1e6}, \Rsqrt {1e7}$%, 

% this one does not work because 10^10.5 > 2^31 causes an arithmetic
% overflow in the "sqrt(J)..1" part. 
% It would not overflow with "sqrt(J)..[-1]..1"
% but then we can wait long time ... 
% from 31622776601 downto
%      10000000000 that's a lot of iterations !
%$\Rsqrt{1e21}$
% The update uses n++ syntax, but this also requires abs(n) to be <2^31
% hence the same remark applies: a "Number too big" error is generated.
% Better actually than to wait the completion of 21622776601 iterations ;-)


% \stop

\pdfsetrandomseed 123456789

% we try with smaller numbers... 1000 replaced by 100...
\begin{multicols}{3}
\xintiloop [1+1]
    \edef\Z {\xinttheiiexpr 
             (\pdfuniformdeviate 100)^2
             *\pdfuniformdeviate 100\relax }%
    $\sqrt{\Z}=\Rsqrt{\Z}$\par
\ifnum\xintiloopindex<51
\repeat
\end{multicols}

\end{document}

---

Third approach: again the faster algorithm, but expandably.

It would be more reasonable to code this à la \numexpr but well. Details in the code comment. The example has now 51 random examples, and funny enough the missing one (from first approach) turned out to be a random square (with the random seed to pdftex random number generator set at 123456789).

\documentclass[a4paper]{article}

\usepackage{geometry}

\usepackage{xintexpr}[2016/03/19]%
% needs xintexpr 1.2g due to 
%  - changed meaning of iter
%  - shift by 1 in [L][n] syntax


% syntax \ExtractRadical {N or <integer expression>} expands to "A, B" with
% N=A^2 B, B square-free
% Algorithm:
% main variable a triple (P, I, J) where
%  - always N = I^2 J
%  - J's prime factors < P have multiplicity one.
% START: (2, 1, N)
% ITER:  (P, I, J)
         % Q=P^2
         % is Q > J ?
         %   - yes: STOP(I, J)
         %   - no:
         %         does Q divide J ?
         %         - yes: I<-I*P, J<-J/Q and repeat until Q does not divide J
         %         - no; continue with (P+2, I, J). Except if P=2 then we go
         %                      on with P=3.
% Also works with N=0 (produces 1, 0) and with N=1 (produces 1, 1)
%

\newcommand*\ExtractRadical [1]{%
  \xinttheiiexpr 
  iter (2, 1, #1; % starting triple P=2, I=1, J=N
  subs(subs(subs(subs(
  % apart from Q=P^2, these substitutions are mainly because [@][n] syntax
  % is cumbersome; and inefficient as it allows n to be itself a complicated
  % expression, hence does some rather unneeded parsing here of n= 0, 1, 2.
  % We really need some better syntax like iter(P=2, I=1, J=#1;...) and then
  % work with P, I, J standing for the last values.
  % Or at least something like subs(..., (Q, P, I, J)=(...)).
  % (not yet with xintexpr 1.2g). 
  (Q>J)?
      {break(I, J)}
      {(J/:Q)?
          {(n)?{P+2}{3}, I, J}
      % must use parentheses here: ([@][1]). Else ]/: will confuse parser.
      % I could have used again subs, but well. 
          {iter(P*I,J//Q;(([@][1])/:Q)?{break((n)?{P+2}{3},@)}
                                       {(P*[@][0],([@][1])//Q)},e=1++)
          }
      }
  , Q=P^2), P=[@][0]), I=[@][1]), J=[@][2]), n=0++)
    \relax
}             
                  
\makeatletter
\def\Rsqrt@ {\expandafter\Rsqrt@@\romannumeral-`0\ExtractRadical\N,}

\def\Rsqrt@@ #1,#2,{\xintiiifOne{#1}{}{#1}\xintiiifOne{#2}{}{\sqrt{#2}}}

\newcommand* \Rsqrt[1]{%
   \begingroup
     \edef\N{\xinttheiexpr #1\relax}%
     \xintiiifSgn \N
        {\pm\edef\N{\xintiiAbs{\N}}\xintiiifOne\N{}{\Rsqrt@}i}
        {0}
        {\xintiiifOne \N{1}{\Rsqrt@}}
   \endgroup
 }

\makeatother
\usepackage{multicol}

\begin{document}

\parindent0pt\def\columnseprule{.4pt}%

% testing

% \xintFor* #1 in {\xintSeq {0}{50}}\do 
% {\ExtractRadical {#1}\par}

% \ExtractRadical {128}

% \ExtractRadical {1024}
% \stop

% $\Rsqrt{5000}$

% \stop

% \begin{multicols}{4}
% \xintFor* #1 in {\xintSeq {10000}{10099}}\do
%        {$\sqrt{#1}=\Rsqrt{#1}$\par}
% \end{multicols}

$\Rsqrt{-10}, \Rsqrt{-1}, \Rsqrt{-16}$

$\Rsqrt {1e6}, \Rsqrt {1e7}, \Rsqrt {1e21}$%, 

\pdfsetrandomseed 123456789

\begin{multicols}{3}
\xintiloop [1+1]
    \edef\Z {\xinttheiiexpr 
             (\pdfuniformdeviate 1000)^2
             *\pdfuniformdeviate 1000\relax }%
    $\sqrt{\Z}=\Rsqrt{\Z}$\par
\ifnum\xintiloopindex<51
\repeat
\end{multicols}

\end{document}

---

Update (2017).

Here is no-package expandable macro, if needed. It expands to I,J where original N is I**2 times J and J is square-free. Uses only \numexpr. Limited to <2**31 integers. Works from low to big, same as elementary factorization method. Stopping criteria should be improved, the comment below applies here too.

\makeatletter

\newcommand\ExtractRadical[1]{%
    \romannumeral0%
    \expandafter
    \ExtractRadical@two@i\expandafter1\expandafter,\the\numexpr#1.%
}%
\def\ExtractRadical@two@i #1,#2.{%
    \ifnum4>#2 \expandafter\ExtractRadical@two@done\fi
    \expandafter\ExtractRadical@two@ii\the\numexpr#2/4;#1,#2.%
}%
\edef\ExtractRadical@two@done #1;#2,#3.%
    {\space#2,#3}% (not sole #2 for readability)
\def\ExtractRadical@two@ii #1;#2,#3.{%
    \ifnum\numexpr#1*4=#3 
      \expandafter\@firstoftwo
    \else
      \expandafter\@secondoftwo
    \fi
   {\expandafter\ExtractRadical@two@i\the\numexpr2*#2,#1.}%
   {\ExtractRadical@i 3;#2,#3.}%
}%
\def\ExtractRadical@i #1;{%
    \expandafter\ExtractRadical@ii\the\numexpr#1*#1.#1;%
}%
\def\ExtractRadical@ii #1.#2;#3,#4.{%
    \ifnum#1>#4 \expandafter\ExtractRadical@done\fi
    \expandafter\ExtractRadical@iii\the\numexpr#4/#1.#1;#4.#2;#3.%
}%
\def\ExtractRadical@iii #1.#2;#3.{%
    \ifnum\numexpr#1*#2=#3
      \expandafter\@firstoftwo
    \else
      \expandafter\@secondoftwo
    \fi
    \ExtractRadical@update
    \ExtractRadical@next
    #1.#2;#3.%
}%
\def\ExtractRadical@update #1.#2;#3.#4;#5.{%
    \expandafter\ExtractRadical@ii
    \the\numexpr#2\expandafter.%
    \the\numexpr#4\expandafter;%
    \the\numexpr#4*#5,#1.%
}%
\def\ExtractRadical@next #1.#2;#3.#4;#5.{%
    \expandafter\ExtractRadical@i\the\numexpr2+#4;#5,#3.%
}%
\edef\ExtractRadical@done #1;#2.#3;#4.{\space#4,#2}%

\makeatother

In expl3:

\documentclass{article}

\usepackage{xparse}

\ExplSyntaxOn
\NewDocumentCommand{\rsqrt}{m}
 {
  \manual_rsqrt:n { #1 }
 }

\int_new:N \l_manual_rsqrt_int

\cs_new_protected:Nn \manual_rsqrt:n
 {
  \int_set:Nn \l_manual_rsqrt_int { \fp_to_decimal:n { trunc(sqrt(#1),0) } }
  \bool_until_do:nn
   {
    \int_compare_p:n { \int_mod:nn { #1 } { \l_manual_rsqrt_int * \l_manual_rsqrt_int } == 0 }
   }
   {
    \int_decr:N \l_manual_rsqrt_int
   }
  \int_compare:nTF { \l_manual_rsqrt_int == 1 }
   {
    \sqrt{#1}
   }
   {
    \int_to_arabic:n { \l_manual_rsqrt_int }
    \int_compare:nF { #1 == \l_manual_rsqrt_int*\l_manual_rsqrt_int }
     {
      \sqrt{ \int_to_arabic:n { #1/(\l_manual_rsqrt_int*\l_manual_rsqrt_int) } }
     }
   }
 }
\ExplSyntaxOff

\begin{document}

$\rsqrt{4}$

$\rsqrt{8}$

$\rsqrt{18}$

$\rsqrt{12}$

$\rsqrt{7}$

\end{document}

If you want to deal with arguments 0 and 1 and also with negative arguments, you can change the main definition into

\NewDocumentCommand{\rsqrt}{m}
 {
  \int_compare:nTF { #1 < 0 }
   {
    \int_compare:nTF { #1 = -1 } { i } { \manual_rsqrt:n { -#1 } i }
   }
   {
    \int_compare:nTF { #1 < 2 } { #1 } { \manual_rsqrt:n { #1 } }
   }
 }

Now the input

$\rsqrt{0}$ and $\rsqrt{1}$

$\rsqrt{-1}$ and $\rsqrt{-4}$ and $\rsqrt{-32}$

would output