4

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?

7
  • What do you want to do with the small numbers? And can you include only the use case you wish to have? IT seems that there is some other problem going on behind the scenes.
    – percusse
    Jun 20 '14 at 12:08
  • @percusse I was afraid the full extent of the problem would cause people running for the hills. Should I post a seperate question or edit this existing one?
    – 1010011010
    Jun 20 '14 at 12:27
  • I was curious about the \qa,\qb,\eval@limit etc. If you can provide just the small number array and what you would like to have, then probably it is easier to get to the bottom of the problem. It is kind of distracting now with all dead ends present and I didn't get the foreach part.
    – percusse
    Jun 20 '14 at 12:40
  • @percusse I'll edit when I get home. Give me an hour or so. :-)
    – 1010011010
    Jun 20 '14 at 12:45
  • @percusse I've updated the question.
    – 1010011010
    Jun 20 '14 at 14:26
8
+100

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

Result

5
  • It is nice, though one question: Will this fully work with existing pdflatex code? As in, this will work with any packages that worked without LuaLaTeX?
    – 1010011010
    Jun 22 '14 at 15:27
  • @1010011010 In principle, lualatex should be mostly compatible with pdflatex, but with more capabilities. However: 1) input should be in utf8 and no inputenc package is allowed. 2) Font management is different in luatex, so perhaps some packages wich install their own fonts don't work. 3) The last word is in the experiment: you better try to compile any of your documents with lualatex instead of pdflatex and see if it works.
    – JLDiaz
    Jun 22 '14 at 15:31
  • @1010011010 more info: tex.stackexchange.com/questions/28642/…
    – JLDiaz
    Jun 22 '14 at 15:38
  • My document compiles as normally with LuaLaTeX. :-) When the bounty expires I'll also check your answer.
    – 1010011010
    Jun 22 '14 at 16:32
  • @JLDiaz The comment about the pasted image is beautiful! Jun 23 '14 at 8:18

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.