Math Alg. Root 数学算法笔记：开根号

from 开根号基础公式 - 百度文库

from 开方（数学术语）_百度百科

对于任意实数的开方，可以使用切线法得到其任意精度的结果，切线法的迭代公式为：

取任意初始值

以上迭代序列将会收敛：

实际应用中一般取初始值为稍微大 

的实数，这样可以加快序列的收敛速度。

e.g. （C）

// 2015-12-24
// By: ChenYu
#include "math.h"
#include "stdio.h"
#define ABS(a) ((a)<0?-(a):(a))
#ifdef _WIN32
    typedef unsigned __int64 uint64;
#else
    typedef unsigned long long uint64;
#endif
// calculate a approximate value
static double calcInitRoot(double x, int n)
{
    const uint64 exptMask=((uint64)1<<11)-1;
     const uint64 fracMask=((uint64)1<<52)-1;
    uint64 xInt=*(uint64*)&x;
    int xExpt=(int)((xInt>>52)&exptMask)-1023;
    xInt=((uint64)((xExpt+1024*n-1)/n)<<52)+(xInt&fracMask)/n;
    return *(double*)&xInt;
}
 
double calcRoot(double x, int n)
{
    int i, j, s=1-((x<0)<<(n&1));
    double a=ABS(x);
    double x1, x0=calcInitRoot(a, n);
    double err=x0*1e-14;
    if(x==0)
        return 0;
    for(i=1; i<50; i++)
    {
        double xn=1;
        for(j=0; j<n-1; j++)
            xn*=x0;
        x1=((n-1)*x0*xn+a)/(xn*n);
        // printf("x%d=%.14f\n", i, x1);
        if(ABS(x1-x0)<=err)
            break;
        x0=x1;
    }
    return s*x1;
}
 
void main()
{
    double x=-31141.592653589793;
    int n=11;
    double y=calcRoot(x, n);
    printf("root(%g,%d)=%+.14f\n", x, n, y);
    printf("root(%g,%d)=%+.14f\n", x, n, pow(ABS(x), 1.0/n));
}