Fukushima, Toshio (2012) Numerical computation of spherical harmonics of arbitrary degree and order by extending exponent of floating point numbers. J Geodesy 86:271–285.Numerical computation of spherical harmonics of arbitrary.pdf

还有什么微分方程组不能数值解？
位错？要多少阶有多少阶......

**************************************************************************************************
**************************************************************************************************
写了一个XNumber 类，计算时可以直接定义XNumber类型的变量。
示例：
XNumber e1(1.0,0)
XNumber e2,e3;
e2=XNumber(2.0,0);
e3=e1+e2;
e3=e1-e2;
e3=e1*e2;
e3=e1/e2;
e3=e3.sqrt2(e1*e2+e2-e1/e2);
e3=e3.Norm(e3);

//////////////////////////////////////////////////////////////////////////

#include "math.h"
class XNumber
{
private: //指数，底，从10进制数换算成BIG进制数，BIG is big
        int IND; 
        double BIG; 
        double BIGI;
        double BIGS; //BIG的平方根
        double BIGSI;
public: 
        double x; //XNumber=x*B**ix
        int ix; //指数部分
public:
        XNumber(); 
        XNumber(double g,int ig);
        XNumber operator+(XNumber &xn); //重载加法
        XNumber operator-(XNumber &xn);  //重载减法 
        XNumber operator*(XNumber &xn); //重载乘法
        XNumber operator/(XNumber &xn); //重载除法 
        XNumber Norm( XNumber xn);  //规格化
        double x2f(XNumber d);  //转换成10进制数
        XNumber sqrt2(XNumber c);  //重载开根号
};
XNumber::XNumber()
{
        IND=960;
        BIG=pow(2.0,IND);
        BIGI=pow(2.0,-IND);
        BIGS=pow(2.0,IND/2);
        BIGSI=pow(2.0,-IND/2);
 }
XNumber::XNumber(double g,int ig)
{        
        IND=960;
        BIG=pow(2.0,IND);
        BIGI=pow(2.0,-IND);
        BIGS=pow(2.0,IND/2);
        BIGSI=pow(2.0,-IND/2);
        x=g;
        ix=ig;
}
XNumber XNumber::operator+(XNumber &xn)
{
        XNumber c;

        int id=ix-xn.ix;
        if(id==0) 
        {
                c.x=x+xn.x;
                c.ix=ix;
        }
        else if(id==1) 
        {
                c.x=x+(xn.x*BIGI);
                c.ix=ix;
        }
        else if(id==-1) 
        {
                c.x=xn.x+x*BIGI;
                c.ix=xn.ix;
        }
        else if(id>1)
        {
                c.x=x;
                c.ix=ix;
        }
        else
        {
                c.x=xn.x;
                c.ix=xn.ix;
        }
        c=c.Norm (c);
        return c;
}
XNumber XNumber::operator-(XNumber &xn)
{
        XNumber c;
        int id=ix-xn.ix;
        if(id==0) 
        {
                c.x=x-xn.x;
                c.ix=ix;
        }
        else if(id==1) 
        {
                c.x=x-(xn.x*BIGI);
                c.ix=ix;
        }
        else if(id==-1) 
        {
                c.x=-xn.x+x*BIGI;
                c.ix=xn.ix;
        }
        else if(id>1)
        {
                c.x=x;
                c.ix=ix;
        }
        else
        {
                c.x=-xn.x;
                c.ix=xn.ix;
        }
        c=c.Norm(c);
        return c;
}

XNumber XNumber:: operator*(XNumber &c)
{        
        XNumber d;
        d.x=x*c.x;
        if(d.x==0.0) 
        {
                d.ix=0;
                return d;
        }
        d.ix=ix+c.ix;
        if(fabs(d.x)>=BIGS)
        {
                d.x=d.x*BIGI;
                d.ix=d.ix+1;
        }
        else if(fabs(d.x)<BIGSI)
        {
                d.x=d.x*BIG;
                d.ix=d.ix-1;
        }
        d=d.Norm (d);
        return d;
}
XNumber XNumber:: operator/(XNumber &c)
{        
        XNumber d;
        d.x=x/c.x;
        if(d.x==0.0) 
        {
                d.ix=0;
                return d;
        }
        d.ix=ix-c.ix;
        if(fabs(d.x)>=BIGS)
        {
                d.x=d.x*BIGI;
                d.ix=d.ix+1;
        }
        else if(fabs(d.x)<BIGSI)
        {
                d.x=d.x*BIG;
                d.ix=d.ix-1;
        }
        d=d.Norm (d);
        return d;
}

XNumber XNumber::sqrt2(XNumber c)
{
        XNumber d;
        d.x=sqrt(c.x);
        d.ix=c.ix/2;
        if(c.ix%2==-1)
        {
                d.x*=BIGSI;
        }
        else if(c.ix%2==1)
        {
                d.x*=BIGS;
        }
        d=d.Norm (d);
        return d;
}

XNumber XNumber::Norm(XNumber c)
{        
        double w;
        w=fabs(c.x);
        if(w>=BIGS)
        {
                c.x=c.x*BIGI;
                c.ix=c.ix+1;
        }
        else if(w<BIGSI)
        {
                c.x=c.x*BIG;
                c.ix=c.ix-1;
        }
        return c;
}

double XNumber::x2f(XNumber c)
{
        double xf;
        if(c.ix==0)
        {
                xf=c.x;
        }
        else if(c.ix<0)
        {
                xf=c.x*BIGI;
        }
        else
        {
                xf=c.x*BIG;
        }
        return xf;
}

转载本文请联系原作者获取授权，同时请注明本文来自周江存科学网博客。
链接地址：https://wap.sciencenet.cn/blog-432384-610926.html

上一篇：高斯滤波系数的计算
下一篇：震后形变