Problems with iterative process

Questions
,
1

For some reasons I need to consider the following iterative process.
We have three positive numbers with constant total such that x+y+z == t
Consider the iterative algorithm:

  • if the sum of numbers is not equal to t - stop (with message “Total error”)

  • if one of the numbers is t/2 - stop (with message “Half”)

  • if one of the numbers is greater than t/2 (let it be z) - returns new numbers 2z-t,2x,2y

  • if all numbers are less than t/2 - return the new numbers t-2z,t-2x,t-2y

      ClearAll[genNumsReals]; genNumsReals[ tot_] :=
   Block[{a = RandomReal[{0., tot}], b},
    b = RandomReal[{0., tot - a}];
    {a, b, tot - a - b}];
  
  ClearAll[genNumsRationals]; genNumsRationals[e_, tot_] :=
   Block[{a = Rationalize[RandomReal[{e, tot - e}], e], b},
    b = Rationalize[RandomReal[{e, tot - a - e}], e];
    {a, b, tot - a - b}];
  
  Clear[nextNums]; nextNums[nums_, tot_] :=
   Which[
    Unequal[Total@nums, tot], {nums, "Error Total"},
    MemberQ[nums, tot/2], {nums, "Half Total"},
    AnyTrue[nums, GreaterThan[tot/2]], Mod[2 nums, tot],
    True, tot - 2 nums];
  
  numsNest[init_, tot_, iterMax_ : 100] :=
    NestWhileList[nextNums[#, tot] &, init,
     (UnsameQ[##] && Length@Last[{##}] == 3) &, All, iterMax];

Firstly, is it equal to Mod[-2 {x,y,z}, t]? Looks like not.
If the numbers are rational, the results are generally expected but interesting.
Either cycles of greater or lesser length are obtained,
either t/2 occurs.

But for reals I see some problems. For any initial data with machine precision
we get “Half Total” sooner or later (here t=1):

{{0.84474, 0.117199, 0.0380613}, {0.68948, 0.234397, 
  0.0761226}, {0.37896, 0.468795, 0.152245}, {0.24208, 0.0624103, 
  0.69551}, {0.48416, 0.124821, 0.39102}, {0.0316802, 0.750359, 
  0.217961}, {0.0633604, 0.500718, 0.435922}, {0.126721, 0.00143543, 
  0.871844}, {0.253441, 0.00287087, 0.743688}, {0.506883, 0.00574174, 
  0.487375}, {0.0137658, 0.0114835, 0.974751}, {0.0275315, 0.0229669, 
  0.949502}, {0.0550631, 0.0459339, 0.899003}, {0.110126, 0.0918678, 
  0.798006}, {0.220252, 0.183736, 0.596012}, {0.440505, 0.367471, 
  0.192024}, {0.118991, 0.265058, 0.615952}, {0.237981, 0.530115, 
  0.231903}, {0.475963, 0.060231, 0.463806}, {0.0480747, 0.879538, 
  0.0723872}, {0.0961494, 0.759076, 0.144774}, {0.192299, 0.518152, 
  0.289549}, {0.384598, 0.0363047, 0.579098}, {0.769195, 0.0726095, 
  0.158195}, {0.538391, 0.145219, 0.31639}, {0.0767815, 0.290438, 
  0.632781}, {0.153563, 0.580876, 0.265561}, {0.307126, 0.161751, 
  0.531122}, {0.614252, 0.323503, 0.062245}, {0.228504, 0.647006, 
  0.12449}, {0.457009, 0.294012, 0.24898}, {0.0859828, 0.411977, 
  0.50204}, {0.171966, 0.823954, 0.00408077}, {0.343931, 0.647907, 
  0.00816154}, {0.687862, 0.295815, 0.0163231}, {0.375725, 0.591629, 
  0.0326462}, {0.75145, 0.183258, 0.0652924}, {0.502899, 0.366516, 
  0.130585}, {0.00579834, 0.733032, 0.261169}, {0.0115967, 0.466064, 
  0.522339}, {0.0231934, 0.932129, 0.0446777}, {0.0463867, 0.864258, 
  0.0893555}, {0.0927734, 0.728516, 0.178711}, {0.185547, 0.457031, 
  0.357422}, {0.628906, 0.0859375, 0.285156}, {0.257813, 0.171875, 
  0.570313}, {0.515625, 0.34375, 0.140625}, {0.03125, 0.6875, 
  0.28125}, {0.0625, 0.375, 0.5625}, {0.125, 0.75, 0.125}, {0.25, 0.5,
   0.25}, {{0.25, 0.5, 0.25}, "Half Total"}}

and for higher precision we get “Total error” (sooner or later too):

ClearAll[genNumsReals]; genNumsReals[tot_] :=
 Block[{a = RandomReal[{0., tot}, WorkingPrecision -> 30], b},
  b = RandomReal[{0., tot - a}, WorkingPrecision -> 30];
  {a, b, tot - a - b}];

Clear[nextNums]; nextNums[nums_, tot_] :=
 Which[
  Unequal[Total@nums, tot], {nums, "Error Total"},
  MemberQ[nums, tot/2], {nums, "Half Total"},
  AnyTrue[nums, GreaterThan[tot/2]],
  SetPrecision[Mod[2 nums, tot], 30],
  True, SetPrecision[tot - 2 nums, 30]];

{{0.659587099030496875713494880795, 0.185853135754329321586581016710,
  0.154559765215173802699924102495}, \
{0.319174198060993674630481109489, 0.371706271508658647739764546714,
  0.309119530430347622118603112540}, \
{0.361651603878012650739037781022, 0.256587456982682704520470906573,
  0.381760939139304755762793774920}, \
{0.276696792243974698521924437955, 0.486825086034634590959058186854,
  0.236478121721390488474412450159}, \
{0.446606415512050602956151124090, 0.0263498279307308180818836262915,
  0.527043756557219023051175099681}, \
{0.893212831024101205912302248180, 0.0526996558614616361637672525831,
  0.0540875131144380461023501993623}, \
{0.786425662048202411824604496360, 0.105399311722923272327534505166,
  0.108175026228876092204700398725}, \
{0.572851324096404823649208992720, 0.210798623445846544655069010332,
  0.216350052457752184409400797449}, \
{0.145702648192809647298417985439, 0.421597246891693089310138020664,
  0.432700104915504368818801594898}, \
{0.708594703614380705403164029121, 0.156805506216613821379723958671,
  0.134599790168991262362396810204}, \
{{0.708594703614380705403164029121, 0.156805506216613821379723958671,
   0.134599790168991262362396810204}, "Error Total"}}

I think the problems are in rounding, but how to understand it?

2

Floating-point analysis

\def\fl{\mathop{\text{fl}}}

For roundoff error in double precision:

Let \fl(u) represent the correctly rounded floating-point number closest to the real number u.
Let x,y,z be floating-point numbers between 0 and 1 (so that x = \fl(x) etc.).
Let x+y+z be true sum and let \fl(x+y+z) the computed sum.
If \fl(x,+y+z)=1 (compared with or without tolerance), then

x+y+z=1+\delta

where \delta is some small roundoff error bounded by

{-\tau_{-} \le \delta \le \tau_{+}}

with the tolerances \tau_{-} \ge \varepsilon/4, \tau_{+} \ge \varepsilon/2 and \varepsilon = 2^{-52} (machine epsilon).

Now 2u, is computed without roundoff error,
as are 2u-1 and 1-2u, if 0.5 \le u < 1.0.
It is not hard to see that x'+y'+z' will be either 1+2\delta or 1-2\delta (3rd or 4th rule resp.),
if (x',y',z') is the iterate following (x,y,z).
If \delta \ne 0, then eventually the iteration stops with “Total Error”, unless a “Half Total” incidentally occurs earlier.

If \delta=0, then the following determines the behavior.
If 0.5 \le u < 1 and the trailing n bits of the mantissa of u are zero,
then the trailing n+1 bits of the mantissa of 2u-1 and 1-2u are zero.
After at most 52 such transformations, u = 0.5 and you get a “Half Total”.
Since each of x, y, and z would experience 52 such transformations
if the iteration never stopped, a “Half Total” is inevitable.

Example: Case \delta \ne 0

seq = numsNest[{0.84474, 0.117199} // Append[#, 1 - Total[#]] &, 1.]
(*
{{0.84474, 0.117199, 0.038061}, {0.68948, 0.234398, 
  0.076122}, {0.37896, 0.468796, 0.152244}, {0.24208, 0.062408, 
  0.695512}, {0.48416, 0.124816, 0.391024}, {0.03168, 0.750368, 
  0.217952}, {0.06336, 0.500736, 0.435904}, {0.12672, 0.001472, 
  0.871808}, {0.25344, 0.002944, 
  0.743616}, {{0.25344, 0.002944, 0.743616}, "Error Total"}}
*)

The error in the total doubles at each step:

Total /@ Most[seq] - 1
Ratios[%]
(*
{0., -1.11022*10^-16, -2.22045*10^-16, 4.44089*10^-16, 
 8.88178*10^-16, -1.77636*10^-15, -3.55271*10^-15, -7.10543*10^-15, 
-1.42109*10^-14}

{ComplexInfinity, 2.`, -2.`, 2.`, -2.`, 2.`, 2.`, 2.`}
*)

The initial zero is because \fl(x+y+z)=1 exactly, but (x+y+z)-1 = -2^{-54}:

SetPrecision[{0.84474, 0.117199} // Append[#, 1 - Total[#]] &, Infinity] //
  Total // N[# - 1] &
(* -5.55112*10^-17  *)

The error after the first step is twice this, namely -2^{53} = -1.11\times10^{-16}.

Example: Case \delta = 0.

Remove the roundoff from an initial value less than 0.5 to get a total that equals 1:

SetPrecision[{0.84474, 0.117199 + 2^-54} // 
    Append[#, 1 - Total[#]] &, Infinity] //
  Total // N[# - 1] &
(*  0.  *)

seq2 = numsNest[{0.84474, 0.117199 + 2^-54} // Append[#, 1 - Total[#]] &, 1.]
(*
{{0.84474, 0.117199, 0.038061}, {0.68948, 0.234398, 0.076122}, 
  ... 
  {0.0136719, 0.244141, 0.742188}, {0.0273438, 0.488281, 0.484375}, 
  {0.945313, 0.0234375, 0.03125}, {0.890625, 0.046875, 0.0625}, 
  {0.78125, 0.09375, 0.125}, {0.5625, 0.1875, 0.25}, {0.125, 0.375, 0.5}, 
  {{0.125, 0.375, 0.5}, "Half Total"}}
*)

A visualization of the effective loss of precision:

MatrixPlot@RealDigits[#, 2][[All, 1]] & /@ 
  Transpose@seq2[[;; -2]] // GraphicsRow

7933Untitled.png
7933Untitled.png1935×601 13 KB

3

Thank you very much!
The problem is that by the results for rationals, there can be deterministic chaos.
And we need to understand that for reals can be natural behavior, and what - errors of rounding.
I see that the subject is not easy.

4

So SetPrecision[tot, Infinity] helps a lot