假设您有一个浮点数列表，这些浮点数大约是一个共同数量的倍数，例如
2.468、3.700、6.1699
大约是 1.234 的所有倍数。您如何描述这种“近似 gcd”，您将如何进行计算或估计？
与我对this question 的回答密切相关。
【问题讨论】：
标签： algorithm language-agnostic math floating-point greatest-common-divisor
您可以运行 Euclid 的 gcd 算法，将小于 0.01（或您选择的少量数字）的任何值作为伪 0。使用您的数字：
3.700 = 1 * 2.468 + 1.232,
2.468 = 2 * 1.232 + 0.004.
所以前两个数字的伪 gcd 是 1.232。现在你用你的最后一个数字来获取这个 gcd：
6.1699 = 5 * 1.232 + 0.0099.
所以 1.232 是伪 gcd，倍数是 2,3,5。为了改善这个结果，您可以对数据点进行线性回归：
(2,2.468), (3,3.7), (5,6.1699).
斜率是改进后的伪gcd。
警告：第一部分是算法在数值上不稳定——如果你从非常脏的数据开始，你就有麻烦了。
【讨论】：
将您的测量值表示为最低值的倍数。因此，您的列表变为 1.00000、1.49919、2.49996。这些值的小数部分将非常接近 1/Nths，因为 N 的某个值取决于您的最低值与基频的接近程度。我建议循环增加 N 直到找到足够精细的匹配。在这种情况下，对于 N=1（即假设 X=2.468 是您的基频），您会发现标准偏差为 0.3333（三个值中有两个是 X * 1 的 0.5），这是不可接受的高。对于 N=2（即假设 2.468/2 是您的基频），您会发现标准偏差几乎为零（所有三个值都在 X/2 的倍数的 0.001 以内），因此 2.468/2 是您的近似值GCD。
我的计划中的主要缺陷是，当最低测量值最准确时效果最好，但情况可能并非如此。这可以通过多次执行整个操作来缓解，每次丢弃测量列表中的最小值，然后使用每次通过的结果列表来确定更精确的结果。另一种改进结果的方法是调整 GCD 以最小化 GCD 的整数倍与测量值之间的标准偏差。
【讨论】：
这让我想起了寻找实数的良好有理数逼近的问题。标准技术是连分数展开式：
def rationalizations(x):
assert 0 <= x
ix = int(x)
yield ix, 1
if x == ix: return
for numer, denom in rationalizations(1.0/(x-ix)):
yield denom + ix * numer, numer
我们可以将此直接应用到 Jonathan Leffler 和 Sparr 的方法中：
>>> a, b, c = 2.468, 3.700, 6.1699
>>> b/a, c/a
(1.4991896272285252, 2.4999594813614263)
>>> list(itertools.islice(rationalizations(b/a), 3))
[(1, 1), (3, 2), (925, 617)]
>>> list(itertools.islice(rationalizations(c/a), 3))
[(2, 1), (5, 2), (30847, 12339)]
从每个序列中挑选出第一个足够好的近似值。 （这里是 3/2 和 5/2。）或者不是直接比较 3.0/2.0 和 1.499189...，您可能会注意到 925/617 使用比 3/2 大得多的整数，使得 3 /2 一个很好的停留地点。
除以哪个数字并不重要。 （例如，使用 a/b 和 c/b 可以得到 2/3 和 5/3。）一旦有了整数比率，就可以使用 shsmurfy 的线性回归来细化基本面的隐含估计。每个人都赢了！
【讨论】：
我假设您的所有数字都是 整数 值的倍数。在我接下来的解释中，A 将表示您尝试查找的“根”频率，B 将是您必须开始的数字数组。
您正在尝试做的事情表面上类似于linear regression。您正在尝试找到一个线性模型 y=mx+b 来最小化线性模型和一组数据之间的平均距离。在您的情况下，b=0，m 是根频率，y 表示给定值。最大的问题是没有明确给出自变量 X。关于 X，我们唯一知道的是它的所有成员都必须是整数。
您的首要任务是尝试确定这些自变量。目前我能想到的最佳方法是假设给定频率具有几乎连续的索引 (x_1=x_0+n)。所以B_0/B_1=(x_0)/(x_0+n) 给出了一个（希望是）小整数 n。然后，您可以利用 x_0 = n/(B_1-B_0) 的事实，从 n=1 开始，并不断提高它，直到 k-rnd(k) 在某个阈值内。在获得 x_0（初始索引）之后，您可以近似根频率 (A = B_0/x_0)。然后您可以通过查找x_n = rnd(B_n/A) 来近似其他索引。这种方法不是很健壮，如果数据中的错误很大，很可能会失败。
如果您想要更好地逼近根频率 A，您可以使用线性回归来最小化线性模型的误差，因为您已经有了相应的因变量。最简单的方法是使用最小二乘拟合。 Wolfram's Mathworld 对此问题进行了深入的数学处理，但通过谷歌搜索可以找到 fairly simple explanation。
【讨论】：
有趣的问题...不容易。
我想我会看看样本值的比率：
然后我会在这些结果中寻找一个简单的整数比率。
我还没有追查过，但是在某个地方，您认为 1:1000 或其他的误差就足够了，然后您回溯以找到基本的近似 GCD。
【讨论】：
我自己看到和使用的解决方案是选择一个常数，比如 1000，将所有数字乘以这个常数，将它们四舍五入，使用标准算法找到这些整数的 GCD，然后将结果除以所述常数（1000）。常数越大，精度越高。
【讨论】：
1.234,2.468 给出了1.234，但1.234,2.467 给出了0.001。
当您先验地选择 3 个正公差（e1、e2、e3）时，这是 shsmurfy 解决方案的重新制定
然后问题是搜索最小的正整数 (n1,n2,n3) 和最大的根频率 f，这样：
f1 = n1*f +/- e1
f2 = n2*f +/- e2
f3 = n3*f +/- e3
我们假设 0 如果我们修复 n1，那么我们会得到这些关系：
f is in interval I1=[(f1-e1)/n1 , (f1+e1)/n1]
n2 is in interval I2=[n1*(f2-e2)/(f1+e1) , n1*(f2+e2)/(f1-e1)]
n3 is in interval I3=[n1*(f3-e3)/(f1+e1) , n1*(f3+e3)/(f1-e1)]
我们从 n1 = 1 开始，然后递增 n1 直到区间 I2 和 I3 包含一个整数 - 即 floor(I2min) different from floor(I2max) 与 I3 相同
然后我们在区间I2中选择最小整数n2，在区间I3中选择最小整数n3。
假设浮点误差为正态分布，根频率 f 的最可能估计是最小化的那个
J = (f1/n1 - f)^2 + (f2/n2 - f)^2 + (f3/n3 - f)^2
那是
f = (f1/n1 + f2/n2 + f3/n3)/3
如果在区间 I2,I3 中有多个整数 n2,n3，我们也可以选择最小化残差的对
min(J)*3/2=(f1/n1)^2+(f2/n2)^2+(f3/n3)^2-(f1/n1)*(f2/n2)-(f1/n1)*(f3/n3)-(f2/n2)*(f3/n3)
另一个变体可能是继续迭代并尝试最小化另一个标准，例如 min(J(n1))*n1，直到 f 低于某个频率（n1 达到上限）...
【讨论】：
我在 MathStackExchange（here 和 here）中找到了这个问题的答案。
我只设法（还）测量了给定谐波频率列表（遵循声音/音乐命名法）的基频的吸引力，如果你有一个减少的选项的数量，并且可以计算每个选项的吸引力，然后选择最合适的选项。
我在 MSE 中的问题中的 C&P（格式更漂亮）：
目标是找到最大化吸引力的 x。这是您的示例 [2.468, 3.700, 6.1699] 的 (gcd_appeal) 图，您发现最佳 GCD 位于 x = 1.2337899957639993
编辑： 您可能会发现这段 JAVA 代码很方便计算除数相对于红利列表的（模糊）divisibility（又名 gcd_appeal）；你可以用它来测试你的哪个候选人是最好的除数。代码看起来很丑，因为我试图优化它的性能。
//returns the mean divisibility of dividend/divisor as a value in the range [0 and 1]
// 0 means no divisibility at all
// 1 means full divisibility
public double divisibility(double divisor, double... dividends) {
double n = dividends.length;
double factor = 2.0 / divisor;
double sum_x = -n;
double sum_y = 0.0;
double[] coord = new double[2];
for (double v : dividends) {
coordinates(v * factor, coord);
sum_x += coord[0];
sum_y += coord[1];
}
double err = 1.0 - Math.sqrt(sum_x * sum_x + sum_y * sum_y) / (2.0 * n);
//Might happen due to approximation error
return err >= 0.0 ? err : 0.0;
}
private void coordinates(double x, double[] out) {
//Bhaskara performant approximation to
//out[0] = Math.cos(Math.PI*x);
//out[1] = Math.sin(Math.PI*x);
long cos_int_part = (long) (x + 0.5);
long sin_int_part = (long) x;
double rem = x - cos_int_part;
if (cos_int_part != sin_int_part) {
double common_s = 4.0 * rem;
double cos_rem_s = common_s * rem - 1.0;
double sin_rem_s = cos_rem_s + common_s + 1.0;
out[0] = (((cos_int_part & 1L) * 8L - 4L) * cos_rem_s) / (cos_rem_s + 5.0);
out[1] = (((sin_int_part & 1L) * 8L - 4L) * sin_rem_s) / (sin_rem_s + 5.0);
} else {
double common_s = 4.0 * rem - 4.0;
double sin_rem_s = common_s * rem;
double cos_rem_s = sin_rem_s + common_s + 3.0;
double common_2 = ((cos_int_part & 1L) * 8L - 4L);
out[0] = (common_2 * cos_rem_s) / (cos_rem_s + 5.0);
out[1] = (common_2 * sin_rem_s) / (sin_rem_s + 5.0);
}
}
【讨论】：