20

Recently, I started to experimenting with expl3 programming language. I found it very powerful for solving problems not only related with typesetting, but also with problems the other programming languages can do.

For good practice, I chose the problems from Project Euler and with the help of tex.se community, I managed to solve some problems that I once solved using a Python.

Below I quote the solution code for one of the problem named Longest Collatz sequence. The code is not perfect, because it uses a lot of estimated time. Perhaps there are more elegant algorithms for solving this problem using expl3.

In this section, I would like to see other solutions to this problem, not only with exp3, but also with other possibilities that can be implemented using LaTeX (e.g. LuaLaTeX and so on)


\documentclass[12pt]{article}

\usepackage{xparse, pgf}
\ExplSyntaxOn

\int_new:N \l_max_int 
\int_zero:N \l_max_int
\int_new:N \l_max_num_int 
\int_zero:N \l_max_num_int
\int_new:N \l_a_int 
\int_zero:N \l_a_int
\int_set:Nn \l_tmpa_int {1}

\NewDocumentCommand{\MaxNumberB}{ mmm }
{
    \int_while_do:nn {\l_tmpa_int <= #3}
    {
    \int_set:Nn \l_a_int { \l_tmpa_int }
    \int_set:Nn \l_tmpb_int { 1 }
        \int_while_do:nn{\l_a_int != 1}
        {
            \int_if_even:nTF {\l_a_int}
                {\int_set:Nn \l_a_int {\l_a_int/2} } 
                {\int_set:Nn \l_a_int {3*\l_a_int + 1} }
                \int_incr:N \l_tmpb_int
        }
        \int_compare:nTF{\l_tmpb_int > \l_max_int}
        {
        \int_set:Nn \l_max_int {\l_tmpb_int}
        \int_set:Nn \l_max_num_int {\l_tmpa_int}
        }{}

    \int_incr:N \l_tmpa_int
    } 

       \cs_new:Npx #1 { \int_eval:n { \l_max_int } }
       \cs_new:Npx #2 { \int_eval:n { \l_max_num_int } }

}

\ExplSyntaxOff

\begin{document}

\pgfmathsetmacro{\num}{1000000}
\MaxNumberB{\MaxNumEl}{\MaxNum}{\num}

For numbers from $0$ to $\num$, the number \MaxNum{} has the maximum number of elements in Collatz sequence~--- \MaxNumEl.


\end{document}
2
  • 1
    i would still use python for this via pythontex. – percusse Jun 1 '17 at 22:01
  • @percusse Of course, I know how to solve a problem with Python. But still, if it's not difficult for you, please bring your pythontex decision, because it is also within the scope of the section – sergiokapone Jun 2 '17 at 10:46
21

Reimplementation of Joseph's answer in LuaTeX. The nice thing about doing computations in Lua is that they are always fully-expandable.

\documentclass{article}
\usepackage{luacode}

\begin{luacode}
function collatz_next(n)
    if ( n \% 2 == 0 ) then
        return .5 * n
    else
        return 3 * n + 1
    end
end

function collatz_count(n, m)
    m = m or 1

    if ( n > 1 ) then
        return collatz_count( collatz_next(n), m + 1 )
    else
        return m
    end
end

function collatz_max_range(start, stop, count, len)
    count = count or 1
    len = len or start
    if start > stop then
        return len
    else
        if collatz_count(start) > count then
            return collatz_max_range(start+1, stop, collatz_count(start), start)
        else
            return collatz_max_range(start+1, stop, count, len)
        end
    end
end
\end{luacode}

\newcommand*\collatzcount[1]{\directlua{tex.sprint(collatz_count(#1))}}
\newcommand*\collatzmaxrange[2]{\directlua{tex.sprint(collatz_max_range(#1,#2))}}

\begin{document}

\edef\x{\collatzcount{13}}
\show\x

\edef\x{\collatzmaxrange{1}{10}}
\show\x

\end{document}

Performance is reasonably good for \collatzmaxrange{1}{1000000} on an Intel(R) Core(TM) i7-4790K.

$ time lualatex --interaction=batchmode test.tex
This is LuaTeX, Version 1.0.4 (TeX Live 2017) 
 restricted system commands enabled.

luaotfload | main : initialization completed in 0.070 seconds
real    0m10.875s
user    0m10.860s
sys 0m0.012s

LuaJIT is quite a bit faster than Lua.

$ time luajittex -fmt=luajitlatex --interaction=batchmode test.tex
This is LuajitTeX, Version 1.0.4 (TeX Live 2017) 
 restricted system commands enabled.

luaotfload | main : initialization completed in 0.048 seconds
real    0m4.675s
user    0m4.644s
sys 0m0.008s

Let's actually turn on the JIT.

$ time luajittex -jiton -fmt=luajitlatex --interaction=batchmode test.tex
This is LuajitTeX, Version 1.0.4 (TeX Live 2017) 
 restricted system commands enabled.

luaotfload | main : initialization completed in 0.048 seconds
real    0m0.807s
user    0m0.780s
sys 0m0.024s

Okay, this is as fast as it can get.

\documentclass{article}

\directlua{collatz = require"collatz"}

\newcommand*\collatzcount[1]{\directlua{tex.sprint(collatz.count(#1))}}
\newcommand*\collatzmaxrange[2]{\directlua{tex.sprint(collatz.max_range(#1,#2))}}

\begin{document}

\collatzmaxrange{1}{1000000}

\end{document}
$ time lualatex -interaction=batchmode test.tex
This is LuaTeX, Version 1.0.4 (TeX Live 2017) 
 restricted system commands enabled.

luaotfload | main : initialization completed in 0.069 seconds
real    0m0.367s
user    0m0.360s
sys 0m0.004s

Compiled with clang -shared -fpic -O3 -I/usr/include/lua5.2 test.c -o collatz.so

#include <lua.h>
#include <lauxlib.h>
#include <assert.h>

typedef unsigned long int_t;

inline int_t collatz_next(int_t n)
{
    if ( n % 2 == 0 )
      return n / 2;
    else
      return 3 * n + 1;
}

int_t collatz_count(int_t n, int_t m)
{
    if ( n > 1 )
      return collatz_count( collatz_next(n), m + 1 );
    else
      return m;
}

int_t collatz_max_range(int_t start, int_t stop, int_t count, int_t len)
{
  if ( start > stop )
    return len;
  else
  {
    if ( collatz_count(start,1) > count )
      return collatz_max_range(start+1, stop, collatz_count(start,1), start);
    else
      return collatz_max_range(start+1, stop, count, len);
  }
}

static int lcollatz_count(lua_State *L)
{
  assert( lua_gettop(L) == 1  && "Number of required arguments is 1" );
  assert( lua_isnumber(L, -1) && "First argument is not a number"    );

  int_t n = lua_tonumber(L, -1);

  int_t result = collatz_count(n,1);
  lua_pushnumber(L, result);

  return 1;
}

static int lcollatz_max_range(lua_State *L)
{
  assert( lua_gettop(L) == 2  && "Number of required arguments is 2" );
  assert( lua_isnumber(L, -2) && "First argument is not a number"    );
  assert( lua_isnumber(L, -1) && "Second argument is not a number"   );

  int_t start = lua_tonumber(L, -2);
  int_t stop  = lua_tonumber(L, -1);

  int_t result = collatz_max_range(start,stop,1,start);
  lua_pushnumber(L, result);

  return 1;
}

static const luaL_Reg lib_collatz[] = {
  {"count", lcollatz_count},
  {"max_range", lcollatz_max_range},
  {NULL, NULL} // sentinel
};

LUALIB_API int luaopen_collatz(lua_State *L)
{
  luaL_newlib(L, lib_collatz);
  return 1;
}
13
  • @jfbu This seems to be a problem with the Lua interpreter itself. On my old laptop I get with LuaTeX 0.95.0 (Lua 5.2) real 0m29.149s, whereas with the bare Lua 5.3 interpreter I get real 0m37.102s. That is a slowdown of 20%! And it is a good example for how the language gets more heavyweight over time. – Henri Menke Jun 2 '17 at 8:24
  • @jfbu Did you always run it multiple times so the cache is hot? – Henri Menke Jun 2 '17 at 8:43
  • the problem appears to be one of the x86_64-darwinlegacy luatex binary. It is twice slower as the x86_64-darwin one. It is the one in (nearly out) TeXLive 2017 for macOSes up to 10.9.x, which is my case. Anyway, I will never ever make my tex files depend upon using an engine in the hands of a very few people behaving as other-wordly gods. So Luatex, which slows down any legacy document, is definitely not on my roadmap anyway. – user4686 Jun 2 '17 at 14:38
  • @jfbu Don't be so harsh on the LuaTeX guys. I think they are doing a good job and LuaTeX is currently the only engine with progressing development. It's not their fault that Lua has become crap recently. – Henri Menke Jun 3 '17 at 2:03
  • tug.org/pipermail/luatex/2017-June/006533.html – user4686 Jun 3 '17 at 5:54
17

I would do this using an expandable implementation: it's easy enough to keep everything on the stack:

\documentclass{article}
\usepackage{expl3}
\usepackage{xparse}
\begin{document}
\ExplSyntaxOn
\cs_new:Npn \collatz_next:n #1
  {
    \int_if_even:nTF {#1}
      { \int_eval:n { #1 / 2 } }
      { \int_eval:n { 3 * #1 + 1 } }
  }
\cs_generate_variant:Nn  \collatz_next:n { f }
\cs_new:Npn \collatz_count:n #1
  { \__collatz_count:nn {#1} { 1 } }
\cs_new:Npn \__collatz_count:nn #1#2
  {
    \int_compare:nNnTF {#1} > { 1 }
      {
        \__collatz_count:ff
          { \collatz_next:n {#1} } { \int_eval:n { #2 + 1 } }
      }
      { #2 }
  }
\cs_generate_variant:Nn \__collatz_count:nn { ff }
\cs_new:Npn \collatz_max_range:nn #1#2
  {
    \__collatz_max_range:nnnn { 1 } {#1} {#1} {#2}
  }
\cs_new:Npn \__collatz_max_range:nnnn #1#2#3#4
  {
    \int_compare:nNnTF {#3} > {#4}
      {#2} % Could also dump length (#1) here
      {
        \__collatz_max_range:fnnnn { \collatz_count:n {#3} } 
          {#1} {#2} {#3} {#4}
      }
  }
\cs_generate_variant:Nn \__collatz_max_range:nnnn { nnf }
\cs_new:Npn \__collatz_max_range:nnnnn #1#2#3#4#5
  {
    \int_compare:nNnTF {#1} > {#2}
      { \__collatz_max_range:nnfn {#1} {#4} }
      { \__collatz_max_range:nnfn {#2} {#3} }
        { \int_eval:n { #4 + 1 } } {#5}
  }
\cs_generate_variant:Nn \__collatz_max_range:nnnnn { f }
\NewExpandableDocumentCommand \collatzcount { m }
  { \collatz_count:n {#1} }
\NewExpandableDocumentCommand \collatzmaxinrange { mm }
  {
    \collatz_max_range:nn {#1} {#2}
  }

\ExplSyntaxOff

\collatzcount{13}
\collatzmaxinrange{1}{13}

\end{document}

Here I'm using the fact that f-type expansion is itself expandable: it stops at the first non-expandable token, which is ideal here as it's the number we want as a result.

As I've noted in the code, one might want to arrange to 'return' two values: that might be done by adding an argument such as

\cs_new:Npn \collatz_max_range:nnN #1#2#3
  {
    \__collatz_max_range:nnnnN { 1 } {#1} {#1} {#2} #3
  }
\cs_new:Npn \__collatz_max_range:nnnnN #1#2#3#4#5
  {
    \int_compare:nNnTF {#3} > {#4}
      { #5 {#1} {#2}} % Could also dump length (#1) here
      {
        \__collatz_max_range:fnnnnN { \collatz_count:n {#3} } 
          {#1} {#2} {#3} {#4}
      }
  }
\cs_generate_variant:Nn \__collatz_max_range:nnnnN { nnf }
\cs_new:Npn \__collatz_max_range:nnnnnN #1#2#3#4#5#6
  {
    \int_compare:nNnTF {#1} > {#2}
      { \__collatz_max_range:nnfnN {#1} {#4} }
      { \__collatz_max_range:nnfnN {#2} {#3} }
        { \int_eval:n { #4 + 1 } } {#5} #6
  }

where the third argument to \collatz_max_range:nnN might be something as simple as \use_i:nn or something much more complex (store the values, typeset both with filler text, etc.)


One might also consider using a 'pre-calculation' approach to save data about each value. For example, we might build up data on which is the 'next' value to pick: this could be used to display either the chain of values or the chain length:

\prop_new:N \g__collatz_next_prop
\int_new:N \l__collatz__tmp_int
\cs_new_protected:Npn \collatz_precalc:n #1
  {
    \int_set:Nn \l__collatz__tmp_int { \collatz_next:n {#1} }
    \prop_gput:NnV \g__collatz_next_prop
      {#1} \l__collatz__tmp_int
    \prop_if_in:NVF \g__collatz_next_prop \l__collatz__tmp_int
      { \collatz_precalc:V  \l__collatz__tmp_int }
  }
\cs_generate_variant:Nn \collatz_precalc:n { V }
\int_step_inline:nnnn { 1 } { 1 } { 100 }
  { \collatz_precalc:n {#1} }

One might consider using one tl (or int) per value, but that gets expensive in terms of csname usage (we are also somewhat limited on int registers.) Depending on the real use case, having either pre-calcuated values or a pre-cacluated set of chains may be useful: depends on whether the data are repeatedly reused.

5
  • One could arrange to output the result as well as the length of the chain, perhaps passing to an auxiliary 'typeset the result' function. Will depend on a real use case ... – Joseph Wright May 31 '17 at 20:08
  • @jfbu Ah, you mean that he wants the longest sequence? That's just a question of doing a stepwise loop over the possible input values: I'd imagine in a real use case one would want to actually typeset all of the results and (probably) the chains too. It's not really clear to me what is wanted. – Joseph Wright Jun 1 '17 at 9:31
  • @jfbu I've addressed that point and added some more possible things one might do (not part of the question but part of a solution to the more general idea). – Joseph Wright Jun 1 '17 at 14:13
  • In my answer I had a comment about using \fontdimen parameters of a dummy font (at 1sp, say) for storage. But I have forgotten the details and limitations of that. And I think this is probably the underlying technique for the l3 arrays, right ? – user4686 Jun 1 '17 at 14:16
  • @jfbu Yes, Bruno has recently added l3intarray for internal use in l3regex using this trick. At this stage we are wary of making it a public interface as there are some restrictions: you have to know the array size up-front and you have to limit values to have of the normal integer limit. On the other hand, it is fast and at least might be useful more generally: raise perhaps on LaTeX-L. – Joseph Wright Jun 1 '17 at 14:19
17

With Plain TeX (and eTeX for convenience)

  1. initial approach limited by TeX arithmetic up to N = 159486 as maximal starting point.

  2. enlarged approach doing big arithmetic and able to master say N = 1,000,000

  3. technical variant using the technique of font dimensions to store all the Collatz lengths with start point from 1 up to N. Maximal N is 5,000,000 (on my TeXLive I could go up to a bit less than 8,000,000. Beyond needs enlarging TeX memory). At the end of the day \number\fontdimen<number>\cz has stored all the Collatz lengths with <number> as starting point from 1 to N=5,000,000.

  4. finally a somewhat simpler algorithm brings close to a 2x speed improvement.


edited my initial post was not at all executing that part of it which handled the storage of known things in control sequences. Hence it was far slower, but on the other hand, when I corrected it, I observed it quickly hit TeX memory limit, due to enclosing expansion in \begingroup and \endgroup. But, after all we are storing absolute data, so no need for a group. And, I decided to store only data with indices up to the maximal N, in order again to escape memory limitations.

% pdftex (or etex)

\newcount\inputNmax

\newcount\intN
\newcount\intNa
\newcount\intNtop
\newcount\intL
\newcount\intLtop

\def\CollOne{%
    \advance\intNa 1
    \ifnum\intNa > \inputNmax
       \CollDone
    \else
      \intL 0 % will store number of steps starting at Na
      \intN = \intNa
      \expandafter\CollTwo
    \fi
}

% store in 
\def\CollTwo{%
    \ifcsname collatz\the\intN\endcsname
       \expandafter\CollThreeA
    \else
       \expandafter\CollThreeB
    \fi
}

\def\CollThreeA{%
    \advance\intL\csname collatz\the\intN\endcsname\relax
    \expandafter\edef\csname collatz\the\intNa\endcsname{\the\intL}%
    \ifnum\intL > \intLtop
          \intLtop = \intL
          \intNtop = \intNa
    \fi
    \intN = \intNa
    \CollUpdate
}

\def\CollUpdate{%
    \advance\intL -1
    \ifodd\intN
      \multiply\intN 3
      \advance\intN  1
    \else
      \divide\intN 2
    \fi
    \let\next\CollOne
    \ifnum\intN>\inputNmax
    \else
     \ifcsname collatz\the\intN\endcsname
     \else
        \expandafter\edef\csname collatz\the\intN\endcsname{\the\intL}%
        \let\next\CollUpdate
     \fi
    \fi
    \next
%% this variant seems to impact a bit negatively execution time
% \ifcsname collatz\the\intN\endcsname
%   \expandafter\CollOne
% \else
%   \ifnum\intN>\inputNmax
%   \else
%     \expandafter\edef\csname collatz\the\intN\endcsname{\the\intL}%
%   \fi
%   \expandafter\CollUpdate
% \fi
} 

\def\CollThreeB{%
    \advance\intL 1
    \ifodd\intN
      \multiply\intN 3
      \advance\intN  1
    \else
      \divide\intN 2
    \fi
    \CollTwo
}

\def\CollDone{%
    From 1 to \the\inputNmax, the longest sequence with smallest starting
    point was observed to start at \the\intNtop, and contained
    \the\intLtop\relax\
    elements.\par
}

\def\CollMax #1{% #1 integer at least 1
    %\begingroup
      \inputNmax=#1\relax
      \intNa   = 1
      \intNtop = 1
      \intLtop = 1
    % I would prefer counting steps to reach 1, so here 0
    % but it seems the question asks for number of elements, so here 1
      \expandafter\def\csname collatz1\endcsname{1}%
      \CollOne
    %\endgroup
}


\hsize10cm

\CollMax {10}

\CollMax {100}

\CollMax {1000}

\CollMax {10000}

\CollMax {100000}

\bye

Output:

enter image description here

Execution time:

$ time pdftex -interaction batchmode pcollatz.tex
This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017) (preloaded format=pdftex)
 restricted \write18 enabled.
entering extended mode

real    0m6.667s
user    0m6.636s
sys 0m0.027s

There is arithmetic overflow for 1000000.

I will add variant using xint at some later stage.

EXTENSION TO BIG INTEGERS

Actually 159487 is the smallest starting point leading to a TeX arithmetic overflow. (its 56th iterate is 2265333694 >= 2**31)

But then, after having modified the code above to use xintcore macros for big arithmetic, there was another problem which was that the "TeX pool" gets filled-up for starting point at N=317893. I did run it with 317892 obtaining

\CollMax {317892}
From 1 to 317892, the longest sequence with smallest starting point was
observed to start at 230631, and contained 443 elements.

Then I modified the code and got it to store the sequence lengths in memory only for the first 100000 integers, thus overcoming that "TeX pool" limitation. Still handling big arithmetic via xintcore macros.

% pdftex (or etex)

\input xintcore.sty % to handle big integers.
% there is also bnumexpr but it is LaTeX interface only
% or xint, xintexpr which are more extensive than xintcore

\newcount\inputNmax

\newcount\intN
\newcount\intNa
\newcount\intNtop
\newcount\intL
\newcount\intLtop

\def\CollLoop{%
    \advance\intNa 1
% \immediate\write-1{\the\intNa}%
    \ifnum\intNa > \inputNmax
       \CollDone
    \else
      \intL 0 % will store number of steps starting at Na
      \intN = \intNa
      \expandafter\CollTwo
    \fi
}

\def\CollTwo{%
    \ifcsname collatz\the\intN\endcsname
       \expandafter\CollThreeA
    \else
       \expandafter\CollThreeB
    \fi
}

% \CollThreeA will be either \CollThreeAwithUpdate or \CollThreeAwithNoUpdate
\def\CollThreeAwithUpdate{%
    \advance\intL\csname collatz\the\intN\endcsname\relax
    \expandafter\edef\csname collatz\the\intNa\endcsname{\the\intL}%
    \ifnum\intL > \intLtop
          \intLtop = \intL
          \intNtop = \intNa
    \fi
    \intN = \intNa
    \CollUpdate
}

\def\CollUpdate{%
    \advance\intL -1
    \ifodd\intN
      \multiply\intN 3
      \advance\intN 1
    \else
      \divide\intN 2
    \fi
    \let\next\CollLoop
    \ifnum\intN>\inputNmax
    \else
     \ifcsname collatz\the\intN\endcsname
     \else
        \expandafter\edef\csname collatz\the\intN\endcsname{\the\intL}%
        \let\next\CollUpdate
     \fi
    \fi
    \next
} 

\def\CollThreeAwithNoUpdate{%
    \advance\intL\csname collatz\the\intN\endcsname\relax
    \ifnum\intL > \intLtop
          \intLtop = \intL
          \intNtop = \intNa
    \fi
    \CollLoop
}

\def\CollThreeB{%
    \let\next\CollTwo
    \ifodd\intN
      \ifnum\intN>\maxdimen
         % notice that necessarily this first happens with previous execution
         % had done (3x+1)/2, so the real antecedent was > "7FFFFFFF
         % and would have created arithmetic overflow if we had done
         % x->3x+1->(3x+1)/2
         \edef\bigintN{\the\intN}%
         \let\next\CollThreeBig
      \else
         \advance\intL 2
         \divide\intN 2
         \multiply\intN 3
         \advance\intN 2
      \fi
    \else
      \advance\intL 1
      \divide\intN 2
    \fi
    \next
}%

% \def\error{\immediate\write-1{\the\intNa, \the\intL}\csname end\endcsname}

\def\CollThreeBig{% 
% 159487 is the smallest starting integer which triggers this, as its
% 56th iterate 2265333694 exceeds 2**31, and the 57th is thus > \maxdimen
% \error
    \advance\intL 1
    % \xintLastItem does no expansion ...
    \ifodd\expandafter\xintLastItem\expandafter{\bigintN}
      \advance\intL 1
      \edef\bigintN{\xintHalf{\xintiiMul{\bigintN}3}}% Half truncates
      % possibly faster to use \xintDouble and an addition, not tested
    \else
      \edef\bigintN{\xintHalf{\bigintN}}%
    \fi
    % \xintLength does no expansion ...
    \ifnum\expandafter\xintLength\expandafter{\bigintN}>9
        \expandafter\CollThreeBig
    \else
        \intN = \bigintN\relax
        \expandafter\CollThreeB
    \fi
}%


\def\CollReport{%
    From 1 to \the\inputNmax, the longest sequence with smallest starting
    point was observed to start at \the\intNtop, and contained
    \the\intLtop\relax\
    elements.\par
}

\let\CollDone\CollReport

\def\CollMaxInitial {%
      \let\CollThreeA\CollThreeAwithUpdate % faster but uses macro storage
      \intNa   = 1
      \intNtop = 1
% I would prefer counting steps to reach 1, so here 0
% but it seems the question asks for number of elements, so 1
% (and not 4 although 1->4->2->1, as I consider that 1 is sequence in itself)
      \expandafter\def\csname collatz1\endcsname{1}%
      \intLtop = 1
      \CollLoop
}

\def\CollMax #1{% #1 integer at least 1
      \ifnum#1>100000
          \inputNmax 100000
          \let\CollDone\empty
          \CollMaxInitial
          \let\CollDone\CollReport
          \let\CollThreeA\CollThreeAwithNoUpdate
          \inputNmax=#1\relax
          \intNa = 100000
          \CollLoop
      \else
          \inputNmax = #1\relax
          \CollMaxInitial
      \fi
}


\hsize10cm

% \CollMax {10}

% \CollMax {100}

% \CollMax {1000}

% \CollMax {10000}

% \CollMax {100000}

% \CollMax {200000}
% From 1 to 200000, the longest sequence with smallest starting
% point was observed to start at 156159, and contained 383 elements.

\CollMax {1000000}
% From 1 to 1000000, the longest sequence with smallest start-
% ing point was observed to start at 837799, and contained 525
% elements.
\bye

This gives

enter image description here

and execution time is

$ time pdftex -interaction batchmode pcollatz-big.tex
This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017) (preloaded format=pdftex)
 restricted \write18 enabled.
entering extended mode

real    0m40.345s
user    0m40.247s
sys 0m0.087s

edit: This is on a computer which is usually about 15% faster than the one used for testing first version of the code with 100000. But trying now this code with N=1000000 and the slower computer I observe about 1m execution time (thus +50%...). For some reasons I don't know pdftex becomes indeed more significantly slower on my laptop when TeX memory is used a lot. I have observed it in the past when testing xint on computations with thousands of digits.

USE OF \fontdimen PARAMETERS

Here we go

% pdftex (or etex)

\input xintcore.sty % to handle big integers.
% there is also bnumexpr but it is LaTeX interface only
% or xint, xintexpr which are more extensive than xintcore

% use fontdimen parameters to have an array where to store
% "Collatz lengths". At the end of the day
% \number\fontdimen<number>\cz gives the length of the
% sequence starting at <number> (for all numbers up to 5,000,000) and
% reaching 1 (included).

\newcount\inputNmax

\newcount\intN
\newcount\intNa
\newcount\intNtop
\newcount\intL
\newcount\intLtop

\newcount\czsize
\czsize 5000000
\font\cz=cmr10 at 1pt
\fontdimen \czsize\cz = 0sp % make room ...
% vz  texmf.cnf
% Words of font info for TeX (total size of all TFM files, approximately).
% Must be >= 20000 and <= 147483647 (without tex.ch changes).
% font_mem_size = 8000000


% make sure array entries are zero (only \fontdimen2 to 7 are populated for cmr10)
\intN 1
\loop
  \fontdimen\intN\cz = 0sp
  \advance\intN 1
\ifnum\intN < 8
\repeat

% do I need to do that for all, aren't they zero except first few ones ?
% no, it's ok (paranoide check here, done once)
% \intN 1
% \loop
%   \ifnum\fontdimen\intN\cz>0 \error\fi
% \ifnum\intN < \czsize
%   \advance\intN 1
% \repeat

\def\CollLoop{%
    \advance\intNa 1
    \ifnum\intNa > \inputNmax
       \CollDone
    \else
      \intL 0 % will store number of steps starting at Na
      \intN = \intNa
      \expandafter\CollTwo
    \fi
}

\def\CollTwo{%
  \let\next\CollThreeB
  \unless\ifnum\intN>\czsize
    \ifnum\fontdimen\intN\cz>0
       \let\next\CollThreeA
    \fi
  \fi
  \next
}

% \CollThreeA will be either \CollThreeAwithUpdate or \CollThreeAwithNoUpdate
% but we are going to use it always with update...
\def\CollThreeAwithUpdate{%
    \advance\intL\fontdimen\intN\cz
    \fontdimen\intNa\cz=\intL sp
    \ifnum\intL > \intLtop
          \intLtop = \intL
          \intNtop = \intNa
    \fi
    \intN = \intNa
    \CollUpdate
}

\def\CollUpdate{%
    \advance\intL -1
    \ifodd\intN
      \multiply\intN 3
      \advance\intN 1
    \else
      \divide\intN 2
    \fi
    \let\next\CollLoop
    \ifnum\intN>\inputNmax % always at most \cssize in this macro
    \else
     \ifnum\fontdimen\intN\cz=0
        \fontdimen\intN\cz=\intL sp
        \let\next\CollUpdate
     \fi
    \fi
    \next
} 

\def\CollThreeAwithNoUpdate{%
    \advance\intL\fontdimen\intN\cz
    \ifnum\intL > \intLtop
          \intLtop = \intL
          \intNtop = \intNa
    \fi
    \CollLoop
}

\def\CollThreeB{%
    \let\next\CollTwo
    \ifodd\intN
      \ifnum\intN>\maxdimen
         % notice that necessarily this first happens with previous execution
         % had done (3x+1)/2, so the real antecedent was > "7FFFFFFF
         % and would have created arithmetic overflow if we had done
         % x->3x+1->(3x+1)/2
         \edef\bigintN{\the\intN}%
         \let\next\CollThreeBig
      \else
         \advance\intL 2
         \divide\intN 2
         \multiply\intN 3
         \advance\intN 2
      \fi
    \else
      \advance\intL 1
      \divide\intN 2
    \fi
    \next
}%

% \def\error{\immediate\write-1{\the\intNa, \the\intL}\csname end\endcsname}

\def\CollThreeBig{% 
% 159487 is the smallest starting integer which triggers this, as its
% 56th iterate 2265333694 exceeds 2**31, and the 57th is thus > \maxdimen
% \error
    \advance\intL 1
    % \xintLastItem does no expansion ...
    \ifodd\expandafter\xintLastItem\expandafter{\bigintN}
      \advance\intL 1
      \edef\bigintN{\xintHalf{\xintiiMul{\bigintN}3}}% Half truncates
      % possibly faster to use \xintDouble and an addition, not tested
    \else
      \edef\bigintN{\xintHalf{\bigintN}}%
    \fi
    % \xintLength does no expansion ...
    \ifnum\expandafter\xintLength\expandafter{\bigintN}>9
        \expandafter\CollThreeBig
    \else
        \intN = \bigintN\relax
        \expandafter\CollThreeB
    \fi
}%


\def\CollReport{%
    From 1 to \the\inputNmax, the longest sequence with smallest starting
    point was observed to start at \the\intNtop, and contained
    \the\intLtop\relax\
    elements.\par
}

\let\CollDone\CollReport

\def\CollMaxInitial {%
      \let\CollThreeA\CollThreeAwithUpdate % limited by font array storage
      \intNa   = 1
      \intNtop = 1
% I would prefer counting steps to reach 1, so here 0
% but it seems the question asks for number of elements, so 1
% (and not 4 although 1->4->2->1, as I consider that 1 is sequence in itself)
      \fontdimen1 \cz= 1sp
      \intLtop = 1
      \CollLoop
}

\def\CollMax #1{% #1 integer at least 1
      \ifnum#1>\czsize
          \inputNmax \czsize
          \let\CollDone\empty
          \CollMaxInitial
          \let\CollDone\CollReport
          \let\CollThreeA\CollThreeAwithNoUpdate
          \inputNmax=#1\relax
          \intNa = \czsize
          \CollLoop
      \else
          \inputNmax = #1\relax
          \CollMaxInitial
      \fi
}


\hsize10cm

% check it does give same results as earlier...
% \CollMax {10}

% \CollMax {100}

% \CollMax {1000}

% \CollMax {10000}

%\CollMax {100000}

%\CollMax {1000000}

\CollMax {5000000}

% From 1 to 5000000, the longest sequence with smallest starting point was
% observed to start at 2929311, and contained 550 elements.
\bye

produces:

enter image description here

For 1000000 it is about 4x-6x faster than method with csname (hard to tell exactly because depends on which computer I test with). Perhaps I can improve more, I am not too much used to that technique.

Here is with my laptop (the "slow" computer) for N=1,000,000:

$ time pdftex --interaction=batchmode pcollatz-array.tex
This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017) (preloaded format=pdftex)
 restricted \write18 enabled.
entering extended mode

real    0m10.408s
user    0m10.313s
sys 0m0.041s

For 5,000,000 it takes about five times as much as for 1,000,000.


An ultimate rewrite brings a 45% speed improvement.

\input xintcore.sty

\newcount\inputNmax
\newcount\inputNstop

\newcount\intN
\newcount\intNa
\newcount\intNb
\newcount\intNtop
\newcount\intL
\newcount\intLtop

\newcount\czsize
\czsize 5000000
\font\cz=cmr10 at 1pt
\fontdimen \czsize\cz = 0sp

\intN 1
\loop
  \fontdimen\intN\cz = 0sp
  \advance\intN 1
\ifnum\intN < 8
\repeat

\catcode`@ 11

\long\def\@gobble#1{}
\edef\sp{\expandafter\@gobble\string\sp}% (no space needed?)

\def\CollLoopA{%
    \advance\intNa 1
    \intNb \intNa
    \divide\intNb 2
    \intL\fontdimen\intNb\cz
    \advance\intL 1
    \fontdimen\intNa\cz=\intL\sp
    \advance\intNa 1
    \intNb \intNa
    \divide\intNb 3
    \multiply\intNb 2
    \advance\intNb 1
    \intL\fontdimen\intNb\cz
    \advance\intL -2
    \fontdimen\intNa\cz=\intL\sp
    \advance\intNa 1        
    \ifnum\intNa > \inputNstop
       \CollDone
    \else
      \let\CollLoopBack\CollLoopB
      \intL 0
      \intN = \intNa
      \expandafter\CollThreeB
    \fi
}

\def\CollLoopB{%
    \advance\intNa 1
    \ifnum\intNa > \inputNstop
       \CollDone
    \else
      \let\CollLoopBack\CollLoopC
      \intL 0
      \intN = \intNa
      \expandafter\CollThreeB
    \fi
}

\def\CollLoopC{%
    \advance\intNa 1
    \intNb \intNa
    \divide\intNb 2
    \intL\fontdimen\intNb\cz
    \advance\intL 1
    \fontdimen\intNa\cz=\intL\sp
    \advance\intNa 1
    \ifnum\intNa > \inputNstop
       \CollDone
    \else
      \let\CollLoopBack\CollLoopA
      \intL 0
      \intN = \intNa
      \expandafter\CollThreeB
    \fi
}

\def\CollTwo{%
  \ifnum\intN<\intNa
    \advance\intL\fontdimen\intN\cz
    \fontdimen\intNa\cz=\intL\sp
    \CollCheckTop
    \expandafter\CollLoopBack
  \else
    \expandafter\CollThreeB 
  \fi
}

\def\CollCheckTopyes {%
    \ifnum\intL > \intLtop
          \intLtop = \intL
          \intNtop = \intNa
    \fi
}%
\let\CollCheckTopno\empty
\let\CollCheckTop\CollCheckTopyes

\def\CollThreeB{%
    \let\next\CollTwo
    \ifodd\intN
      \ifnum\intN>\maxdimen
         \edef\bigintN{\the\intN}%
         \let\next\CollThreeBig
      \else
         \advance\intL 2
         \divide\intN 2
         \multiply\intN 3
         \advance\intN 2
      \fi
    \else
      \advance\intL 1
      \divide\intN 2
    \fi
    \next
}%

\def\CollThreeBig{% 
    \advance\intL 1
    % \xintLastItem does no expansion ...
    \ifodd\expandafter\xintLastItem\expandafter{\bigintN}
      \advance\intL 1
      \edef\bigintN{\xintHalf{\xintiiMul{\bigintN}3}}% Half truncates
      % possibly faster to use \xintDouble and an addition, not tested
    \else
      \edef\bigintN{\xintHalf{\bigintN}}%
    \fi
    % \xintLength does no expansion ...
    \ifnum\expandafter\xintLength\expandafter{\bigintN}>9
        \expandafter\CollThreeBig
    \else
        \intN = \bigintN\relax
        \expandafter\CollThreeB
    \fi
}%


\def\CollReport{%
    From 1 to \the\inputNmax, the longest sequence with smallest starting
    point was observed to start at \the\intNtop, and contained
    \the\intLtop\relax\
    elements.\par
}

\let\CollDone\empty

\def\CollMax #1{%
      \inputNmax #1\relax
      \inputNstop \inputNmax
      \divide\inputNstop 2
      \intNa = 1
      \intNtop = 1
      \fontdimen1 \cz= 1sp
      \intLtop = 1
      \let\CollCheckTop\CollCheckTopno
      \CollLoopC
      \inputNstop \inputNmax
      \let\CollCheckTop\CollCheckTopyes
      \ifcase \numexpr3+\intNa - 6*(\intNa/6)\relax
         \advance\intNa-2 \let\next\CollLoopC\or
         \error\or
         \error\or
         \advance\intNa-3 \let\next\CollLoopA\or
         \advance\intNa-1 \let\next\CollLoopB\else
         \error
      \fi 
      \next
      \CollReport
}
\catcode`@ 12

\hsize10cm

\CollMax {1000000}

\bye

With time pdftex --interaction=batchmode pcollatz-arrayIII.tex I get a user time of 0m5.780s compared to the former 0m10.313s. On @HenriMenke's machine, execution time of this pure eTeX approach should be around 2s (for N=1000000).

8
  • If one is only interested into obtaining the iterates of the 3x+1 sequence, in a fully expandable way, there is a one-liner in the xint manual: \def\syr #1{\xinttheiiexpr rseq(#1; @<=1)?{break(i)}{odd(@)?{3@+1}{@//2}},i=0++)\relax}. One needs \usepackage{xintexpr} in LaTeX or \input xintexpr.sty in TeX. Then for example \edef\x{\syr{1000}} and \message{\meaning\x} or just \x. – user4686 Jun 1 '17 at 14:11
  • it is very strange that the code with fontdimen's has about the same execution time both on my laptop and the other machine, be it with N=1,000,000 or N=5,000,000, but the one with csname's takes about 40s on the other machine for N=1,000,000 to be compared with more than 1mn on my laptop. I have tested and re-tested... – user4686 Jun 2 '17 at 7:41
  • Thanks for the perfect example, I even unsuspected that this can be done on plain TeX – sergiokapone Jun 2 '17 at 10:36
  • @jfbu As you were interested in the stats of your fontdimen solution on my desktop machine, here some outputs: gist.github.com/anonymous/163040b9139e8324ec63ad79f4c80331 – Henri Menke Jun 3 '17 at 2:01
  • @HenriMenke thanks for the timing data. This seems to indicate that your CPU is of the order of 2.5x--3x faster than mine (2 GHz Intel Core i7 with 8 Go 1600 MHz DDR3, 2012 mac laptop running under 10.9.5). Just an order of magnitude indication. (it might be rather 2.5x in general, but is a 2.8x in this specific instance with pdftex; assuming the pdftex binary is compiled with same flags on both machines) – user4686 Jun 7 '17 at 7:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.