禁止使用sqrt等返回浮点数的函数，如何最高效的得到最小的不小于给定正整数的完全平方数？第1页

gao-tian-50 网友的相关建议:

为了方便讨论，我们假设输入值为i，输出值为result^2。同时我们不妨默认i是一个大于等于1的整数，忽略corner case们。

题目的要求就是找到一个result，使得(result-1)^2 < i <= result^2。同时根据题目的要求，我们在整个计算过程中不引入任何浮点数（否则我可以手写一个等价于sqrt的浮点数函数）。

这道题目的naive解是O(N^0.5)的，从1开始一个一个试，每试一次的代价是固定的，试到ceil(sqrt(i))，也就是其平方根的上取整的时候，就搞定了，一共需要试N^0.5次

       uint32_t naive(uint32_t i) {     for (uint32_t result = 1; result <= i; result ++) {         if ((result - 1) * (result - 1) < i &&                  i <= result * result) {             return result * result;         }     } }

如果你在面试的时候遇到这道题，那上面的方法就显得太low了。这道题面试的正解是binary search。我们知道需要的result一定在1和i之间，因此通过二分法测试，我们可以把测试的次数减小到O(log N)

       uint32_t bin_search(uint32_t i) {     uint32_t head   = 1;     uint32_t tail   = i;     uint32_t result = 0;          if (tail >= (1 << 16)) {         // 这里记得检查一下overflow         tail = (1 << 16) - 1;     }      result = (head + tail) / 2;      while (1) {         if ((result - 1) * (result - 1) >= i) {             tail = result - 1;             result = (head + tail) / 2;         } else if ( i > result * result) {             head = result + 1;             result = (head + tail) / 2;         } else {             return result * result;         }     }

一般来说，程序员面试答到这里就过关了，但是作为一个掌握了数学这个伟大工具的程序员，我们还知道一个方法叫做牛顿法。这道题本质上还是求sqrt，我们可以利用牛顿法逼近，同时把数据处理依然控制在整型范围（避免使用浮点数）。由于牛顿法逼近sqrt(i)可能从左也可能从右，同时因为我们做的都是整型计算，可能会出现误差问题，在测试的时候我们同时检查两个数。

       uint32_t test(uint32_t input, uint32_t result) {     if (result >= (1 << 16)) {         // 这里也要注意overflow的问题         return 0;     } else if ((input > (result-1) * (result-1)) && (input <= result * result)) {         return result;     } else if (input > result * result && input <= (result + 1) * (result + 1)) {         return result + 1;     }     return 0; }  uint32_t newton(uint32_t i) {     uint32_t guess_val = 1;     uint32_t loop = 0;     while (1) {         uint32_t result = test(i, guess_val);         loop += 1;         if (result) {             return result * result;         } else {             guess_val = (i / guess_val + guess_val) / 2;         }     } }

既然我们函数都写出来了，不妨跑一个benchmark测试吧！为了不等太久，我们先把计算数量级设置到(1<<22)。我们对1到(1<22)之间的所有整数去计算，看看三个算法分别需要多久。

首先出场的是我们的naive选手，经过了一通狂跑，他取得了13.559s的好成绩

       gaotian@cercis:/tmp/funprog$ time ./test   real 0m13.559s user 0m13.559s sys 0m0.000s

接下来大摇大摆目中无人地出场的是我们的工业界领袖bin_search选手，他以碾压的姿态，取得了0.198s的成绩，完爆naive选手。

       gaotian@cercis:/tmp/funprog$ time ./test   real 0m0.199s user 0m0.199s sys 0m0.000s

最后出场的就是万众瞩目的学术界新型，newton选手，他在摩拳擦掌之后，取得了——0.513秒的成绩。

       gaotian@cercis:/tmp/funprog$ time ./test   real 0m0.513s user 0m0.513s sys 0m0.000s

什么！我不信！不阔能！怎么可能牛顿法没有愚蠢的二分法好呢？

别着急，我们慢慢分析。

首先，应newton选手的要求，我们把数据级别提升到了(1<<25)的级别，然而bin_search依然以相似的优势领先

       gaotian@cercis:/tmp/funprog$ time ./test   real 0m1.574s user 0m1.570s sys 0m0.004s  gaotian@cercis:/tmp/funprog$ time ./test   real 0m4.657s user 0m4.657s sys 0m0.000s

newton选手思考了一下，他觉得自己的差距可能是在test函数上。这个函数本身有消耗，同时在大部分情况下都要走过所有的判断，所以他决定改写一下这个函数

       static inline uint32_t test(uint32_t input, uint32_t result) {     if ((result >= (1 << 16)) ||              (result+1) * (result+1) < input ||             (result-1) * (result-1) >= input) {         return 0;     } else if (input <= result * result) {         return result;     } else {         return result + 1;     } }

沾沾自喜的牛顿选手重新提交了程序，结果是——

       gaotian@cercis:/tmp/funprog$ time ./test   real 0m4.671s user 0m4.670s sys 0m0.000s

甚至还慢了一点……

看着哇哇大哭的newton选手，来自Green Hills的MULTI老师缓缓走来，对他说，孩子，听我一句话，叫做“Never optimize before profilng”。当你不知道你程序的瓶颈在哪儿的时候，你要先找到自己程序的问题。来，请欣赏我的表演。

Green Hills的MULTI把newton选手和bin_search选手的程序同时编译了，在编译的时候打开了coverage功能，然后运行了一下。

可以看到，newton选手的方法在收敛次数上和bin_search是相似的，516m次 vs 503m次，甚至newton还多一些。这和我们的幻想并不一致。

其原因是什么呢？很简单，尽管newton法可以不断提升精度，即结果离理论值的距离，但是在离结果较远的时候，它的收敛速度并不算优异。而对这个问题来说，远方的收敛速度才是决定性的，在收敛到+-1之后，结果就出来了，而这里才是newton法开始发挥的空间。

我们可以简单地做个数学计算。bin_search在开始阶段从x逼近sqrt(x)的时候，趋势是x/2, x/4, x/8。而newton法是(i/x + x/2)，也就是除二之后还要加一个数。这个逼近速度是要逊于bin_search的（当然实际上因为这个数比较小，是差不多的）。只有当bin_search overshoot了之后，newton的渐进优势才能体现出来。

同时其实我们可以知道，用newton method的时候，如果可以事先取得一个比较好的猜测，那逼近就会变得容易很多。所以我们可以通过如下的方式去先“猜”一个数。我们把i不断除2，直到除到它的平方比i小为止。然后我们以此作为第一个猜测。

       uint32_t guess(uint32_t i) {     int len = 0;     uint32_t temp = i;     while (temp * temp > i) {         temp >>= 1;     }      return temp; }

       gaotian@cercis:/tmp/funprog$ time ./test   real 0m1.611s user 0m1.611s sys 0m0.000s

于是newton在数学的帮助下得到了一个和bin_search速度接近的代码。为什么速度接近呢？仔细想一下，这个猜测的过程，其实就是bin_search的前半部分啊！

天真的我曾以为，故事到了这里就有了结尾。

我没有去思考为什么在复制了bin_search的前半部分之后，newton还是慢一些。我想当然地把因素归结到newton本身的复杂机制，但是——我错了，我把代码写错了，而且是一个我在前面强调了好几次的地方……overflow……

我们把guess函数重新写一下，考虑进去overflow的因素

       uint32_t guess(uint32_t i) {     int len = 0;     uint32_t temp = i;     if (temp >= (1<<16)) {         temp = (1<<16) - 1;     }     while (temp * temp > i) {         temp >>= 1;     }      return temp; }

然后再跑一下。

       gaotian@cercis:/tmp/funprog$ time ./test  real    0m1.051s user    0m1.051s sys     0m0.000s

牛顿威武！

不过细心的小伙伴应该看出来了，我们这里其实耍赖皮了，我们直接默认开方之后的数最大就是(1<<16)。如果newton可以这么玩赖，那bin_search当然也可以。正当我回去准备把bin_search的代码也加入这个赖皮因素再比较的时候，我才意识到——我已经做了这件事了……

bin_search里已经把tail值在i太大的时候锁定在(1<<16)，而我没意识到这件事！而这个优化很可能才是bin_search之前领先的核心因素，它少做了好多个iteration啊！

这个故事告诉我们，即便你觉得自己花了很多时间和精力把一个程序研究透了，你依然可能有不经意间留下的bug和错误，等着你未来打自己的脸。

认同（划掉）。

不认同！！

禁止使用sqrt等返回浮点数的函数，如何最高效的得到最小的不小于给定正整数的完全平方数？的其他答案点击这里

禁止使用sqrt等返回浮点数的函数，如何最高效的得到最小的不小于给定正整数的完全平方数？第1页

相关话题

前一个讨论

下一个讨论

相关的话题

禁止使用sqrt等返回浮点数的函数，如何最高效的得到最小的不小于给定正整数的完全平方数？ 第1页

相关话题

前一个讨论

下一个讨论

相关的话题

禁止使用sqrt等返回浮点数的函数，如何最高效的得到最小的不小于给定正整数的完全平方数？第1页