計算機是怎樣進行開方和冪運算的？

01-04

比如 $sqrt{2.5}$ 或 $5^{12.34}$ 等。

那個說迭代的，你弄錯了一個非常重要的東西：pow 的實現必須是 O(1) 的，而且最好能用電路完成，所以我們需要利用一些數值技術。這裡不考慮邊界情況，只處理「普通」正浮點數的指數運算。

先明確個概念：

正浮點數的 IEEE-754 表示：IEEE-754 中的所有正浮點數都是用如下形式表示： $x=2^{E_x}(1+F_x)$ ，其中整數 $E_x$ 為其指數部分， $F_xin[0,1)$ 為小數部分。

然後，因為 $x^y=2^{ylog_2 x}=2^{yE_x+ylog_2(1+F_x)}$ ，因此我們只需要算出 $yE_x+ylog_2(1+F_x)$ 就可以用泰勒級數算出 $x^y$ 了。

因為 $F_xin[0,1)$ ， $log_2(1+F_x)in[0,1)$ ，此結果可以用泰勒級數算出，我們用 $J$ 表示之。將 $y$ 表為 $2^{E_y}(1+F_y)$ 之後我們有 $y[E_x+log_2(1+F_x)]=2^{E_y}(1+F_y)(E_x+J)$ ，後面兩個因子的乘積是個比較好處理的量（範圍確定），於是我們就能算出它來，繼而得到冪 $x^y$ 的數值。

Julia 所用的 openlibm 中是這樣實現的：https://github.com/JuliaLang/openlibm/blob/master/src/e_pow.c

如果是x86/x87架構的CPU，是用FSQRT指令實現32位浮點數的開平方，用FYL2X和F2XM1指令來實現任意次冪的。FSQRT指令內部採用牛頓迭代法原理實現[1]。如果用軟體實現，當然可以採用更簡潔的方法，但人看起來簡潔不等於用指令實現簡潔[2]。硬體的實現思路和數學方法還是有一定差距的，這一點請 @Belleve 注意。

一、牛頓法求根的原理

我們要求給定數 $a$ 的平方根 $x$ （ $sqrt{a}=x$ ），就相當於解方程： $x^2-a=0$ 的根。其中 $a$ 為已知數， $x$ 為未知數。令 $f(x)=x^2-a$ ，也就是要求使得 $f(x)=0$ 的 $a$ 值，即該函數與 $x$ 軸的交點。這個函數顯然在 $x>0$ 的情況下為單調遞增函數：

我們可以先假設一個 $x$ 的解 $x_0$ ，然後通過一次次的修正這個值，也就是逼近來獲得最終的 $x$ 值。

$f(x)$ 在任意一點 $x$ 的斜率為其導函數： $f$ （求導公式），故 $x_0$ 處的斜率為 $2x_0$ 。用點斜式過 $P_0$ 作直線交 $x$ 軸於 $x_1$ 點，此時 $x_1$ 更接近於要求的 $x$ 點了。同樣的方法再求得 $x_2$ 、 $x_3$ ……。

當然到最後求得的解之間的差越來越小，意味著與我們要求的點已經十分接近了。根據精度的需要我們可以設定一個值，當兩個解之間的差小於這個值我們就停止運算，此時算得的解就是在滿足精度要求下的解。

二，求任意次冪pow(x, p)

由於 $x^p=left(e^{ln x} ight)^p=e^{pln x}$ ，然後用FYL2X指令計算 $pln x$ ，再用F2XM1指令計算 $e^{pln x}$ 即可得解[3]。FYL2X指令內部則是令 $ln x=frac{log_2x}{log_2e}$ ，而 $log_2x$ 的值可用迭代求解，過程參見[4]。又因為 $log_2e$ 是個常數， $p$ 為輸入值，所以 $pln x$ 即用除法和乘法求出。F2XM1指令則又是一次迭代的過程。

[1] http://en.wikipedia.org/wiki/Computational_complexity_of_mathematical_operations

[2] http://www4.wittenberg.edu/academics/mathcomp/bjsdir/zusez3squarerootalgorithm.htm

[3] http://pastebin.com/VWfE9CZT

[4] https://en.wikipedia.org/wiki/Binary_logarithm#Recursive_approximation

實際上是用牛頓迭代法求的.

牛頓迭代法求平方根

恕不贅述

在數值計算權威Mathematica中：

（參考Some Notes on Internal Implementation ）

如果是整數的power 運算，進行所謂的binary exponentiation。也就是把任意次數的power 轉化成平方計算。

例如：

$5^{81}=5*5^{80}=5*(5^{40})^2=5*((5^{20})^2)^2 =5*(((5^{10})^2)^2)^2=5*((((5^5)^2)^2)^2)^2=5*((((5*5^4)^2)^2)^2)^2=5*((((5*(5^2)^2)^2)^2)^2)^2$

這樣乘法的次數基本縮減到了log(n),n是指數的大小。

def exponen(base, exponent): if base==1: return 1 if base==0: return 0


    result=1
    while exponent&>0:

        if exponent1:

            result=result*base

        exponent=exponent&>&>1

        base=base**2

return result

如果是開方或者非整數次方，採用的是牛頓法。

這種事都是數學解決的

就題目而言，開方使用的是【牛頓迭代法】，原理你可以搜索

幾月前為了回答某個問題，我曾做過一個玩具，當時就用【牛頓迭代】實現了開方

原答案已經刪了，土豆還留了幾個視頻，感興趣可以一看

畫線下部分_土豆

開方函數如下截圖

關於開根號的最偉大的傳說便是《雷神之錘》裡面的方法了，參見：

一個Sqrt函數引發的血案

論文參見： http://www.matrix67.com/data/InvSqrt.pdf

牛頓太慢了，能算出(1,2)就能算出來任意IEEE754的double，所以只寫一個(1,2)範圍內的開方，供參考

#1&

def mysqrt(p):

x = 0.44333414+p*(0.62925810-p*0.07211624)

d = 0.25*( p / (x*x) - 1)

return x*(1+2*d*(1-d*(1-d*(2-5*d))))

結果

$ python sqrt.py
p pow(p,0.5) mysqrt(p)

1.00 1.000000000000000 1.000000000000000

1.10 1.048808848170152 1.048808848170151

1.20 1.095445115010332 1.095445115010332

1.30 1.140175425099138 1.140175425099138

1.40 1.183215956619923 1.183215956619923

1.50 1.224744871391589 1.224744871391589

1.60 1.264911064067352 1.264911064067352

1.70 1.303840481040530 1.303840481040530

1.80 1.341640786499874 1.341640786499874

1.90 1.378404875209022 1.378404875209022

64位linux下的數學庫尤其是double版本為了達到ieee754的最後一位精度要求都是用c演算法實現的，包括intel平台的。一般是什麼泰勒展開之類。

你捜fdlibm行了。

不過比較複雜，不容易看明白。

看了你就明白了，一個工業級的數學庫是多麼的複雜。

我是來把 bhuztez 的python代碼轉換成公式的：

$a^{b} = (e^{ln(a)})^{b} = e^{b imes ln(a)}$

2**2=4