2 * Copyright (C) 1999,2000 Erik Walthinsen <omega@cse.ogi.edu>
3 * 2000 Wim Taymans <wtay@chello.be>
4 * 2004 Wim Taymans <wim@fluendo.com>
5 * 2015 Jan Schmidt <jan@centricular.com>
7 * gstclock-linreg.c: Linear regression implementation, used in clock slaving
9 * This library is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Library General Public
11 * License as published by the Free Software Foundation; either
12 * version 2 of the License, or (at your option) any later version.
14 * This library is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Library General Public License for more details.
19 * You should have received a copy of the GNU Library General Public
20 * License along with this library; if not, write to the
21 * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
22 * Boston, MA 02110-1301, USA.
25 #include "gst_private.h"
31 #include "glib-compat-private.h"
33 /* Compute log2 of the passed 64-bit number by finding the highest set bit */
35 gst_log2 (GstClockTime in)
38 { 0x2, 0xC, 0xF0, 0xFF00, 0xFFFF0000, 0xFFFFFFFF00000000LL };
39 const guint64 S[] = { 1, 2, 4, 8, 16, 32 };
43 for (i = 5; i >= 0; i--) {
53 /* http://mathworld.wolfram.com/LeastSquaresFitting.html
57 _priv_gst_do_linear_regression (GstClockTime * times, guint n,
58 GstClockTime * m_num, GstClockTime * m_denom, GstClockTime * b,
59 GstClockTime * xbase, gdouble * r_squared)
61 GstClockTime *newx, *newy;
62 GstClockTime xmin, ymin, xbar, ybar, xbar4, ybar4;
63 GstClockTime xmax, ymax;
64 GstClockTimeDiff sxx, sxy, syy;
70 xbar = ybar = sxx = syy = sxy = 0;
75 xmin = ymin = G_MAXUINT64;
77 for (i = j = 0; i < n; i++, j += 4) {
78 xmin = MIN (xmin, x[j]);
79 ymin = MIN (ymin, y[j]);
81 xmax = MAX (xmax, x[j]);
82 ymax = MAX (ymax, y[j]);
88 /* strip off unnecessary bits of precision */
89 for (i = j = 0; i < n; i++, j += 4) {
90 newx[j] = x[j] - xmin;
91 newy[j] = y[j] - ymin;
94 #ifdef DEBUGGING_ENABLED
95 GST_CAT_DEBUG (GST_CAT_CLOCK, "reduced numbers:");
96 for (i = j = 0; i < n; i++, j += 4)
97 GST_CAT_DEBUG (GST_CAT_CLOCK,
98 " %" G_GUINT64_FORMAT " %" G_GUINT64_FORMAT, newx[j], newy[j]);
101 /* have to do this precisely otherwise the results are pretty much useless.
102 * should guarantee that none of these accumulators can overflow */
104 /* quantities on the order of 1e10 to 1e13 -> 30-35 bits;
105 * window size a max of 2^10, so
106 this addition could end up around 2^45 or so -- ample headroom */
107 for (i = j = 0; i < n; i++, j += 4) {
108 /* Just in case assumptions about headroom prove false, let's check */
109 if ((newx[j] > 0 && G_MAXUINT64 - xbar <= newx[j]) ||
110 (newy[j] > 0 && G_MAXUINT64 - ybar <= newy[j])) {
111 GST_CAT_WARNING (GST_CAT_CLOCK,
112 "Regression overflowed in clock slaving! xbar %"
113 G_GUINT64_FORMAT " newx[j] %" G_GUINT64_FORMAT " ybar %"
114 G_GUINT64_FORMAT " newy[j] %" G_GUINT64_FORMAT, xbar, newx[j], ybar,
125 /* multiplying directly would give quantities on the order of 1e20-1e26 ->
126 * 60 bits to 70 bits times the window size that's 80 which is too much.
127 * Instead we (1) subtract off the xbar*ybar in the loop instead of after,
128 * to avoid accumulation; (2) shift off some estimated number of bits from
129 * each multiplicand to limit the expected ceiling. For strange
130 * distributions of input values, things can still overflow, in which
131 * case we drop precision and retry - at most a few times, in practice rarely
134 /* Guess how many bits we might need for the usual distribution of input,
135 * with a fallback loop that drops precision if things go pear-shaped */
136 max_bits = gst_log2 (MAX (xmax - xmin, ymax - ymin)) * 7 / 8 + gst_log2 (n);
138 pshift = max_bits - 64;
142 #ifdef DEBUGGING_ENABLED
143 GST_CAT_DEBUG (GST_CAT_CLOCK,
144 "Restarting regression with precision shift %u", pshift);
147 xbar4 = xbar >> pshift;
148 ybar4 = ybar >> pshift;
150 for (i = j = 0; i < n; i++, j += 4) {
151 GstClockTime newx4, newy4;
152 GstClockTimeDiff tmp;
154 newx4 = newx[j] >> pshift;
155 newy4 = newy[j] >> pshift;
157 tmp = (newx4 + xbar4) * (newx4 - xbar4);
158 if (G_UNLIKELY (tmp > 0 && sxx > 0 && (G_MAXINT64 - sxx <= tmp))) {
160 /* Drop some precision and restart */
164 } while (G_MAXINT64 - sxx <= tmp);
166 } else if (G_UNLIKELY (tmp < 0 && sxx < 0 && (G_MAXINT64 - sxx >= tmp))) {
168 /* Drop some precision and restart */
172 } while (G_MININT64 - sxx >= tmp);
177 tmp = newy4 * newy4 - ybar4 * ybar4;
178 if (G_UNLIKELY (tmp > 0 && syy > 0 && (G_MAXINT64 - syy <= tmp))) {
183 } while (G_MAXINT64 - syy <= tmp);
185 } else if (G_UNLIKELY (tmp < 0 && syy < 0 && (G_MAXINT64 - syy >= tmp))) {
190 } while (G_MININT64 - syy >= tmp);
195 tmp = newx4 * newy4 - xbar4 * ybar4;
196 if (G_UNLIKELY (tmp > 0 && sxy > 0 && (G_MAXINT64 - sxy <= tmp))) {
201 } while (G_MAXINT64 - sxy <= tmp);
203 } else if (G_UNLIKELY (tmp < 0 && sxy < 0 && (G_MININT64 - sxy >= tmp))) {
208 } while (G_MININT64 - sxy >= tmp);
215 if (G_UNLIKELY (sxx == 0))
221 *b = (ybar + ymin) - gst_util_uint64_scale (xbar, *m_num, *m_denom);
222 *r_squared = ((double) sxy * (double) sxy) / ((double) sxx * (double) syy);
224 #ifdef DEBUGGING_ENABLED
225 GST_CAT_DEBUG (GST_CAT_CLOCK, " m = %g", ((double) *m_num) / *m_denom);
226 GST_CAT_DEBUG (GST_CAT_CLOCK, " b = %" G_GUINT64_FORMAT, *b);
227 GST_CAT_DEBUG (GST_CAT_CLOCK, " xbase = %" G_GUINT64_FORMAT, *xbase);
228 GST_CAT_DEBUG (GST_CAT_CLOCK, " r2 = %g", *r_squared);
235 GST_CAT_DEBUG (GST_CAT_CLOCK, "sxx == 0, regression failed");