A TeX translation of old C code - Evaluating elements of an array *inside* a foreach loop

Question

I'm trying to translate a conceptually simple part of an algorithm, from the original source for the R-engine (written in C) to TeX. This shouldn't take more than some logic trees and for/while loops.

Per @percusse's comment I'll explain the full use case to understand what I want to do.

The code I wish to mimic is the following:

int attribute_hidden chebyshev_init(double *dos, int nos, double eta)
{
    int i, ii;
    double err;

    if (nos < 1)
        return 0;

    err = 0.0;
    i = 0;  /* just to avoid compiler warnings */
    for (ii=1; ii<=nos; ii++) {
        i = nos - ii;
        err += fabs(dos[i]);
        if (err > eta) {
        return i;
        }
    }
    return i;
}

Source: https://github.com/wch/r-source/blob/776708efe6003e36f02587ad47b2eaaaa19e2f69/src/nmath/chebyshev.c#L48 (Lines 48-66)

The real use case can be found at https://github.com/wch/r-source/blob/0b493ddbefbc3cc06f675841c4672b3197dbfbbb/src/nmath/gamma.c#L103. The arguments show:

ngam = chebyshev_init(gamcs, 42, DBL_EPSILON/20)

Gamcs is defined as:

const static double gamcs[42] = {
+.8571195590989331421920062399942e-2,
+.4415381324841006757191315771652e-2,
+.5685043681599363378632664588789e-1,
-.4219835396418560501012500186624e-2,
+.1326808181212460220584006796352e-2,
-.1893024529798880432523947023886e-3,
+.3606925327441245256578082217225e-4,
-.6056761904460864218485548290365e-5,
+.1055829546302283344731823509093e-5,
-.1811967365542384048291855891166e-6,
+.3117724964715322277790254593169e-7,
-.5354219639019687140874081024347e-8,
+.9193275519859588946887786825940e-9,
-.1577941280288339761767423273953e-9,
+.2707980622934954543266540433089e-10,
-.4646818653825730144081661058933e-11,
+.7973350192007419656460767175359e-12,
-.1368078209830916025799499172309e-12,
+.2347319486563800657233471771688e-13,
-.4027432614949066932766570534699e-14,
+.6910051747372100912138336975257e-15,
-.1185584500221992907052387126192e-15,
+.2034148542496373955201026051932e-16,
-.3490054341717405849274012949108e-17,
+.5987993856485305567135051066026e-18,
-.1027378057872228074490069778431e-18,
+.1762702816060529824942759660748e-19,
-.3024320653735306260958772112042e-20,
+.5188914660218397839717833550506e-21,
-.8902770842456576692449251601066e-22,
+.1527474068493342602274596891306e-22,
-.2620731256187362900257328332799e-23,
+.4496464047830538670331046570666e-24,
-.7714712731336877911703901525333e-25,
+.1323635453126044036486572714666e-25,
-.2270999412942928816702313813333e-26,
+.3896418998003991449320816639999e-27,
-.6685198115125953327792127999999e-28,
+.1146998663140024384347613866666e-28,
-.1967938586345134677295103999999e-29,
+.3376448816585338090334890666666e-30,
-.5793070335782135784625493333333e-31
    }

DBL_EPSILON is defined as DBL_EPSILON = 2^(-52).

To know what result I'm looking for, we can mimic the entire calculation of chebyshev_init. To do that, you can run the following code in R:

DBL_EPSILON = 2^(-52)
err = 0.0
i = 0
dos = c(+.8571195590989331421920062399942e-2,
+.4415381324841006757191315771652e-2,
+.5685043681599363378632664588789e-1,
-.4219835396418560501012500186624e-2,
+.1326808181212460220584006796352e-2,
-.1893024529798880432523947023886e-3,
+.3606925327441245256578082217225e-4,
-.6056761904460864218485548290365e-5,
+.1055829546302283344731823509093e-5,
-.1811967365542384048291855891166e-6,
+.3117724964715322277790254593169e-7,
-.5354219639019687140874081024347e-8,
+.9193275519859588946887786825940e-9,
-.1577941280288339761767423273953e-9,
+.2707980622934954543266540433089e-10,
-.4646818653825730144081661058933e-11,
+.7973350192007419656460767175359e-12,
-.1368078209830916025799499172309e-12,
+.2347319486563800657233471771688e-13,
-.4027432614949066932766570534699e-14,
+.6910051747372100912138336975257e-15,
-.1185584500221992907052387126192e-15,
+.2034148542496373955201026051932e-16,
-.3490054341717405849274012949108e-17,
+.5987993856485305567135051066026e-18,
-.1027378057872228074490069778431e-18,
+.1762702816060529824942759660748e-19,
-.3024320653735306260958772112042e-20,
+.5188914660218397839717833550506e-21,
-.8902770842456576692449251601066e-22,
+.1527474068493342602274596891306e-22,
-.2620731256187362900257328332799e-23,
+.4496464047830538670331046570666e-24,
-.7714712731336877911703901525333e-25,
+.1323635453126044036486572714666e-25,
-.2270999412942928816702313813333e-26,
+.3896418998003991449320816639999e-27,
-.6685198115125953327792127999999e-28,
+.1146998663140024384347613866666e-28,
-.1967938586345134677295103999999e-29,
+.3376448816585338090334890666666e-30,
-.5793070335782135784625493333333e-31)
nos = 42
eta = DBL_EPSILON/20
for (ii in 1:nos) {
i = nos - ii
err = err + abs(dos[i])
if (err>eta)
{
break
}
}

The resulting value of i is 23, which is indeed the case if err=>eta. We've namely crossed the intended value by one by setting it larger (>) instead of larger or equal to (=>). The intended result also shows at https://github.com/wch/r-source/blob/0b493ddbefbc3cc06f675841c4672b3197dbfbbb/src/nmath/gamma.c#L115. Yes, I'm aware that knowing the result makes this a rather useless excercise, but this case is explanatory and the algorithm is required in more than just this case.

So R in fact does give the correct result. I've made my own version of this algorithm, which would ideally would look a little like this:

\documentclass{article}
\usepackage{tikz}
\usetikzlibrary{fpu}
\usepackage{xintexpr}
\def\gamcs{{
+.8571195590989331421920062399942e-2,
+.4415381324841006757191315771652e-2,
+.5685043681599363378632664588789e-1,
-.4219835396418560501012500186624e-2,
+.1326808181212460220584006796352e-2,
-.1893024529798880432523947023886e-3,
+.3606925327441245256578082217225e-4,
-.6056761904460864218485548290365e-5,
+.1055829546302283344731823509093e-5,
-.1811967365542384048291855891166e-6,
+.3117724964715322277790254593169e-7,
-.5354219639019687140874081024347e-8,
+.9193275519859588946887786825940e-9,
-.1577941280288339761767423273953e-9,
+.2707980622934954543266540433089e-10,
-.4646818653825730144081661058933e-11,
+.7973350192007419656460767175359e-12,
-.1368078209830916025799499172309e-12,
+.2347319486563800657233471771688e-13,
-.4027432614949066932766570534699e-14,
+.6910051747372100912138336975257e-15,
-.1185584500221992907052387126192e-15,
+.2034148542496373955201026051932e-16,
-.3490054341717405849274012949108e-17,
+.5987993856485305567135051066026e-18,
-.1027378057872228074490069778431e-18,
+.1762702816060529824942759660748e-19,
-.3024320653735306260958772112042e-20,
+.5188914660218397839717833550506e-21,
-.8902770842456576692449251601066e-22,
+.1527474068493342602274596891306e-22,
-.2620731256187362900257328332799e-23,
+.4496464047830538670331046570666e-24,
-.7714712731336877911703901525333e-25,
+.1323635453126044036486572714666e-25,
-.2270999412942928816702313813333e-26,
+.3896418998003991449320816639999e-27,
-.6685198115125953327792127999999e-28,
+.1146998663140024384347613866666e-28,
-.1967938586345134677295103999999e-29,
+.3376448816585338090334890666666e-30,
-.5793070335782135784625493333333e-31
    }}
\xintNewExpr\DBLEPSILON{\xintPow{2}{-52}}
\ExplSyntaxOn
\cs_new:Npn \breakloop #1 #2
  {
    \fp_compare:nNnT { #1 } > { #2 }
      { \breakforeach }
  }
\ExplSyntaxOff
\makeatletter
\pgfmathdeclarefunction{chebyshevinit}{3}{%
  \begingroup
    \pgfkeys{/pgf/fpu,/pgf/fpu/output format=fixed}
      \pgfmathparse{#1}%
      \edef\dos{\pgfmathresult}%
      \pgfmathparse{#2-1}% This will allow to refer to the first element as 1 instead of 0
      \edef\nos{\pgfmathresult}%
      \pgfmathparse{#3}%
      \edef\eta{\pgfmathresult}%
        \pgfmathparse{0}%
        \edef\err{\pgfmathresult}%
        \edef\i{\pgfmathresult}%
          \foreach \ii in {1,...,\nos}
          {
            \pgfmathparse{\nos-\ii}%
            \xdef\i{\pgfmathresult}%
            \pgfmathparse{\err}%
            \xdef\err@{\pgfmathresult}%
            \pgfmathparse{\err@ + abs(\dos[\i])}%
            \xdef\err{\pgfmathresult}%
            \breakloop{\err}{\eta}%
          }
        \pgfmathparse{\i}%
    \pgfmath@smuggleone\pgfmathresult%
  \endgroup
}
\makeatother
\begin{document}
\pgfmathparse{chebyshevinit(\gamcs,42,\DBLEPSILON/20)}\pgfmathresult    %
\end{document}

Obviously the above example does not work: Arrays are not an actual "value" but hold multiple values which are expanded in various ways, making it considerably hard to use in a \pgfmathparse setting, even more so inside an actual function (Carlisle, 2014).

I've had a simpler case of Andrew Swann's solution at PGF loop - how to use an array inside pgfmathdeclare function?, but this caused underflow. The fpu library cannot handle the smaller numbers!

Andrew Swann also had the idea to move parts of the loop outside of the loop to be evaluated first in my previous question.

This, however, does not answer my specific question now. Evaluating the value beforehand would not work if I had to evaluate a range of elements (which would be the case now).

So how can I do an evaluation like this one in pgfmath? Solutions outside pgfmath will also work, as long as I can reuse the final result inside pgfmath, tikzpictures, etc., and as long as the end result is at least similar to some digits to the R output. (10^{-4}? I'll go with anything at this point. :-/.)

---

EDIT I have brought this problem back to the essence a bit. This is another example with the same issue: how to evaluate the array inside the loop? Consider this code:

\documentclass{article}
\usepackage{tikz}
\usetikzlibrary{fpu}    
\gdef\myarray{{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42}}

\makeatletter
\pgfmathdeclarefunction{blah}{2}{%
  \begingroup
    \pgfkeys{/pgf/fpu,/pgf/fpu/output format=fixed}
    \pgfmathparse{#1}%
    \edef\x{\pgfmathresult}%
    \pgfmathparse{#2}%
    \edef\n{\pgfmathresult}%
    \pgfmathparse{\x * 2}%
    \edef\twox{\pgfmathresult}%
    \pgfmathparse{0}%
    \edef\b{\pgfmathresult}%
    \edef\b@{\pgfmathresult}%
    \edef\b@@{\pgfmathresult}%
      \foreach \i in {1,...,\n}
      {
      \pgfmathparse{\b@}%
      \xdef\b@@{\pgfmathresult}%
      \pgfmathparse{\b}%
      \xdef\b@{\pgfmathresult}%
      \pgfmathparse{\myarray[\n-\i]}%
      \xdef\arrayvalue{\pgfmathresult}%
      \pgfmathparse{\twox * \b@ - \b@@ + \arrayvalue}%
      \xdef\b{\pgfmathresult}%
      }
    \pgfmathparse{(\b-\b@@) * .5}%
    \pgfmath@smuggleone\pgfmathresult
  \endgroup
}
\makeatother

\begin{document}

\pgfmathparse{blah(5,22)}\pgfmathresult %

\end{document}

So as you can see, Andrew Swann's trick doesn't apply here and I'm struggling to find a workaround. How to use myarray inside my foreach?

Answer 1

Here you have a luatex implementation. It is very straightforward:

\documentclass{article}
% Define the function in lua, and the constant DBL_EPSILON
\directlua{
DBL_EPSILON = math.pow(2,-52)
function chebyshev_init(dos, nos, eta)
    local i, ii, err;
    if nos<1 then
        return 0
    end

    err = 0.0
    i = 0
    for ii=0,nos do
        i = nos - ii
        err = err + math.abs(dos[i])
        if err > eta then
            return i
        end
    end
    return i
end
}

% Define a tex macro to call the lua function
% and "print" the result
\def\chebyshevInit#1#2#3{%
   \directlua{tex.print(chebyshev_init(#1,#2,#3))}
}

% Define example array to test if it works
\xdef\gmcs{{
     .8571195590989331421920062399942e-2,
     .4415381324841006757191315771652e-2,
     .5685043681599363378632664588789e-1,
    -.4219835396418560501012500186624e-2,
     .1326808181212460220584006796352e-2,
    -.1893024529798880432523947023886e-3,
     .3606925327441245256578082217225e-4,
    -.6056761904460864218485548290365e-5,
     .1055829546302283344731823509093e-5,
    -.1811967365542384048291855891166e-6,
     .3117724964715322277790254593169e-7,
    -.5354219639019687140874081024347e-8,
     .9193275519859588946887786825940e-9,
    -.1577941280288339761767423273953e-9,
     .2707980622934954543266540433089e-10,
    -.4646818653825730144081661058933e-11,
     .7973350192007419656460767175359e-12,
    -.1368078209830916025799499172309e-12,
     .2347319486563800657233471771688e-13,
    -.4027432614949066932766570534699e-14,
     .6910051747372100912138336975257e-15,
    -.1185584500221992907052387126192e-15,
     .2034148542496373955201026051932e-16,
    -.3490054341717405849274012949108e-17,
     .5987993856485305567135051066026e-18,
    -.1027378057872228074490069778431e-18,
     .1762702816060529824942759660748e-19,
    -.3024320653735306260958772112042e-20,
     .5188914660218397839717833550506e-21,
    -.8902770842456576692449251601066e-22,
     .1527474068493342602274596891306e-22,
    -.2620731256187362900257328332799e-23,
     .4496464047830538670331046570666e-24,
    -.7714712731336877911703901525333e-25,
     .1323635453126044036486572714666e-25,
    -.2270999412942928816702313813333e-26,
     .3896418998003991449320816639999e-27,
    -.6685198115125953327792127999999e-28,
     .1146998663140024384347613866666e-28,
    -.1967938586345134677295103999999e-29,
     .3376448816585338090334890666666e-30,
    -.5793070335782135784625493333333e-31
}}

\begin{document}
The result is: \chebyshevInit{\gmcs}{42}{DBL_EPSILON/20}
\end{document}

Compile it with lualatex, and you get the following result (sorry, the pasted image is visually dull, but hopefully exciting for the OP :-)