0.math庫裡的實現（庫裡？）

在寫這個函式之前，我想先查查c語言math庫裡自帶的pow是怎麼實現的
我先是在math.h裡找，結果只找到了這個

_CRTIMP	double __cdecl pow (double, double);

和這個

/* 7.12.7.4 The pow functions. Double in C89 */
extern float __cdecl powf (float, float);
#ifndef __NO_INLINE__
__CRT_INLINE float __cdecl powf (float x, float y)
  {return (float) pow (x, y);}
#endif
extern long double __cdecl powl (long double, long double);
 

您這隻給個宣告不給我函式定義是幾個意思啊
既然庫檔案查不到，我就上谷歌裡搜搜吧
突然，不知怎的，我有一種不好的預感
果然，我啥也沒查到，倒是有一大堆人在那裡研究pow是怎麼實現的
更有甚者，只寫了一個for迴圈，然後給出示例

double value3 = mypow(5.21,4.11);
printf("value3 = %f\n", value3);

我還特意用手機算了一下5.21的4.11次冪

您說您這不是自欺欺人嗎？
但我就好奇了，math.h的原始碼到底定義在哪裡呢？
終於在知乎上找到了答案，感謝Belleve大佬，點選這裡檢視

math.h 裡的函式都是定義在 libm 裡，而每個 libm 實現都不同
gcc 的 glibm 中數學函式的實現完全是平臺依存的，在 x86 機器上，能呼叫 FPU 指令的就用 FPU（比如 sqrt() 就實際上呼叫 FSQRT，log() 呼叫的是 FYL2X），否則再自己實現
如果需要軟體實現，方法基本上是泰勒級數。當然對於 sin 這類，可以用專門的優化演算法，比如 CORDIC
CPU 中的電路基本上也就是泰勒級數。

原來是個叫libm的東西，gcc的就叫glibm
那glibm在哪裡呢？在glibc，也就是c標準庫裡
是在是不想再在mac裡找各種檔案了，我再重新下載一個吧
libm.so位於 libc6-dev 這個包裡，那就到pkgs.org上找找
結果網速感人，5MB的檔案可能要下十分鐘，在這個空檔，我又搜了一下x86浮點運算指令集
果然不好好學彙編是要付出代價的，這些用軟體實現費勁的函式，幾條彙編指令就搞定了，比如說以下幾條

我又想到儲存浮點數的IEEE標準，這個計算乘方和log以2為底真是引人深思啊
既然浮點數在暫存器裡就是指數形式儲存的，直接在[1:8]位操作不就行了？再加上上面的浮點運算指令，搞定！
這時候檔案也下載好了，開始唸咒語（Windows的同學別學我哈，容易被室友當成瘋子）

ar -vx libc6-dev_2.28-2_i386.deb   
tar -xzvf data.tar.gz  

然後肉眼搜尋（一個一個資料夾點開……）
終於在usr/lib/i386-linux-gnu裡找到了libm.so
幹啊，.so的檔案我還打不開，.so惹不起.a總能開啟吧
敲程式碼！

#include <stdio.h>
int main() {
    FILE *fp = fopen("libm.a","rb");
    char element;
    
    while(!feof(fp)){
        element = getc(fp);
        printf("%2x",element);
    }
    
    fclose(fp);
    
    return 0;
}

然而也僅僅是打開了而已，有啥用呢，讀起來要累死人
還是用synalysis吧，Synalyze It!
結果對著一大堆二進位制檔案synalyze半個小時，體驗極其糟糕
看到了一大堆ieee754，和各種powf32，powf64
因為最後實在讀不下去了，就姑且認為我的猜測就是正確的吧（我咋怎麼懶）
要是有哪位大佬看到了pow在libm.a裡的實現，請務必告訴我！！！
其實我還想知道像f2xm1，fyl2x這樣的指令在底層（閘電路層次）是如何實現的，計算系統基礎紮實的同學請務必帶帶我！！！

不扯淡了，我們現在切入正題，嘗試用不同的方法實現pow函式

1.牛頓迭代法實現pow

其實這個方法是很容易想到的，特別是在嘗試用牛頓迭代法實現sqrt函式之後，很自然的聯想到也可以用這個方法解決pow函式的實現
首先考慮計算pow(2.17,3.14)
我們可以用for迴圈實現pow(2.17,3)
接下來只要把pow（2.17，0.14）也計算出來，再把兩部分乘到一起，問題就解決了
後面的那部分，也就是一個數的小數次冪怎麼算呢？
我再分
將0.14分成0.1和0.04，0.04就是0.01*4
由牛頓迭代法，我們就可以求出pow（2.17,0.1）和pow（2.17,0.01），以及整數次冪pow（2.17,3）,這樣我們就可以把這個值算出來
還是稍微寫一下原理吧
計算2.17^0.1，即解出方程x¹⁰=2.17的根，然後百度百科吧……

用牛頓迭代法解非線性方程，是把非線性方程 f(x)=0線性化的一種近似方法。把 f(x)在點x0的某鄰域內展開成泰勒級數

取其線性部分（即泰勒展開的前兩項），並令其等於0，即 f(x0)+f’(x0)(x-x0)=0，以此作為非線性方程f(x)=0的近似方程，若f’(x0)!=0，則其解為

這樣，得到牛頓迭代法的一個迭代關係式

下面還有一些細節問題，例如當pow(x,y)中的x為負數，y為整數時，函式會正常工作，但y若是小數，就會返回nan
那我們根據x和y的取值列一個表格

但是實際上我們不能製造出nan，於是可以將0.000000作為nan輸出
以及pow函式的精度問題，我們可以通過控制y具體計算到小數點後幾位來決定pow保留的精度
先寫個小函式求絕對值

double Absolute(double a){
	if(a<0){
		a=-1*a;	
	}
	return a;
}

再寫個函式，計算x的y次冪，這個y是像0.1，0.01，0.001這樣的小數，實際上，這樣的小數可以寫成10^-n的形式，我們只需要傳整數n到函式中

double func1(double x, int n){
    int number=1;
    int i,j;
    double x1n=1,x1=1;                              	//x1n表示x1的number次冪
    
    for(i=0;i<n;i++){                           		//number即10^n
        number=number*10;
    }
    
    while(Absolute(x1n-x)>0.0001){           	//通過牛頓迭代法計算t^number=x(t為未知數)的根
        x1n=1;
        for(j=0;j<number;j++){
            x1n=x1n*x1;
        }

        x1=((number-1)*x1n+x)/(number*x1n/x1);
   
        x1n=1;
        for(j=0;j<number;j++){
            x1n=x1n*x1;
        }                                          		 	//經過多次迭代，x1n的值逼近方程的根t
    }
    
    return x1;
}

接下來只要再寫一個double pow1(double x, double y)函式，對x和y對取值各種if，然後用兩個for迴圈分別處理y的整數部分和小數部分，其餘的工作交給我們之前寫的func1
下面是程式碼

double pow1(double x, double y){
    int positive = 1;                                   //用來儲存y的正負
    double answer=1;                                  //即結果
    
    if(x>0&&y>0){
        positive=1;
    }
    else if(x!=0&&y==0){
        return 1;
    }
    else if(x>0&&y<0){
        y=-y;
        positive=0;
    }
    else if(x==0&&y>0){
        return 0;
    }
    else if(x==0&&y<=0){
        return 1/0;
    }
    else if(x<0&&(y-(int)y<0.0001||(y-(int)y>0.999))){
        if(y>0){
            positive=1;
        }
        else{
            y=-y;
            positive=0;
        }
    }
    else {
        return 1/0;
    }
    
    int integer=(int)y;
    int i,j;
    
    for(i=0;i<integer;i++){
        answer=answer*x;
    }
    
    if((y-(int)y<0.0001||y-(int)y>0.999)&&positive==1){                   //如果y為整數，跳過小數運算部分
        return answer;
    }
    if((y-(int)y<0.0001||y-(int)y>0.999)&&positive==0){
        return 1/answer;
    }
    
    double decimal=y-(int)y;
    
    for(i=1;i<=ACCURACY;i++){      							//ACCURACY表示計算精度
        decimal=decimal*10;
   
        for(j=0;j<(int)decimal;j++){
            answer=answer*func1(x, i);
        }
        decimal=decimal-(int)decimal;
    }
    
    if(positive==1){
        return answer;
    }
    else {
        return 1/answer;
    }
}

把這三個函式拼起來就完工了！
這裡還有一個小細節是

x<0&&(y-(int)y<0.0001||(y-(int)y>0.999))

這是因為double型儲存形式與int有本質上的區別（ieee754）
比如說1這個double，在記憶體中以二進位制形式儲存，轉化為十進位制可能是1.000000001也可能是0.99999999
而強制型別轉換是隻保留double型小數點前的部分，不是四捨五入
因此很可能(int)1的結果是0
所以我們要用(y-(int)y<0.0001||(y-(int)y>0.999)讓兩種可能都通過

下面讓我們計算以下正版的5.21的4.11次冪

printf("%lf",pow1(5.21,4.11));

快看控制檯！

883.492894Program ended with exit code: 0

和手機算出來的資料基本吻合，增大ACCURACY的值還會提高精度
（我還看到有的博主計算y為整數的情形時用了分治的方法，我懶得寫了）

2.冪級數展開實現pow

為了避免歧義，這次pow接受的引數改為a,b
即double pow2(double a ,double b)
我採用麥克勞林展開的方法，用下面的公式

在pow(a,b)中，令t=a-1，設f(t)=(t+1)^b
就可以得到
f(t)=(t+1)^b=1+bt+b(b-1)t²/2+…+b(b-1)(b-2)…(b-n+1)tⁿ/n!+…
這就把小數次冪的問題轉化成了計算整數次冪和計算階乘
然而，還存在一個問題，那就是當|t|>1時冪級數是發散的，隨著我們設定的ACCURACY的增大，計算結果反而更加不精確
可以這樣，當a>2時，對a取倒數，對應的t=a-1就滿足|t|<=1了
然而，直覺告訴我這事不靠譜
其實我們可以對a不停的除2，直到0<a<2為止，再把相應的指數乘到b上，這樣算精度應該會高一些
最後就是針對a,b的一些極端情況用if else縫縫補補
直接上程式碼！

#include <stdio.h>
#define ACCURACY 100
double func1(double t,int n);
double func2(double b,int n);
double pow2(double a,double b);
int main() {
    printf("%lf",pow2(5.21,4.11));
    return 0;
}

double pow2(double a,double b){
    if(a==0&&b>0){
        return 0;
    }
    else if(a==0&&b<=0){
        return 1/0;
    }
    else if(a<0&&!(b-(int)b<0.0001||(b-(int)b>0.999))){
        return 1/0;
    }

    if(a<=2&&a>=0){
        double t=a-1;
        double answer=1;
        for(int i=1;i<ACCURACY;i++){
            answer=answer+func1(t,i)*func2(b,i);
        }
        return answer;
    }
    
    else if(a>2){
        int time=0;
        
        while(a>2){
            a=a/2;
            time++;
        }
        
        return pow2(a,b)*pow2(2,b*time);
    }
    
    else{
        if((int)b%2==0){
            return pow2(-a,b);
        }
        else {
            return -pow2(-a,b);
        }
    }
}
double func1(double t,int n){
    double answer=1;
    for(int i=0;i<n;i++){
        answer=answer*t;
    }
    
    return answer;
}
double func2(double b,int n){
    double answer=1;
    for(int i=1;i<=n;i++){
        answer=answer*(b-i+1)/i;
    }
    
    return answer;
}

最後再試一下5.21的4.11次冪

printf("%lf",pow2(5.21,4.11));

見證奇蹟的時刻

883.492857Program ended with exit code: 0

嗷嗷，終於完成了，這個的精確度更高，原因是ACCURACY取了100