texdimens (10nov21)

git-svn-id: svn://tug.org/texlive/trunk@61018 c570f23f-e606-0410-a88d-b1316a301751

1 files changed, 673 insertions, 264 deletions

diff --git a/Master/texmf-dist/tex/generic/texdimens/texdimens.tex b/Master/texmf-dist/tex/generic/texdimens/texdimens.tex
index 61086d138dd..de6f20e58fe 100644
--- a/Master/texmf-dist/tex/generic/texdimens/texdimens.tex
+++ b/Master/texmf-dist/tex/generic/texdimens/texdimens.tex
@@ -1,170 +1,207 @@
 % This is file texdimens.tex, part of texdimens package, which
 % is distributed under the LPPL 1.3c. Copyright (c) 2021 Jean-François Burnol
-% 2021/11/04 v0.99d
+% 2021/11/10 v1.0
 \edef\texdimensendinput{\endlinechar\the\endlinechar%
 \catcode`\noexpand _=\the\catcode`\_%
 \catcode`\noexpand @=\the\catcode`\@\relax\noexpand\endinput}%
 \endlinechar13\relax%
-\catcode`\_=11 \catcode`\@=11 % only for using \p@ of Plain. Check exists?
+\catcode`\_=11 \catcode`\@=11 % only for using \p@ (also \z@ now) of Plain. Check exists?
 %
-% Mathematics ("down" and "up" macros)
+% Mathematics (the challenge of the "down" and "up" macros)
 % ===========
 %
 % In the entire discussion here, "uu" stands for some core unit,
-% or some unit corresponding to a dimension > 1pt. For the case
-% of a unit corresponding to a dimension < 1pt, i.e. to
-% \texdimenwithunit macro added at 0.99, refer to the
-% comments of issue #2 on the tracker site.
+% or some unit corresponding to an internal dimension > 1pt.
 %
 % Is T sp attainable from unit "uu"?.
 % If not, what is largest dimension < Tsp which is attainable?
 % Here we suppose T>0.
 %
-% phi>1, psi=1/phi, psi<1.
+% D uu --> N = round(D * 65536) --> T = trunc (N * phi)
 %
-%     U(N,phi)=trunc(N phi) is the strictly increasing sequence,
-%     indexed by non-negative integers, of attainable dimensions.
-%     (in sp unit)
+% phi>1 is the conversion factor associated to "uu"
+% psi=1/phi, psi<1. Define U(N, phi) = trunc (N * phi)
+% 
+%     U(N,phi) is thus the strictly increasing sequence,
+%     indexed by non-negative integers, of non-negative
+%     attainable dimensions. (in sp unit)
+%
+% T>0, then:
 %
 %     U(N)<= T <  U(N+1)    iff    N = ceil((T+1)psi) - 1
 %     U(M)<  T <= U(M+1)    iff    M = ceil(T psi)    - 1
 %
-% Stumbling block
-% ---------------
+% In other words:
+%
+% - the largest attainable dimension not exceeding T sp
+%   is obtained via the integer "Zd = ceil((T+1)psi) - 1 = N",
+%   (i.e. find D with Zd=round(65536 D) then "D uu" is "down" approximation)
+%
+% - the smallest attainable dimension at least equal to T sp
+%   is obtained from the integer "Zu = ceil(T psi) = M + 1"
+%
+% - the two "Z"'s are either equal (i.e. T is attained) or Zu=Zd+1.
+%
+% The round((T+0.5)*psi) based approach (basis of the "\texdimenUU" macros)
+% -------------------------------------
+%
+% Recall in all of this T > 0. And phi>1, psi=1/phi<1.
+%
+% Let's return to our analysis of the
+%
+%         U(N)<= T < U(N+1) and U(M)< T <= U(M+1)
+%
+% equations. We will also use the N=Zd, and M+1=Zu notations.
+%
+% case1:  M = N, i.e. Zd<Zu, i.e. T is not attainable:
+%         M=N=Zd < T psi < (T+1) psi <= N+1=Zu
+%
+%         Then clearly R = round((T+0.5)psi) is either Zd or Zu.
+%         We will not know which one before computing trunc(R phi)
+%         and check if it is < T or > T.
+%
+%         As will be explained later trunc(R phi) can be computed very
+%         easily by hijacking TeX's handling of dimensions, we don't
+%         have to launch into \numexpr evaluations for that.
+%
+% case2:  M = N - 1, i.e. T = Zd = Zu is attained:
+%         T psi <= N < (T+1) psi, T = trunc(N phi)
+%
+%         Let v=(T+0.5)psi. As v = T psi + 0.5 psi it is < N+0.5
+%         And as v = (T+1)psi - 0.5psi it is > N - 0.5.
+%         So R = round(v) = N.
+%
+% We thus have the initial observation which was at the core of this
+% package initial release:
+%
+% - compute R = round((T+0.5) psi)
+%
+% - if T is attained, then T = trunc(R * phi)
+%
+% - if T is not attained then either { Zd = R and Zu = R+1 }  or 
+%    {Zd = R-1 and Zu = R}.
+%
+%   How do we check if R = Zd or Zu? We need to evaluate trunc(R phi) and
+%   compare it with T. This trunc(R phi) can be computed the following way:
 %
-% The stumbling block is that computing "ceil((T+1)psi) - 1" without
-% overflow is not obvious: yes \numexpr/\dimexpr allow so-called
-% "scaling operations" but only in the "rounding up" variant.
+%   - obtain D pt from \the\dimexpr R sp. Knuth's algorithm guarantees
+%   that R = round(D * 65536)
 %
-% If we attempt computing the ceil(x) function via round(x+0.5),
-% for example with psi=100/7227 which corresponds to the unit "in",
-% this necessitates evaluating:
+%   - then D uu where uu is the unit with conversion factor phi is
+%     converted by TeX into "trunc(R phi) sp", i.e.  trunc(R phi) =
+%     \number\dimexpr Duu\relax, where D pt = \the\dimexpr Rsp\relax.
 %
-%     round((((T+1)*200)+7227)/14454)
+% Conclusion:
 %
-% But as far as I can tell currently, for this we need to be able
-% to evaluate without overflow (T+1)*200+7227 and this limits to
-% T's which are (roughly) such that 100 T is less than \maxdimen.
+% 1. the macro \texdimenuu does the one-liner R=round((T+0.5) psi)
+%    then \the\dimexpr Rsp\relax and strips the "pt" unit
 %
-% A work-around
-% -------------
+% 2. macros \texdimenuuup and \texdimenuudown go further and check
+%    which one of Zd or Zu is R, obtaining thus Zd or Zu.
+%    [update: this describes the macros prior to 1.0 release]
 %
-% The rest of the discussion is about an algorithm providing an
-% alternative route to N, using \numexpr/\dimexpr/TeX facilities,
-% and working with (almost, as we will see) the full range of allowed
-% T's, 0 < T <= \maxdimen. (that the algorithm works for T=0 is to be
-% checked manually after the main discussion).
+% For units with conversion factor phi>2, a simplification is possible.
+% In that case let X = round(T psi) (it has the advantage compared to 
+% R that we can apply the formula without checking the sign of T).
 %
-% Let's return to the U(N)<= T < U(N+1) and U(M)< T <= U(M+1) equations.
+% Going back to our earlier analyis, now with psi < 0.5 (1uu>2pt)
 %
-% Either (recall in all of this T > 0):
+% case1: T is not attainable
+%        M=N=Zd < T psi < (T+1) psi <= N+1=Zu 
+%        As Zd < T psi < Zu, we have round(T psi) = Zd or Zu
 %
-% case1:  M = N, i.e. T is not attainable, M=N < T psi < (T+1) psi <= N+1
-% case2:  M = N - 1, i.e. T is attained, T psi <= N < (T+1) psi, T = trunc(N phi)
+% case2: T is attained, i.e. T psi <= N < (T+1) psi.
+%        As psi<0.5, and T psi + psi > N, we have T psi > N - 0.5.
+%        And T psi <= N so N = round(T psi).
 %
-% Let X = round(T psi). And let Y = trunc(X phi). We will explain later
-% how X and Y can be computed using \numexpr/\dimexpr/TeX.
+% So, for psi < 0.5, the X=round(T psi) can play the same role as
+% R=round((T+0.5)psi). If T is attained, we get the decimal D from this
+% X and if T is not attained we know that X is either Zd or Zu.
 %
-% case1: X can be N or N+1. It will be N+1 iff Y > T.
-% case2: X can be N or N-1. It will be N iff trunc((X+1)phi)>T.
+% The computations of X and Y=trunc(X phi) can be done independently of
+% sign of T.  But the final test has to be changed to Y < T if T < 0 and
+% then one must replace X by X+1. So we must filter out the sign of the
+% input.
 %
-% This is not convenient: if Y < T it could be that we are in case 2
-% but to decide we must check if trunc((X+1) phi) = T or not, so
-% this means a second computation.
+% Going back to the 1<phi<2 case, psi>0.5, then it would be slightly
+% less costly to compute X = round(T psi) than R = round((T + 0.5) psi),
+% but if we then realize that trunc(X phi) < T we do not yet know if
+% trunc((X+1) phi) = T or is > T, i.e. we don't know if Zd =X or X+1,
+% and we can not tell yet if T is attained or not.
 %
-% If psi < 0.5
-% ------------
+% In contrast if we find out that trunc(R phi)<T, we then know for sure
+% that Zd=R, Zu=R+1 and that T is not attained.
 %
-% The situation then simplifies:
+% Problems with \maxdimen in the obtention of Zu and Zd
+% -----------------------------------------------------
 %
-% case1: X can be N or N+1. It will be N+1 iff Y = trunc(X phi) > T.
-% case2: X is necessarily N.
+% Obtaining R = round((T+0.5)psi) has no risk of overflow.
+% But checking as described above which one of Zd or Zu (or both)
+% is R goes via a test computation which will cause overflow
+% if by bad luck R = Zu and Zu will give rise to a decimal D
+% such that D uu > \maxdimen.
 %
-% Thus:
-% a) compute X = round(T psi)
-% b) compute Y = trunc(X phi) and test if Y > T. If true, we
-%    were in case 1, replace X by X - 1, else we were either
-%    in case 1 or case 2, and we leave X as it is.
-% We have thus found N.
+% For T=\maxdimen (or very close) this is what happens for the units
+% "dd", "nc", and "in".
 %
-% The operation Y = trunc(X phi) can be achieved this way:
-% i) use \the\dimexpr to convert X sp into D pt,
-% ii) use \the\numexpr\dimexpr  to convert "D uu" into sp.
-% These steps give Y.
+% Besides, it turns out that this test which is done to decide whether
+% R=Zu or R=Zd, and on which the initial implementation of the macros
+% "up" and "down" was done at 0.9 gamma release is a bit costly.
 %
-% This way we find the maximal dimension at most T sp exactly
-% representable in "uu" unit.
+% At 1.0 release, all the "up" and "down" macros were re-implemented
+% via a more stubborn usage of the ceil() based formulae for Zd and Zu.
+% This made all usable even with \maxdimen input and besided, proved
+% on average slightly faster.
 %
-% The computations of X and Y can be done independently of sign of T.
-% But the final test has to be changed to Y < T if T < 0 and then
-% one must replace X by X+1. So we must filter out the sign of the input.
+% Overcoming the ceil() stumbling block (i.e. the 1.0 manner for "up/down")
+% ------------------------------------- 
 %
-% If the goal is only to find a decimal D such that "D uu" is
-% exactly T sp in the case this is possible, then things are simpler
-% because from X = round(T psi) we get D such as X sp is same as D pt
-% and "D uu" will work.
-% We don't have to take sign into account for this computation.
-% But if T sp was not attainable we don't know if this X will give
-% a D such that D uu < T sp or D uu > T sp.
+% I will in what follows refer to trunc(), floor() or ceil() only for
+% positive arguments, obtained as ratios x/y or sometimes as a numexpr
+% "scaling" operation" x*y/z which uses temporarily use doubled
+% precision.
 %
-% If psi > 0.5
-% ------------
+% As \numexpr's x/y is round(x/y), with rounding away from zero, we have
+% access to floor(t) for t>=0 as round(t+0.5)-1 and for t>0 also as
+% round(t-0.5). The former may cause overflow as it involves
+% (2x+y)/(2y) but the latter (2x-y)/(2y) will not overflow if x comes
+% from a dimension as 2x<2**31 then.
 %
-% For example unit "bp" has phi=803/800.
+% ceil(t) is more complex as it is floor(t)+1 only for t not an integer.
+% Let's explain how to overcome the challenge for Zd and the "in" unit,
+% i.e. a conversion factor of 7227/100.
 %
-% It is then not true that if T sp is attainable, the X = round(T psi)
-% will always work.
+% We want Zd = ceil((T+1)*100/7227) - 1, with T assumed positive.
 %
-% But it is true that R = round((T + 0.5) psi) will always work.
-% Here we must use -0.5 if T < 0, though.
+% Let T = k*7227 + r with 0<= r < 7227, 0<=k, and r>0 if k=0.
 %
-% This R=round((T+0.5) psi) can always be computed via \numexpr because 2T+1
-% will not trigger arithmetic overflow.
+% (T+1)*100/7227 becomes 100*k + (r+1)*100/7227 and thus
 %
-% So this gives an approach to find a D such that "D uu" is exactly
-% T sp when this is possible.
+% Zd = 100 * k + ceil(x) - 1
 %
-% If Tsp (positive) is not attainable, this R however can produce
-% either N or N+1.
+% with  x = n*100/7227, and n = 1+r, so 0<n<=7227
 %
-% But we can decide what happened by computing Z = trunc(R phi).
-% If and only if Z > T this means R was N+1.
+% Here we have a nice situation 0 < x <= 100. Then:
+% 
+% ceil(x) = 100 - floor(100 - x)
+%         = 100 - (round(100 - x + 0.5) - 1)
+%         = 101 - round(100 * (1 - n/7227) + 0.5)
+%         = 101 - round((200 * (7227 - n) + 7227)/14454)
 %
-% It is slightly less costly to compute X = round(T psi) than
-% R = round((T + 0.5) psi),
-% but if we then realize that trunc(X phi) < T  we do not yet know
-% if trunc((X+1) phi) = T  or is > T. So we proceed via R, not X,
-% to not have to make a second computation if a dimension comparison
-% test goes awry.
+% We can thus achieve the computation of Zd = ceil((T+1)*100/7227) - 1
+% for T>0 without overflow in \numexpr this way:
 %
-% To recapitulate: we have our algorithm for all units to find out
-% maximal dimension exactly attainable in "uu" unit and at most equal
-% to (positive) T sp.
+%     k = floor(T/7227) = round(T/7227 - 0.5)
+%                       = round((2*T - 7227) / 14454)  (T>0 used here)
 %
-% Unfortunately the check that Y (in case psi < 0.5) or Z (in case psi >
-% 0.5) verifies or not Y > T may trigger a Dimension too large error if
-% T sp was near non-attainable \maxdimen. It turns out this sad
-% situation happens only for the units `dd`, `nc`, and `in`, and T sp
-% very close to \maxdimen (like for all units apart from `pt`, `bp`,
-% `nd`, the \maxdimen is not attainable, and by bad luck for `dd`, `nc`,
-% and `in`, the X will correspond to a decimal D such that Duu>\maxdimen
-% is the nearest virtually attaible dimensions from above not from
-% below; see the README.md for the tabulation of the maximal usable inputs).
-%
-% Regarding the \texdimen<uu> macros, and units with phi > 2, I
-% hesitated using either the round((T+0.5)psi) or round(T psi), but for
-% Tsp = \maxdimen, both formulas turned out to give the same result for
-% all such units, so I chose for these \texdimen<uu> macros and the
-% units with phi>2 to use the simpler round(T psi) which does not need
-% to check the sign of T.
-%
-% For the "up" and "down" macros, we again use the round(T psi), but do
-% have to check the sign anyhow. We could also have used the
-% round((T+0.5)psi) which requires a sign check too, but it costs a bit
-% more. It would have allowed though to share the same codebase for all
-% units, here we have to prepare some slightly different shared macros
-% for the first batch bp, nd, dd and the second batch mm, pc, nc, cc,
-% cm, in.
+%     r = T - 7227 * k  = T modulo 7227
+%
+%     Zd = 100 * k + 100 - round( (201*7227 - 200*(r+1))/14454 )
+%
+% Everything here is computable within \numexpr and has absolutely no
+% potential overflow problem at all. The same analysis can be done for
+% Zu = ceil(T*100/7227) and for all core TeX units. See the comments
+% below for all obtained formulae and some additional details.
 %
 % Implementation
 % ==============
@@ -173,45 +210,6 @@
 {\catcode`p 12\catcode`t 12
  \csname expandafter\endcsname\gdef\csname texdimenstrippt\endcsname#1pt{#1}}%
 %
-% down macros:
-% for units with phi < 2:
-\def\texdimendown_A#1{\if-#1\texdimendown_neg\fi\texdimendown_B#1}%
-\def\texdimendown_B#1;#2;{\expandafter\texdimendown_c\the\numexpr(2*#1+1)#2;#1;}%
-% for units with phi > 2:
-\def\texdimendown_a#1{\if-#1\texdimendown_neg\fi\texdimendown_b#1}%
-\def\texdimendown_b#1;#2;{\expandafter\texdimendown_c\the\numexpr#1#2;#1;}%
-% shared macros:
-\def\texdimendown_c#1;{\expandafter\texdimendown_d\the\dimexpr#1sp;#1;}%
-{\catcode`P 12\catcode`T 12\lowercase{\gdef\texdimendown_d#1PT};#2;#3;#4;%
-   {\ifdim#1#4>#3sp \texdimendown_e{#2}\fi\texdimenfirstofone{#1}}%
-}%
-% this #2 will be \fi
-\def\texdimendown_e#1#2#3#4{#2\expandafter\texdimenstrippt\the\dimexpr\numexpr#1-1sp\relax}%
-% negative branch:
-% The problem here is that if input very small, output can be 0.0, and we
-% do not want -0.0 as output.
-% So let's do this somewhat brutally and non-efficiently.
-% Anyhow, negative inputs are not our priority.
-% #1 is \fi here and #2 is \texdimendown_b or _B:
-\def\texdimendown_neg#1#2-#3;#4;#5;{#1\expandafter\texdimenstrippt\the\dimexpr-#2#3;#4;#5;pt\relax}%
-%
-% up macros:
-\def\texdimenup_A#1{\if-#1\texdimenup_neg\fi\texdimenup_B#1}%
-\def\texdimenup_B#1;#2;{\expandafter\texdimenup_c\the\numexpr(2*#1+1)#2;#1;}%
-\def\texdimenup_a#1{\if-#1\texdimenup_neg\fi\texdimenup_b#1}%
-\def\texdimenup_b#1;#2;{\expandafter\texdimenup_c\the\numexpr#1#2;#1;}%
-\def\texdimenup_c#1;{\expandafter\texdimenup_d\the\dimexpr#1sp;#1;}%
-{\catcode`P 12\catcode`T 12\lowercase{\gdef\texdimenup_d#1PT};#2;#3;#4;%
-   {\ifdim#1#4<#3sp \texdimenup_e{#2}\fi\texdimenfirstofone{#1}}%
-}%
-% this #2 will be \fi
-\def\texdimenup_e#1#2#3#4{#2\expandafter\texdimenstrippt\the\dimexpr\numexpr#1+1sp\relax}%
-% negative branch:
-% Here we can me more expeditive than for the "down" macros.
-% But this breaks f-expandability.
-% #1 will be \fi and #2 is \texdimenup_b or _B:
-\def\texdimenup_neg#1#2-{#1-#2}%
-%
 % pt
 %
 \def\texdimenpt#1{\expandafter\texdimenstrippt\the\dimexpr#1\relax}%
@@ -222,10 +220,6 @@
 \def\texdimenbp_#1#2;{%
     \expandafter\texdimenstrippt\the\dimexpr\numexpr(2*#1#2+\if-#1-\fi1)*400/803sp\relax
 }%
-% \texdimenbpdown: maximal dim exactly expressible in bp and at most equal to input
-\def\texdimenbpdown#1{\expandafter\texdimendown_A\the\numexpr\dimexpr#1;*400/803;bp;}%
-% \texdimenbpup: minimal dim exactly expressible in bp and at least equal to input
-\def\texdimenbpup#1{\expandafter\texdimenup_A\the\numexpr\dimexpr#1;*400/803;bp;}%
 %
 % nd 685/642
 %
@@ -233,10 +227,6 @@
 \def\texdimennd_#1#2;{%
     \expandafter\texdimenstrippt\the\dimexpr\numexpr(2*#1#2+\if-#1-\fi1)*321/685sp\relax
 }%
-% \texdimennddown: maximal dim exactly expressible in nd and at most equal to input
-\def\texdimennddown#1{\expandafter\texdimendown_A\the\numexpr\dimexpr#1;*321/685;nd;}%
-% \texdimenndup: minimal dim exactly expressible in nd and at least equal to input
-\def\texdimenndup#1{\expandafter\texdimenup_A\the\numexpr\dimexpr#1;*321/685;nd;}%
 %
 % dd 1238/1157
 %
@@ -244,60 +234,296 @@
 \def\texdimendd_#1#2;{%
     \expandafter\texdimenstrippt\the\dimexpr\numexpr(2*#1#2+\if-#1-\fi1)*1157/2476sp\relax
 }%
-% \texdimendddown: maximal dim exactly expressible in dd and at most equal to input
-\def\texdimendddown#1{\expandafter\texdimendown_A\the\numexpr\dimexpr#1;*1157/2476;dd;}%
-% \texdimenddup: minimal dim exactly expressible in dd and at least equal to input
-\def\texdimenddup#1{\expandafter\texdimenup_A\the\numexpr\dimexpr#1;*1157/2476;dd;}%
 %
-% mm 7227/2540 phi now >2, use from here on the simpler approach
+% mm 7227/2540 phi now >2, use from here on the X = round(T psi) approach
 %
 \def\texdimenmm#1{\expandafter\texdimenstrippt\the\dimexpr(#1)*2540/7227\relax}%
-% \texdimenmmdown: maximal dim exactly expressible in mm and at most equal to input
-\def\texdimenmmdown#1{\expandafter\texdimendown_a\the\numexpr\dimexpr#1;*2540/7227;mm;}%
-% \texdimenmmup: minimal dim exactly expressible in mm and at least equal to input
-\def\texdimenmmup#1{\expandafter\texdimenup_a\the\numexpr\dimexpr#1;*2540/7227;mm;}%
 %
 % pc 12/1
 %
 \def\texdimenpc#1{\expandafter\texdimenstrippt\the\dimexpr(#1)/12\relax}%
-% \texdimenpcdown: maximal dim exactly expressible in pc and at most equal to input
-\def\texdimenpcdown#1{\expandafter\texdimendown_a\the\numexpr\dimexpr#1;/12;pc;}%
-% \texdimenpcup: minimal dim exactly expressible in pc and at least equal to input
-\def\texdimenpcup#1{\expandafter\texdimenup_a\the\numexpr\dimexpr#1;/12;pc;}%
 %
 % nc 1370/107
 %
 \def\texdimennc#1{\expandafter\texdimenstrippt\the\dimexpr(#1)*107/1370\relax}%
-% \texdimenncdown: maximal dim exactly expressible in nc and at most equal to input
-\def\texdimenncdown#1{\expandafter\texdimendown_a\the\numexpr\dimexpr#1;*107/1370;nc;}%
-% \texdimenncup: minimal dim exactly expressible in nc and at least equal to input
-\def\texdimenncup#1{\expandafter\texdimenup_a\the\numexpr\dimexpr#1;*107/1370;nc;}%
 %
 % cc 14856/1157
 %
 \def\texdimencc#1{\expandafter\texdimenstrippt\the\dimexpr(#1)*1157/14856\relax}%
-% \texdimenccdown: maximal dim exactly expressible in cc and at most equal to input
-\def\texdimenccdown#1{\expandafter\texdimendown_a\the\numexpr\dimexpr#1;*1157/14856;cc;}%
-% \texdimenccup: minimal dim exactly expressible in cc and at least equal to input
-\def\texdimenccup#1{\expandafter\texdimenup_a\the\numexpr\dimexpr#1;*1157/14856;cc;}%
 %
 % cm 7227/254
 %
 \def\texdimencm#1{\expandafter\texdimenstrippt\the\dimexpr(#1)*254/7227\relax}%
-% \texdimencmdown: maximal dim exactly expressible in cm and at most equal to input
-\def\texdimencmdown#1{\expandafter\texdimendown_a\the\numexpr\dimexpr#1;*254/7227;cm;}%
-% \texdimencmup: minimal dim exactly expressible in cm and at least equal to input
-\def\texdimencmup#1{\expandafter\texdimenup_a\the\numexpr\dimexpr#1;*254/7227;cm;}%
 %
 % in 7227/100
 %
 \def\texdimenin#1{\expandafter\texdimenstrippt\the\dimexpr(#1)*100/7227\relax}%
-% \texdimenindown: maximal dim exactly expressible in in and at most equal to input
-\def\texdimenindown#1{\expandafter\texdimendown_a\the\numexpr\dimexpr#1;*100/7227;in;}%
-% \texdimeninup: minimal dim exactly expressible in in and at least equal to input
-\def\texdimeninup#1{\expandafter\texdimenup_a\the\numexpr\dimexpr#1;*100/7227;in;}%
-% both in and cm
-% Mathematics ("both" macros)
+%
+% "no test" ceil()-based approach to up and down macros (release 1.0)
+% =====================================================
+%
+% The notation <u/v> means u/v in numexpr, which does rounding
+% away from zero. It is essential that the argument be >-0.5 else <x+1>
+% not same as <x>+1. All formulae are overflow free.
+%
+% The comments are for T > 0.
+%
+% Roughly such an approach works for phi = a/b > 1, such that:
+%
+%     a*(2b+1)<2**31 if a is odd, <2**32 if a is even
+%
+% This is true for all core units with quite some margin, the one with
+% largest a*b being phi=7227/2540 for "mm".
+%
+% Note: for a unit such as "ex" or "em" where morally b=65536=2**16,
+% this limits to a<=16383 if a is odd and to a<=32766 if a is even.
+% Thus the general \texdimenwithunit{dim1}{dim2} (which for dim2<1pt
+% computes basically an "up" value) can *not imitate fully* this scheme.
+%
+% The macros and formulas in the comments were obtained from a template
+% (see file generateupdownmacros.py at the project repository),
+% and we could actually combine them into a generic macro handling
+% general a/b (assuming above bounds are verified).
+% But for the the sake of efficiency, this is "rolled-out" here unit per unit.
+%
+\def\texdimenuudownup_zero#1;{\z@\relax}%
+\def\texdimenuudownup_neg#1-{-#1}%
+% bp 803/800
+% T = 803 k + r
+% Zd = 800 k + 800 - <(1284003 - 1600 r)/1606>
+\def\texdimenbpdown#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimenbpdown_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimenbpdown_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                        \krof \texdimenbpdown_b#1}%
+\def\texdimenbpdown_b#1;{\expandafter\texdimenbpdown_c\the\numexpr(2*#1-803)/1606;#1;}%
+\def\texdimenbpdown_c#1;#2;{\expandafter\texdimenbpdown_d\the\numexpr#2-803*#1;#1;}%
+\def\texdimenbpdown_d#1;#2;{\numexpr800*#2+800-(1284003-1600*#1)/1606sp\relax}%
+% Zu = 800 k + 800 + 1 - <(1285603 - 1600 r)/1606>
+\def\texdimenbpup#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimenbpup_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimenbpup_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                      \krof \texdimenbpup_b#1}%
+\def\texdimenbpup_b#1;{\expandafter\texdimenbpup_c\the\numexpr(2*#1-803)/1606;#1;}%
+\def\texdimenbpup_c#1;#2;{\expandafter\texdimenbpup_d\the\numexpr#2-803*#1;#1;}%
+\def\texdimenbpup_d#1;#2;{\numexpr800*#2+801-(1285603-1600*#1)/1606sp\relax}%
+% nd 685/642
+% T = 685 k + r
+% Zd = 642 k + 642 - <(878941 - 1284 r)/1370>
+\def\texdimennddown#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimennddown_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimennddown_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                        \krof \texdimennddown_b#1}%
+\def\texdimennddown_b#1;{\expandafter\texdimennddown_c\the\numexpr(2*#1-685)/1370;#1;}%
+\def\texdimennddown_c#1;#2;{\expandafter\texdimennddown_d\the\numexpr#2-685*#1;#1;}%
+\def\texdimennddown_d#1;#2;{\numexpr642*#2+642-(878941-1284*#1)/1370sp\relax}%
+% Zu = 642 k + 642 + 1 - <(880225 - 1284 r)/1370>
+\def\texdimenndup#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimenndup_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimenndup_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                      \krof \texdimenndup_b#1}%
+\def\texdimenndup_b#1;{\expandafter\texdimenndup_c\the\numexpr(2*#1-685)/1370;#1;}%
+\def\texdimenndup_c#1;#2;{\expandafter\texdimenndup_d\the\numexpr#2-685*#1;#1;}%
+\def\texdimenndup_d#1;#2;{\numexpr642*#2+643-(880225-1284*#1)/1370sp\relax}%
+% dd 1238/1157
+% T = 1238 k + r
+% Zd = 1157 k + 1157 - <(1431828 - 1157 r)/1238>
+\def\texdimendddown#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimendddown_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimendddown_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                        \krof \texdimendddown_b#1}%
+\def\texdimendddown_b#1;{\expandafter\texdimendddown_c\the\numexpr(#1-619)/1238;#1;}%
+\def\texdimendddown_c#1;#2;{\expandafter\texdimendddown_d\the\numexpr#2-1238*#1;#1;}%
+\def\texdimendddown_d#1;#2;{\numexpr1157*#2+1157-(1431828-1157*#1)/1238sp\relax}%
+% Zu = 1157 k + 1157 + 1 - <(1432985 - 1157 r)/1238>
+\def\texdimenddup#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimenddup_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimenddup_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                      \krof \texdimenddup_b#1}%
+\def\texdimenddup_b#1;{\expandafter\texdimenddup_c\the\numexpr(#1-619)/1238;#1;}%
+\def\texdimenddup_c#1;#2;{\expandafter\texdimenddup_d\the\numexpr#2-1238*#1;#1;}%
+\def\texdimenddup_d#1;#2;{\numexpr1157*#2+1158-(1432985-1157*#1)/1238sp\relax}%
+% mm 7227/2540
+% T = 7227 k + r
+% Zd = 2540 k + 2540 - <(36715307 - 5080 r)/14454>
+\def\texdimenmmdown#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimenmmdown_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimenmmdown_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                        \krof \texdimenmmdown_b#1}%
+\def\texdimenmmdown_b#1;{\expandafter\texdimenmmdown_c\the\numexpr(2*#1-7227)/14454;#1;}%
+\def\texdimenmmdown_c#1;#2;{\expandafter\texdimenmmdown_d\the\numexpr#2-7227*#1;#1;}%
+\def\texdimenmmdown_d#1;#2;{\numexpr2540*#2+2540-(36715307-5080*#1)/14454sp\relax}%
+% Zu = 2540 k + 2540 + 1 - <(36720387 - 5080 r)/14454>
+\def\texdimenmmup#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimenmmup_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimenmmup_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                      \krof \texdimenmmup_b#1}%
+\def\texdimenmmup_b#1;{\expandafter\texdimenmmup_c\the\numexpr(2*#1-7227)/14454;#1;}%
+\def\texdimenmmup_c#1;#2;{\expandafter\texdimenmmup_d\the\numexpr#2-7227*#1;#1;}%
+\def\texdimenmmup_d#1;#2;{\numexpr2540*#2+2541-(36720387-5080*#1)/14454sp\relax}%
+% pc 12/1
+% T = 12 k + r
+% Zd = 1 k + 1 - <(17 - 1 r)/12>
+\def\texdimenpcdown#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimenpcdown_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimenpcdown_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                        \krof \texdimenpcdown_b#1}%
+\def\texdimenpcdown_b#1;{\expandafter\texdimenpcdown_c\the\numexpr(#1-6)/12;#1;}%
+\def\texdimenpcdown_c#1;#2;{\expandafter\texdimenpcdown_d\the\numexpr#2-12*#1;#1;}%
+\def\texdimenpcdown_d#1;#2;{\numexpr#2+1-(17-#1)/12sp\relax}%
+% Zu = 1 k + 1 + 1 - <(18 - 1 r)/12>
+\def\texdimenpcup#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimenpcup_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimenpcup_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                      \krof \texdimenpcup_b#1}%
+\def\texdimenpcup_b#1;{\expandafter\texdimenpcup_c\the\numexpr(#1-6)/12;#1;}%
+\def\texdimenpcup_c#1;#2;{\expandafter\texdimenpcup_d\the\numexpr#2-12*#1;#1;}%
+\def\texdimenpcup_d#1;#2;{\numexpr#2+2-(18-#1)/12sp\relax}%
+% nc 1370/107
+% T = 1370 k + r
+% Zd = 107 k + 107 - <(147168 - 107 r)/1370>
+\def\texdimenncdown#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimenncdown_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimenncdown_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                        \krof \texdimenncdown_b#1}%
+\def\texdimenncdown_b#1;{\expandafter\texdimenncdown_c\the\numexpr(#1-685)/1370;#1;}%
+\def\texdimenncdown_c#1;#2;{\expandafter\texdimenncdown_d\the\numexpr#2-1370*#1;#1;}%
+\def\texdimenncdown_d#1;#2;{\numexpr107*#2+107-(147168-107*#1)/1370sp\relax}%
+% Zu = 107 k + 107 + 1 - <(147275 - 107 r)/1370>
+\def\texdimenncup#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimenncup_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimenncup_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                      \krof \texdimenncup_b#1}%
+\def\texdimenncup_b#1;{\expandafter\texdimenncup_c\the\numexpr(#1-685)/1370;#1;}%
+\def\texdimenncup_c#1;#2;{\expandafter\texdimenncup_d\the\numexpr#2-1370*#1;#1;}%
+\def\texdimenncup_d#1;#2;{\numexpr107*#2+108-(147275-107*#1)/1370sp\relax}%
+% cc 14856/1157
+% T = 14856 k + r
+% Zd = 1157 k + 1157 - <(17194663 - 1157 r)/14856>
+\def\texdimenccdown#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimenccdown_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimenccdown_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                        \krof \texdimenccdown_b#1}%
+\def\texdimenccdown_b#1;{\expandafter\texdimenccdown_c\the\numexpr(#1-7428)/14856;#1;}%
+\def\texdimenccdown_c#1;#2;{\expandafter\texdimenccdown_d\the\numexpr#2-14856*#1;#1;}%
+\def\texdimenccdown_d#1;#2;{\numexpr1157*#2+1157-(17194663-1157*#1)/14856sp\relax}%
+% Zu = 1157 k + 1157 + 1 - <(17195820 - 1157 r)/14856>
+\def\texdimenccup#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimenccup_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimenccup_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                      \krof \texdimenccup_b#1}%
+\def\texdimenccup_b#1;{\expandafter\texdimenccup_c\the\numexpr(#1-7428)/14856;#1;}%
+\def\texdimenccup_c#1;#2;{\expandafter\texdimenccup_d\the\numexpr#2-14856*#1;#1;}%
+\def\texdimenccup_d#1;#2;{\numexpr1157*#2+1158-(17195820-1157*#1)/14856sp\relax}%
+% cm 7227/254
+% T = 7227 k + r
+% Zd = 254 k + 254 - <(3678035 - 508 r)/14454>
+\def\texdimencmdown#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimencmdown_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimencmdown_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                        \krof \texdimencmdown_b#1}%
+\def\texdimencmdown_b#1;{\expandafter\texdimencmdown_c\the\numexpr(2*#1-7227)/14454;#1;}%
+\def\texdimencmdown_c#1;#2;{\expandafter\texdimencmdown_d\the\numexpr#2-7227*#1;#1;}%
+\def\texdimencmdown_d#1;#2;{\numexpr254*#2+254-(3678035-508*#1)/14454sp\relax}%
+% Zu = 254 k + 254 + 1 - <(3678543 - 508 r)/14454>
+\def\texdimencmup#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimencmup_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimencmup_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                      \krof \texdimencmup_b#1}%
+\def\texdimencmup_b#1;{\expandafter\texdimencmup_c\the\numexpr(2*#1-7227)/14454;#1;}%
+\def\texdimencmup_c#1;#2;{\expandafter\texdimencmup_d\the\numexpr#2-7227*#1;#1;}%
+\def\texdimencmup_d#1;#2;{\numexpr254*#2+255-(3678543-508*#1)/14454sp\relax}%
+% in 7227/100
+% T = 7227 k + r
+% Zd = 100 k + 100 - <(1452427 - 200 r)/14454>
+\def\texdimenindown#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimenindown_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimenindown_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                        \krof \texdimenindown_b#1}%
+\def\texdimenindown_b#1;{\expandafter\texdimenindown_c\the\numexpr(2*#1-7227)/14454;#1;}%
+\def\texdimenindown_c#1;#2;{\expandafter\texdimenindown_d\the\numexpr#2-7227*#1;#1;}%
+\def\texdimenindown_d#1;#2;{\numexpr#200+100-(1452427-2*#100)/14454sp\relax}%
+% Zu = 100 k + 100 + 1 - <(1452627 - 200 r)/14454>
+\def\texdimeninup#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimeninup_a\the\numexpr\dimexpr#1;%
+}%
+\def\texdimeninup_a#1{\texdimenzerominusfork
+                        #1-\texdimenuudownup_zero
+                        0#1\texdimenuudownup_neg
+                         0-{}%
+                      \krof \texdimeninup_b#1}%
+\def\texdimeninup_b#1;{\expandafter\texdimeninup_c\the\numexpr(2*#1-7227)/14454;#1;}%
+\def\texdimeninup_c#1;#2;{\expandafter\texdimeninup_d\the\numexpr#2-7227*#1;#1;}%
+\def\texdimeninup_d#1;#2;{\numexpr#200+101-(1452627-2*#100)/14454sp\relax}%
+%
+% Mathematics ("both in and cm", added at release 0.9)
 % ===========
 %
 % Let a and b be two non-negative integers such that U = floor(a 7227/100) = 
@@ -374,8 +600,110 @@
 \def\texdimenbothsp_neg#1#2-#3;{#1-\numexpr(((2*#3-3612)/7227)*7227+1)/2-1\relax\relax}%
 %
 \let\texdimenbothcminsp\texdimenbothincmsp
-% (\texdimenwithunit, added at 0.99)
-% Mathematics
+%
+% Mathematics and Algorithm ("both mm and bp", added at release 1.0)
+% =========================
+% We start from a dimension expressed in sp unit, "T sp". Assume T positive.
+% We know how to get largest "X sp <= T sp" which is exactly expressible
+% in mm unit
+% i.e. can be written X=trunc(a 7227/2540) for some non-negative integer a.
+% We want to achieve X=trunc(b 803/800) for some b.
+%
+%    Only the congruence of X modulo 803 matters for this.
+%    It turns out that the mod 803 impossible values are 267, 535, 802.
+%    As pointed out by Ruixi Zhang on the package repo issue #10,
+%    when a<--a+2540, X increases by 7227=9*803 hence the value
+%    modulo 803 does not change. Thus only "a modulo 2540" matters
+%    to check if X(a) is attainable with bp unit. Ruixi Zhang found by
+%    brute force that there are modulo 2540 nine excluded a-values
+%
+% Rather than checking if "a mod. 2540" avoids the 9 Ruixi Zhang values
+% or if "X mod. 803" avoids  267, 535, 802, we will simply basically
+% check if X sp = \texdimenbp{X sp}bp, as this approach is probably
+% about the same cost or even less than computing "X mod. 803" and
+% correspondingly branching.
+%
+% The key is that if "a" is bad, then "a-1" is automatically good as
+% pointed out by R.Z. on #10, which can be seen without knowing the 9
+% bad congruences, simply by noticing that a<--a-1 modifies X either to
+% X-2 or X-3, so if X was bad certainly the new one is not.
+%
+% Once "a" has gotten its final value, we apply "\the\dimexpr a sp
+% = D pt" trick to recover the D such that "D mm" gives rise to the found
+% dimension.  We go via this "Dmm" intermediary also to express the final
+% result as "X sp", because anyhow the "X" we worked with and had in
+% our token stream has to be recomputed if a<--a-1, so lets always
+% recompute it from final "a", and this goes via "D mm" (but see
+% the paragraph MEMO for alternative for this trunc(a 7227/2540) step).
+%
+% I will copy here the style I used for bothincm expansion triggering
+% via an already positioned \dimexpr waiting to output final result.
+\def\texdimenbothbpmm#1{\expandafter\texdimenstrippt\the\dimexpr
+                        \expandafter\texdimenbothbpmm_fork\the\numexpr\dimexpr#1;}%
+% the \texdimenzerominusfork is defined in the \texdimenwithunit section
+\def\texdimenbothbpmm_fork#1{\texdimenzerominusfork
+                             #1-\texdimenbothbpmm_zero
+                             0#1\texdimenbothbpmm_neg
+                             0-\texdimenbothbpmm_a
+                             \krof#1}%
+% because this is *inside* a pre-positioned \dimexpr, we don't have
+% to worry about zero output ending up as -0.0
+\def\texdimenbothbpmm_neg-{-\texdimenbothbpmm_a}%
+\def\texdimenbothbpmm_zero#1;{\z@\relax}%
+% now, find X sp <= T sp maximal and expressible in mm unit
+% it will be X=trunc(a 7227/2540), we first get a candidate for "a"
+\def\texdimenbothbpmm_a#1;%
+    {\expandafter\texdimenbothbpmm_b\the\numexpr#1*2540/7227;#1;}%
+% we get in a single line the X from this candidate, hijacking TeX's
+% built-in *7227/2540... the "MEMO" above explains one could do this
+% purely within \numexpr, working around its division rounds, and
+% avoiding overflow, but I suspect this would be more costly.
+\def\texdimenbothbpmm_b#1;{\expandafter\texdimenbothbpmm_c
+    \the\numexpr\dimexpr\expandafter\texdimenstrippt\the\dimexpr#1spmm;#1;}%
+% now we have X;a;T;
+\def\texdimenbothbpmm_c#1;#2;#3;{%
+% If X>T, our candidate "a=#2" must be decreased by 1 and we go to _ca
+% The original #3 is not needed anymore
+    \ifnum#1>#3 \expandafter\texdimenbothbpmm_ca\fi
+% Else we decide whether it is "a" or "a-1" we must use. I preferred
+% to induce a re-grabbing cost here, rather than have \texdimenbothbpmm_ca
+% re-grab its arguments from \texdimenbothbpmm_d replacement text.
+    \texdimenbothbpmm_d#1;#2;%
+}%
+% Here, dynamically at the time of the concluding \dimexpr, we
+% check if X sp is expressible in bp unit and then use "a" or "a-1"
+% accordingly
+\def\texdimenbothbpmm_d#1;#2;{#2sp%
+    \ifnum\dimexpr
+    \expandafter\texdimenstrippt\the\dimexpr\numexpr(2*#1+1)*400/803spbp=#1
+    \else-1sp\fi
+% and a \relax to stop the concluding \dimexpr
+    \relax
+}%
+% Here we must decrease "a=#2" by 1, recompute X=#1, then loop
+% back to \texdimenbothbpmm_d. Hesitation between forcing a
+% re-grab or doing it in one step with the subtraction of 1 done twice
+\def\texdimenbothbpmm_ca\texdimenbothbpmm_d#1;#2;%
+   {\expandafter\texdimenbothbpmm_cb\the\numexpr#2-1;}%
+\def\texdimenbothbpmm_cb#1;{%
+    \expandafter\texdimenbothbpmm_d
+    \the\numexpr\dimexpr\expandafter\texdimenstrippt\the\dimexpr#1spmm;#1;%
+}%
+% done...
+% now the lazy way for \texdimenbothmmbp
+\def\texdimenbothmmbp#1{\expandafter\texdimenstrippt\the\dimexpr
+    \expandafter\texdimenbothmmbp_a\the\numexpr\dimexpr\texdimenbothbpmm{#1}mm;}%
+% or remove the + and do \if-#1-\else+\fi1 ?
+% If zero at this stage, we will correctly get 0.0 in the end
+\def\texdimenbothmmbp_a#1#2;{\numexpr(2*#1#2+\if-#1-\fi1)*400/803sp\relax}%
+% \texdimenbothbpmmpt and its alias \texdimenbothmmbppt
+\def\texdimenbothbpmmpt#1{\texdimenpt{\texdimenbothbpmm{#1}mm}}%
+\let\texdimenbothmmbppt\texdimenbothbpmmpt
+% \texdimenbothbpmmsp and its alias \texdimenbothmmbpsp
+\def\texdimenbothbpmmsp#1{\the\numexpr\dimexpr\texdimenbothbpmm{#1}mm\relax\relax}%
+\let\texdimenbothmmbpsp\texdimenbothbpmmsp
+%
+% Mathematics (\texdimenwithunit, added at release 0.99, modified at 1.0)
 % ===========
 %
 % As explained in the README.md, the ex and em units are
@@ -404,96 +732,177 @@
 % If f/65536<1, all dimensions Tsp are attainable as D uu. Indeed
 % D uu is parsed by TeX via N=round(D*65536), then T=trunc(N*phi),
 % with phi=f/65536. Starting from T we need to find an N such that
-% T/phi <= N< (T+1)/phi. We first consider v=(T+0.5)/phi. As its
-% distance to the extremities is 0.5/phi>0.5,(*) its rounding M
+% T/phi <= N< (T+1)/phi. 
+%
+% This is equivalent to ceil(T/phi)<= N < ceil((T+1)/phi)
+%
+% Now obsolete remark: let v=(T+0.5)/phi. As its
+% distance to the extremities is 0.5/phi>0.5, (phi>1) its rounding M
 % to an integer verifies automatically T/phi < M < (T+1)/phi, so
-% is a candidate. The TeX core conversion of Msp to a Dpt with
-% D a decimal of at most 5 fractional digits will provide a D
-% such that indeed M=round(D*65536).
+% is a valid candidate. This was used at 0.99 release.
+% (it is funny that N=round((2T+1)*32768/f) works for all f>0
+%  *except* f=65536).
 %
-% (*) attention that this fails if phi=1, we definitely can not
-% use the N=round((2T+1)*32768/f) formula for f=65536, it gives
-% N=T+1 in place of N=T... It is funny the formula works for all
-% f apart from f=65536...
+% The 1.0 release chooses to implement the ceil(T/phi) formula rather as
+% it is closer to naive expectation "dim1/dim2" of a division.
 %
-% Attention that if T=0, the M=round(0.5/phi) will be at least 1,
-% so the produced decimal D will not be 0.0, although it will
-% be true that D \dimexpr f sp\relax produces the zero dimension,
-% the above analysis being correct even for T=0.
+% It is not obvious to compute this ceil(T/phi) without overflow.
 %
 % Implementation
 % ==============
 %
-% \texdimenwithunit{dim1}{dim2}. dim2>0 assumed.
+% \texdimenwithunit{dim1}{dim2}
+%
+% First done at 0.99, then refactored at 1.0:
+% - to add support for dim2<0pt
+% - to handle differently the dim2<1pt case and make the output
+%   closer to mathematical dim1/dim2
+%
+% To handle dim2<0pt, we simply simultaneously do
+% dim1<-- (-dim1) and dim2<-- (-dim2).
+%
+% dim2=0pt is not intercepted and will cause division by zero low-level
+% error.  Code comments below were not adjusted and handle only
+% dim2>0pt.
+%
 % We first get f from dim2 and branch according to whether f>65536,
-% f=65536 or f<65536.
+% or f<=65536.
 % We will also need to check the sign of T (dim1=T sp).
-% We then compute in both branches round((2T+1)*32768/f), but
-% in an indirect way in the f<65536 branch to avoid overflow.
-% We will need also to intercept T=0 to avoid producing a non
-% pleasing positive D in that case (still correct for f<65536;
-% this is issue #3).
-%
-% In the f<65536 branch we first do the Euclidean division
-% 2T+1 = k*2*f + R with 0<=R<2f. The k is obtained as round((2T+1-f)/(2f))
-% which can be computed in a numexpr (and never gives -1 even for T=0)
-% Then (2T+1)*32768/f=65536*k + R*32768/f
-% Then R*32768/f<=65536-32768/f<65536-32768/65536=65536-0.5
-% Hence the numexpr evaluation gives an integer B<65536.
-%
-% N.B.: if f>=65536, we still have R*32768/f<65536 as R<2f
-% so the only difference is that the B could be here 65536
-%
-% From \the\dimexpr Bsp, we get a decimal E < 1, so E=0.ddddd
-% (or less digits) and the looked for D will be the contatenation
-% k.ddddd with k as obtained earlier. This procedure has no possible
-% arithmetic overflow.
-%
-% #2 is assumed positive. The case #2=1pt=65536sp is special, and
-% must be filtered out (issues #4, #6).
-% pre-multiplication of #1 by 2 here
-\def\texdimenwithunit#1#2{\expandafter\texdimenwithunit_
-    \the\numexpr\dimexpr#2\expandafter;\the\numexpr2*\dimexpr#1;}%
+% f>65536: we compute round((2T+1)*32768/f)
+% f=65536: merged with f<65536 branch
+% f<65536: we compute ceil(T*65536/f)
+%
+%          rationale: round((2T+1)*32768/f) which was used at 0.99
+%          would be ok [if f<65536 not f=65536 ! cf #3, #4] 
+%          also for this branch
+%
+%          BUT it diverges noticeably from naive expectation
+%          dim1/dim2 "=" T*65536/f the more so when f is small.
+%          See issue #16 and also the discussion at #13.
+%
+%          As was explained in comments to issue #2 which asked for a
+%          \texdimenwithunit the ceil(T*65536/f) is the smallest
+%          allowable choice
+%
+%          To avoid arithmetic overflow issues we first do the euclidean
+%          division T = k f + r, 0<= r < f
+%
+%          The final result in sp units will be k*65536 + ceil(r * 65536/f)
+%          We don't do this k*65536 explicitly as it may overflow
+%          but output the decimal k + E where E is the conversion
+%          to a decimal 0.ddddd or 1.0 of "ceil(r * 65536/f) sp"
+%          (which is at most 65536sp=1pt so E is at most 1.0).
+%
+%          To compute the Euclidean quotient k in \numexpr we use there
+%          (2T-f)/(2f) i.e. round((2T-f)/2f) = trunc(T/f)
+%          as we are careful to never have T=0 inthere...
+%
+%          Computing C = ceil(r * 65536/f) in \numexpr is the delicate
+%          part, as r can be as large as f-1 hence 65535 and r*65536 would
+%          overflow. We could compute R=round(r*65536/f) ("scaling operation")
+%          then C=R+1 if R*f-65536*r<0, else C=R.
+%
+%          The problem is then: how to get the sign of R*f-65536*r without
+%          overflow? I considered various ways.
+%
+%          But in the end, opted for simply this:
+%          - get R=round(r*65536/f) as \the\numexpr r*65536/f ("scaling" no overflow)
+%          - hence get E pt=\the\dimexpr R sp
+%          - let TeX compute E<f sp>: if it turns out < r sp,
+%            then C=R+1,
+%            else C=R. Done.
+%
+\def\texdimenwithunit#1#2{\expandafter\texdimenwithunit_i
+% no premultiplication of dim1 by 2 as was done for technical
+% reasons when dim2<1pt branch used round((2T+1)*32768/f)
+    \the\numexpr\dimexpr#2\expandafter;\the\numexpr\dimexpr#1;}%
+%
+% so tempted to do \input xintkernel.sty to have some utilities...
+% not even a \@gobble in Plain...
 \let\texdimenorthat\texdimenfirstofone
 \def\texdimendothis#1#2\texdimenorthat#3{\fi#1}%
-\def\texdimenwithunit_#1;#2{%
-        \ifnum#1=\p@\texdimendothis\texdimenwithunit_p@\fi
+\def\texdimengobtominus#1-{}%
+\def\texdimenwithunit_i#1{%
+     \texdimengobtominus#1\texdimenwithunit_switchsigns-%
+     \texdimenwithunit_j#1%
+}%
+\def\texdimenwithunit_switchsigns-\texdimenwithunit_j-#1;#2%
+{%
+% due to \texdimenwithunit_Bneg we can not simply prefix dim1
+% with -, as -0 is bad there. So let's check also if #2 is 0
+    \texdimenzerominusfork
+      #2-\texdimenwithunit_Bzero % abusive double usage
+      0#2\texdimenwithunit_j     % abusive shortcut
+       0-{\texdimenwithunit_ic#2}%
+    \krof
+    #1;%
+}%
+\def\texdimenwithunit_ic#1#2;{\texdimenwithunit_j#2;-#1}%
+\def\texdimenwithunit_j#1;#2{%
+        % \ifnum#1=\p@\texdimendothis\texdimenwithunit_p@\fi
         \ifnum#1>\p@\texdimendothis\texdimenwithunit_A\fi
         \texdimenorthat\texdimenwithunit_B#2#1;%
 }%
-% we premultiplied the first argument by 2... must be undone now!
-% and we must avoid overflow. Very underoptimal, but user is not
-% supposed to do something as silly as actually using this unit=1pt case!
-\def\texdimenwithunit_p@#1#2;#3;{\expandafter\texdimenstrippt\the\dimexpr\numexpr#1#3/2sp\relax}%
-% unit>1pt, handle this as for bp. Attention it would be wrong for unit=1pt!
-\def\texdimenwithunit_A#1#2;#3;{\expandafter\texdimenstrippt
-    \the\dimexpr\numexpr(#1#3+\if-#1-\fi1)*32768/#2sp\relax
+% not needed, as the "ceil" branch is fine to use for f = 65536
+% \def\texdimenwithunit_p@#1#2;#3;{%
+%     \expandafter\texdimenstrippt\the\dimexpr#1#3sp/2\relax
+% }%
+% unit>1pt, handle this as for bp. Attention it would be wrong for
+% unit=1pt!
+\def\texdimenwithunit_A#1#2;#3;{%
+    \expandafter\texdimenstrippt
+    \the\dimexpr\numexpr(2*#1#3+\if-#1-\fi1)*32768/#2sp\relax
+    % - fine if dim1>0, <0, or =0
+    % - with *\p@ better but an early doubled dim2 would complicate 1pt
+    % test and not sure if doing \p@/(2*#2) here advantageous
 }%
-% unit<1pt. Attention it would be wrong for unit=1pt!
+% unit<=1pt.
 % if dim1<0, simply negate result for dim1>0 as it can not possibly be 0.0
-% Indeed (2T+1)*32768/f will be at least 3*32768/f so its rounding at least 2.
-% Add a check for dim1=0 to fix issue #3
+% Indeed T*65536/f will be at least 1 so its ceil also (in fact ceil
+% will even be at least 2 if f<65536).
+% The dim1=0 case must get filtered out due to way of calculating the
+% "ceil" in \numexpr
 \def\texdimenzerominusfork #10-#2#3\krof {#2}%
 \def\texdimenwithunit_B#1{\texdimenzerominusfork
                            #1-\texdimenwithunit_Bzero
                            0#1\texdimenwithunit_Bneg
                            0-\texdimenwithunit_Ba
                           \krof#1}%
-\def\texdimenwithunit_Ba#1#2;#3;{\expandafter\texdimenwithunit_Bb\the\numexpr#1#3+1;#2;}%
-\def\texdimenwithunit_Bb#1;#2;{\expandafter\texdimenwithunit_Bc\the\numexpr(#1-#2)/(2*#2);#1;#2;}%
-% I was hesitating between leaving k in the stream (breaking f-expandability)
-% and then remove the "0" and trailing "pt" from 0.ddddd pt, but opted
-% finally for doing \the\numexpr..+0.ddddd which is f-expandable and allows
-% recycling strippt here with no need of another utility.
-%
-% This means (see the nota bene above) that we could apply this procedure
-% also for f>65536, because at worst we will get a \the\numexprk+1.0, which
-% gives the correct result. I tested and found about 39% longer execution time
-% if dim2>1pt does same calculations as for dim2<1pt, and at the same time
-% dim2<1pt was improved about 11% from skipping the conditional
-\def\texdimenwithunit_Bc#1;#2;#3;{\the\numexpr#1+\expandafter\texdimenstrippt
-                                  \the\dimexpr\numexpr(#2-#1*2*#3)*32768/#3sp\relax}%
-% Here, definitely not caring about f-expandability. Or efficiency.
+\def\texdimenwithunit_Ba#1#2;#3;{%
+    % no overflow possible from 2*#1#3=2*dim1
+    \expandafter\texdimenwithunit_Bb\the\numexpr(2*#1#3-#2)/(2*#2);#1#3;#2;%
+}%
+% now k;T;f;. Get the remainder r=T-k*f
+\def\texdimenwithunit_Bb#1;#2;#3;{%
+    \expandafter\texdimenwithunit_Bc\the\numexpr#2-#1*#3;#1;#3;%
+}%
+% now r;k;f;. We can start \the\numexpr k+ ....
+% and there we will need to get R=round(r*65536/f), Ept=Rsp,
+% check if E"f sp"<"r sp", if yes replace R by R+1 else keep E etc.
+\def\texdimenwithunit_Bc#1;#2;#3;{%
+% \the\numexpr k+0.ddddd is handy because it can well be actually
+% \the\numexpr k+1.0, now that we use ceil approach in this branch
+    \the\numexpr#2+\expandafter\texdimenwithunit_Bd
+                   \the\numexpr #1*\p@/#3;#1;#3;%
+}%
+% R;r;f;
+\def\texdimenwithunit_Bd#1;{%
+    \expandafter\texdimenwithunit_Be\the\dimexpr#1sp;#1;%
+}%
+% Ept;R;r;f;
+% #1=0.ddd... or 1.0 but has no end marker hence the
+% \texdimenfirstofone{#1} as in \texdimendown_d and \texdimenup_d
+{\catcode`P 12\catcode`T 12
+\lowercase{\gdef\texdimenwithunit_Be#1PT};#2;#3;#4;{%
+    \ifdim#1\dimexpr#4sp<#3sp \texdimenwithunit_Bf{#2}\fi
+    \texdimenfirstofone{#1}%
+    }%
+}%
+% #2 is \fi. Add a dimen storage \onesp for 1sp?
+\def\texdimenwithunit_Bf#1#2\texdimenfirstofone#3{#2%
+    \expandafter\texdimenstrippt\the\dimexpr#1sp+1sp\relax
+}%
+% Here definitely not caring about f-expandability. Or efficiency.
 \def\texdimenwithunit_Bneg-{-\texdimenwithunit_Ba{}}%
 \def\texdimenwithunit_Bzero#1;#2;{0.0}%
 \texdimensendinput