From 289f92da7e1a3c4433143dce4765d046e53c0a55 Mon Sep 17 00:00:00 2001 From: Ichthyostega Date: Thu, 1 Dec 2022 23:23:50 +0100 Subject: [PATCH] Timeline: safely calculate sum/difference of large fractional times MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ...in a similar vein as done for the product calculation. In this case, we need to check the dimensions carefully and pick the best calculation path, but as long as the overall result can be represented, it should be possible to carry out the calculation with fractional values, albeit introducing a small error. As a follow-up, I have now also refactored the re-quantisation functions, to be usable for general requantisation to another grid, and I used these to replace the *naive* implementation of the conversion FSecs -> µ-Grid, which caused a lot of integer-wrap-around However, while the test now works basically without glitch or wrap, the window position is still numerically of by 1e-6, which becomes quite noticeably here due to the large overall span used for the test. --- src/lib/rational.hpp | 44 ++++-- src/lib/time/time.cpp | 6 +- src/lib/util-quant.hpp | 7 + src/stage/model/zoom-window.hpp | 64 +++++++- tests/library/rational-test.cpp | 5 +- tests/stage/model/zoom-window-test.cpp | 11 +- wiki/thinkPad.ichthyo.mm | 210 +++++++++++++++++++++++++ 7 files changed, 326 insertions(+), 21 deletions(-) diff --git a/src/lib/rational.hpp b/src/lib/rational.hpp index 03b6d2a12..f1bd84851 100644 --- a/src/lib/rational.hpp +++ b/src/lib/rational.hpp @@ -73,6 +73,7 @@ #include #include +#include "lib/util-quant.hpp" namespace util { @@ -154,6 +155,39 @@ namespace util { } + /** + * Re-Quantise a number into a new grid, truncating to the next lower grid point. + * @remark Grid-aligned values can be interpreted as rational numbers (integer fractions), + * where the quantiser corresponds to the denominator and the numerator counts + * the number of grid steps. To work both around precision problems and the + * danger of integer wrap-around, the integer division is performed on the + * old value and then the re-quantisation done on the remainder, using + * 128bit floating point for maximum precision. This operation can + * also be used to re-form a fraction to be cast in terms of the + * new quantiser; this introduces a tiny error, but typically + * allows for safe or simplified calculations. + * @param num the count in old grid-steps (#den) or the numerator + * @param den the old quantiser or the denominator of a fraction + * @param u the new quantiser or the new denominator to use + * @return the adjusted numerator, so that the fraction with u + * will be almost the same than dividing #num by #den + */ + inline int64_t + reQuant (int64_t num, int64_t den, int64_t u) + { + u = 0!=u? u:1; + auto [d,r] = util::iDiv (num, den); + using f128 = long double; + // round to smallest integer fraction, to shake off "number dust" + f128 const ROUND_ULP = 1 + 1/(f128(std::numeric_limits::max()) * 2); + + // construct approximation quantised to 1/u + f128 frac = f128(r) / den; + int64_t res = d*u + int64_t(frac*u * ROUND_ULP); + ENSURE (abs (f128(res)/u - rational_cast(Rat{num,den})) <= 1.0/abs(u)); + return res; + } + /** * re-Quantise a rational number to a (typically smaller) denominator. * @param u the new denominator to use @@ -170,15 +204,7 @@ namespace util { inline Rat reQuant (Rat src, int64_t u) { - int64_t d = rational_cast (src); - int64_t r = src.numerator() % src.denominator(); - using f128 = long double; - - // construct approximation quantised to 1/u - f128 frac = rational_cast (Rat{r, src.denominator()}); - Rat res = d + int64_t(frac*u) / Rat(u); - ENSURE (abs (rational_cast(src) - rational_cast(res)) <= 1.0/abs(u)); - return res; + return Rat{reQuant (src.numerator(), src.denominator(), u), u}; } } // namespace util diff --git a/src/lib/time/time.cpp b/src/lib/time/time.cpp index c25ba27c7..cb38b3a2a 100644 --- a/src/lib/time/time.cpp +++ b/src/lib/time/time.cpp @@ -46,6 +46,7 @@ #include "lib/error.hpp" #include "lib/time.h" #include "lib/time/timevalue.hpp" +#include "lib/rational.hpp" #include "lib/util-quant.hpp" #include "lib/format-string.hpp" @@ -336,7 +337,10 @@ lumiera_tmpbuf_print_time (gavl_time_t time) gavl_time_t lumiera_rational_to_time (FSecs const& fractionalSeconds) { - return rational_cast (fractionalSeconds * lib::time::TimeValue::SCALE); + return gavl_time_t(util::reQuant (fractionalSeconds.numerator() + ,fractionalSeconds.denominator() + ,lib::time::TimeValue::SCALE + )); } gavl_time_t diff --git a/src/lib/util-quant.hpp b/src/lib/util-quant.hpp index bca65cd68..1f5d94cd0 100644 --- a/src/lib/util-quant.hpp +++ b/src/lib/util-quant.hpp @@ -78,6 +78,13 @@ namespace util { { } }; + template + inline IDiv + iDiv (I num, I den) ///< support type inference and auto typing... + { + return IDiv{num,den}; + } + /** floor function for integer arithmetics. * Unlike the built-in integer division, this function diff --git a/src/stage/model/zoom-window.hpp b/src/stage/model/zoom-window.hpp index 5d59ff262..e08f57baa 100644 --- a/src/stage/model/zoom-window.hpp +++ b/src/stage/model/zoom-window.hpp @@ -116,6 +116,7 @@ namespace model { using util::min; using util::max; + using util::sgn; namespace { ///////////////////////////////////////////////////////////////////////////////////////////////TICKET #1259 : reorganise raw time base datatypes : need conversion path into FSecs /** @@ -143,6 +144,12 @@ namespace model { % duration.denominator(); } + inline double + approx (Rat r) + { + return util::rational_cast (r); + } + /////////////////////////////////////////////////////////////////////////////////////////////////////////TICKET #1261 : why the hell did I define time entities to be immutable. Looks like a "functional programming" fad in hindsight /** @todo we need these only because the blurry distinction between * lib::time::TimeValue and lib::time::Time, which in turn is caused @@ -521,12 +528,13 @@ namespace model { static FSecs scaleSafe (FSecs duration, Rat factor) { - auto approx = [](Rat r){ return rational_cast (r); }; - if (not util::can_represent_Product(duration, factor)) { - if (approx(MAX_TIMESPAN) < approx(duration) * approx (factor)) - return MAX_TIMESPAN; // exceeds limits of time representation => cap the result + auto guess{approx(duration) * approx (factor)}; + if (approx(MAX_TIMESPAN) < abs(guess)) + return MAX_TIMESPAN * sgn(guess); // exceeds limits of time representation => cap the result + if (0 == guess) + return 0; // slightly adjust the factor so that the time-base denominator cancels out, // allowing to calculate the product without dangerous multiplication of large numbers @@ -538,6 +546,52 @@ namespace model { return duration * factor; } + /** + * Calculate sum (or difference) of possibly large time durations, avoiding integer wrap-around. + * Again, this is a heuristics, based on re-quantisation to a smaller common denominator. + * @return exact result if representable, otherwise approximation + * @note result is capped to MAX_TIMESPAN when exceeding domain + */ + static FSecs + addSafe (FSecs t1, FSecs t2) + { + if (not util::can_represent_Sum (t1,t2)) + { + auto guess{approx(t1) + approx(t2)}; + if (approx(MAX_TIMESPAN) < abs(guess)) + return MAX_TIMESPAN * sgn(guess); // exceeds limits => cap the result + + // re-Quantise numbers to achieve a common denominator, + // thus avoiding to multiply numerators for normalisation + int64_t n1 = t1.numerator(); + int64_t d1 = t1.denominator(); + int s1 = sgn(n1)*sgn(d1); + n1 = abs(n1); d1 = abs(d1); + int64_t n2 = t2.numerator(); + int64_t d2 = t2.denominator(); + int s2 = sgn(n2)*sgn(d2); + n2 = abs(n2); d2 = abs(d2); + // quantise to smaller denominator to avoid increasing any numerator + int64_t u = d1 0 // check numerators to detect danger of wrap-around + and (62>= 1; // danger zone! wrap-around imminent + + n1 = d1==u? n1 : reQuant (n1,d1, u); + n2 = d2==u? n2 : reQuant (n2,d2, u); + FSecs res{s1*n1 + s2*n2, u}; + ENSURE (abs (guess - approx(res)) < 1.0/u); + return detox (res); + } + else + // directly calculate ordinary numbers... + return t1 + t2; + } + static Rat establishMetric (uint pxWidth, Time startWin, Time afterWin) @@ -765,7 +819,7 @@ namespace model { Rat posFactor = canvasOffset / FSecs{afterAll_-startAll_}; posFactor = parabolicAnchorRule (posFactor); // also limited 0...1 FSecs partBeforeAnchor = scaleSafe (duration, posFactor); - startWin_ = startAll_ + (canvasOffset - partBeforeAnchor); + startWin_ = startAll_ + addSafe (canvasOffset, -partBeforeAnchor); establishWindowDuration (duration); startAll_ = min (startAll_, startWin_); afterAll_ = max (afterAll_, afterWin_); diff --git a/tests/library/rational-test.cpp b/tests/library/rational-test.cpp index fe8f8e127..9a1460398 100644 --- a/tests/library/rational-test.cpp +++ b/tests/library/rational-test.cpp @@ -359,11 +359,12 @@ namespace test { CHECK (util::toString (sleazy+1) == "134217727/16777216sec"); // also works towards larger denominator, or with negative numbers... - CHECK (reQuant (poison, MAX-7) == 104811045873349724_r/14973006553335675); + CHECK (reQuant (1/poison, MAX) == 1317624576693539413_r/9223372036854775807); CHECK (reQuant (-poison, 7777) == -54438_r/ 7777); CHECK (reQuant (poison, -7777) == -54438_r/-7777); - CHECK (approx (reQuant (poison, MAX-7)) == 7); + CHECK (approx ( 1/poison ) == 0.142857149f); + CHECK (approx (reQuant (1/poison, MAX)) == 0.142857149f); CHECK (approx (reQuant (poison, 7777)) == 6.99987125f); } }; diff --git a/tests/stage/model/zoom-window-test.cpp b/tests/stage/model/zoom-window-test.cpp index 08323ad11..837128d5d 100644 --- a/tests/stage/model/zoom-window-test.cpp +++ b/tests/stage/model/zoom-window-test.cpp @@ -554,8 +554,10 @@ namespace test { CHECK (poison == 206435633551724850_r/307445734561825883); CHECK (2_r/3 < poison and poison < 1); // looks innocuous... CHECK (poison + Time::SCALE < 0); // simple calculations fail due to numeric overflow - CHECK (Time(FSecs(poison)) < Time::ZERO); // conversion to µ-ticks also leads to overflow - CHECK (-6 == _raw(Time(FSecs(poison)))); + CHECK (poison * Time::SCALE < 0); + CHECK (-6 == rational_cast(poison * Time::SCALE)); // naive conversion to µ-ticks would leads to overflow + CHECK (671453 == _raw(Time(FSecs(poison)))); // however the actual conversion routine is safeguarded + CHECK (671453.812f == rational_cast(poison)*Time::SCALE); using util::ilog2; CHECK (40 == ilog2 (LIM_HAZARD)); // LIM_HAZARD is based on MAX_INT / Time::Scale @@ -624,7 +626,7 @@ namespace test { Rat poisonousDuration = win.pxWidth() / poison; // Now, to demonstrate this "poison" was actually dangerous CHECK (poisonousDuration == 7071251894921995309_r/8257425342068994); // ...when we attempt to calculate the new duration directly.... - CHECK (Time(poisonousDuration) < Time::ZERO); // ...then a conversion to TimeValue will cause integer wrap + CHECK (poisonousDuration * Time::SCALE < 0); // ...then a conversion to TimeValue will cause integer wrap CHECK(856.350708f == rational_cast (poisonousDuration)); // yet numerically the duration actually established is almost the same CHECK(856.350708f == rational_cast (_FSecs(win.visible().duration()))); CHECK (win.px_per_sec() == 575000000_r/856350691); // the new metric however is comprised of sanitised fractional numbers @@ -642,7 +644,8 @@ namespace test { SHOW_EXPR(win.px_per_sec()); SHOW_EXPR(win.pxWidth()); SHOW_EXPR(_raw(win.overallSpan().duration()) * rational_cast (poison)) - Time targetPos{TimeValue{gavl_time_t(_raw(win.overallSpan().duration()) * rational_cast (poison))}}; + TimeValue targetPos{gavl_time_t(_raw(win.overallSpan().duration()) + * rational_cast (poison))}; SHOW_EXPR(targetPos); SHOW_EXPR(_raw(targetPos)); SHOW_EXPR(_raw(win.visible().start())) diff --git a/wiki/thinkPad.ichthyo.mm b/wiki/thinkPad.ichthyo.mm index cf66a25cd..d8b16744d 100644 --- a/wiki/thinkPad.ichthyo.mm +++ b/wiki/thinkPad.ichthyo.mm @@ -40386,6 +40386,122 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

+ man muß Floating-point runden wenn man glatte Werte will +

+ + + + + + + + + + + + +

+ das ist essentiell wichtig. "Negative" Zeiten dürfen sich keinesfalls  anders verhalten. Eine andere Quantsierungs-Regel kann man dann ggfs. high-level auf ein left-Truncate aufsetzen (z.B. Mitte Frame-Intervall) +

+ + + + + + + + + + +

+ ⟹ also muß man hier technisch runden +

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

+ ...das sind 4ms +

+ + + + + + + + + +

+ parsing '1/250sec' resulted in 0:00:00.003 instead of 0:00:00.004 +

+ + + + + + + @@ -40434,6 +40550,100 @@ + + + + + + + + + + + +

+ ...tut sie zwar nicht in dem Beispiel hier, aber mit genügend krimineller Energie ließe sich ein valides Beispiel konstruieren, wobei +

+
    +
  • + die Ziel-Position dann außerhalb des legalen Bereichs liegen würde +
    • +
  • + bei korrekter Behandlung daher das Ergebnis-Fenster in den legalen Bereich geschoben werden müßte +
    • +
  • + aber ohne weitere Schutzmaßname hier die Berechnung der Zielposition entgleist +
    • +
+ + + + + + + + + + + + + + + + + + + + +

+ 0111 + 0111 = 1110 +

+

+ 0011 + 0101 = 1000 +

+

+ aber... +

+

+ 0010 + 0101 = 0111 +

+ + + + + + + + + + + + + + + +

+ Gefahr besteht... +

+
    +
  • + wenn beide Vorzeichen gleichgerichtet sind +
    • +
  • + wenn mindestens einer der Zähler das 63te Bit gesetzt hat (2^62) +
    • +
+ + + + + + + + + + +