大数计算：e^1000/300！

原创已于 2023-10-17 00:17:36 修改 · 442 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#大数计算 #阶乘 #幂运算 #C++

于 2023-10-13 17:33:34 首次发布

分析基础专栏收录该内容

7 篇文章

订阅专栏

文章讨论了大数计算中可能遇到的数据类型限制问题，提出使用科学计数法解决，并给出了一个基于模板的BigNumeric类实现，包括乘法、除法、加法、减法、开方和幂运算的代码示例。

1.问题：大数计算可能超出数据类型范围

当单独计算 $e^{1000}$ ，因为 $e^{1000} > e^{700} \approx 1.01432e304$ ，double的最大取值为1.79769e+308，所以 $e^{1000}$ 肯定超过了double的表示范围。

同样，对于300！也是如此。

那我们应该怎么去计算和存储结果呢？

2.解决方案

从数学角度出发，对于超级大的数，运算方式、运算规律等肯定保持不变的。

很多时候，我们主要是利用相关的定理、公式进行简化或者极限处理。

由于我项目里的精度要求不是很高，于是，可以采用相对宽松的方式解决这个问题：

科学计数法！

3.代码实现

/*
* @file     BigNumeric.hpp
* @brief    模板类，用于大数计算，以科学计数法为实现形式，尾数部分有效位数12，指数部分完全准确.
* 直接在工程中引入该hpp文件即可使用该类，重载了+-/*以及幂运算和开方等.
* 由于指数部分采用的市 int ，故限制了本类使用时科学计数法的最大表示范围为：
*                       a*10^2147483647 (-1.0<=a<=1.0)         
**
* @author   DongFengChuiLiu
* @email    308021782@qq.com 
* @version  1.0 
* @date     2023/10/15
*/


#ifndef __BigNumeric_h
#define __BigNumeric_h
#include <cmath>
#include <limits>


template<typename _T>
class BigNumeric
{
private:
    _T mantissa;        //尾数
    int exponent;       //指数

public:
    BigNumeric(_T num)
    {
        if (num == (_T)0)
        {
            exponent = 0;
            mantissa = 0;
        }
        else
        {
            exponent = std::floor(std::log10(std::abs(num)));
            mantissa = num / std::pow(10.0, exponent);
        }
    }

    BigNumeric(_T mnt,int exp)
    {
        this->mantissa = mnt;
        this->exponent = exp;
    }

    ~BigNumeric() {}

    BigNumeric(const BigNumeric& other)
    {
        if (this == &other)
            return;
        this->exponent = other.exponent;
        this->mantissa = other.mantissa;
    }

    BigNumeric& operator=(const BigNumeric& other)
    {
        if(this == &other)
            return *this;
        this->exponent = other.exponent;
        this->mantissa = other.mantissa;
        return *this;
    }

public:
    //取 BigNumeric 对象的实际值
    _T value()
    {
        if (this->exponent >= 308) return std::numeric_limits<_T>::infinity();
        return this->mantissa * std::pow(10.0, this->exponent);
    }


    //乘法运算
    // BigNumeric * BigNumeric
    BigNumeric operator*(const BigNumeric opr)
    {
        BigNumeric<_T> resmnt(this->mantissa * opr.mantissa);
        resmnt.exponent = CheckRange((_T)resmnt.exponent + this->exponent + opr.exponent);
        return resmnt;
    }
    // BigNumeric * _T
    BigNumeric operator*(const _T opr)
    {
        return *this * BigNumeric<_T>(opr);
    }
    //_T * BigNumeric
    friend    BigNumeric operator*(const _T opr, BigNumeric opr1)
    {
        return opr1 * opr;
    }


    //除法
    // BigNumeric / BigNumeric
    BigNumeric operator/(const BigNumeric opr)
    {
        BigNumeric<_T> resmnt(this->mantissa / opr.mantissa);
        resmnt.exponent = CheckRange((_T)resmnt.exponent + this->exponent - opr.exponent);
        return resmnt;
    }
    // BigNumeric / _T
    BigNumeric operator/(const _T opr)
    {
        return *this / BigNumeric<_T>(opr);
    }
    // _T / BigNumeric
    friend     BigNumeric operator/(const _T opr, const BigNumeric opr1)
    {
        return BigNumeric<_T>(opr) / opr1;
    }


    //加法
    // BigNumeric + BigNumeric
    BigNumeric operator+(const BigNumeric opr)
    {
        if (this->exponent - opr.exponent > 15) return *this;
        else if (this->exponent - opr.exponent < -15) return opr;
        int min = this->exponent > opr.exponent ? opr.exponent : this->exponent;
        BigNumeric<_T> resmnt(this->mantissa * std::pow(10.0, this->exponent - min) + opr.mantissa * std::pow(10.0, opr.exponent - min));
        resmnt.exponent = CheckRange((_T)resmnt.exponent + min);
        return resmnt;
    }
    // BigNumeric + _T
    BigNumeric operator+(const _T opr)
    {
        return *this + BigNumeric<_T>(opr);
    }
    // _T + BigNumeric
    friend    BigNumeric operator+(const _T opr, BigNumeric opr1)
    {
        return opr1 + opr;
    }


    //减法
    // BigNumeric - BigNumeric
    BigNumeric operator-(const BigNumeric opr)
    {
        return *this + BigNumeric<_T>(opr) * (-1.0);
    }
    // BigNumeric - _T
    BigNumeric operator-(const _T opr)
    {
        return *this - BigNumeric<_T>(opr);
    }
    // _T - BigNumeric
    friend    BigNumeric operator-(const _T opr, BigNumeric opr1)
    {
        return opr1 - opr;
    }


    //开方
    // BigNumeric^0.5
    BigNumeric Sqrt()
    {
        BigNumeric<_T> result(0.0);
        _T bgnmant = std::sqrt(this->mantissa);
        int bgnexp = this->exponent;
        if (bgnexp % 2 == 0)
        {
            result.mantissa = bgnmant;
            result.exponent = bgnexp / 2;
        }
        else
        {
            BigNumeric temp(bgnmant * std::sqrt(10.0));
            result.mantissa = temp.mantissa;
            result.exponent = temp.exponent + bgnexp / 2;
        }
        return result;
    }


    //幂
    // BigNumeric^_T
    BigNumeric Pow(_T exp)
    {
        BigNumeric temp(Vpow(this->mantissa, exp) * Vpow(10.0, this->exponent * (exp - (int)exp)));
        temp.exponent = CheckRange((_T)temp.exponent + (_T)this->exponent * exp);
        return temp;
    }

public:
    //整数阶乘 int!
    static BigNumeric Factorial(const int opr)
    {
        if (opr < 0) return 1.0 / Factorial(-1.0 * opr + 1);
        else if (opr == 0) return BigNumeric(1.0);
        return Factorial(opr - 1) * opr;
    }


    //e指数幂运算 e^_T
    static BigNumeric Epow(const _T exp)
    {
        BigNumeric res(1.0);
        double e = 2.71828182845904523536;
        if (exp <= 700) return BigNumeric(std::pow(e, exp));
        int count = exp / 700;
        BigNumeric bgn(std::pow(e, 700));
        for (size_t i = 0; i < count; i++)
            res = res * bgn;
        BigNumeric bgn1(std::pow(e, exp - count * 700));
        res = res * bgn1;
        return res;
    }
    //e^BigNumeric
    static BigNumeric Epow(const BigNumeric opr)
    {
        BigNumeric res(1.0);
        double e = 2.71828182845904523536;
        BigNumeric mnt(std::pow(e, opr.mantissa));
        int count = std::abs(opr.exponent) > 1000 ? 1000 : std::abs(opr.exponent);
        if (count == 0) return mnt;
        if(opr.exponent < 0) res = mnt.Pow(0.1);
        else res = mnt.Pow(10.0);
        for (size_t i = 0; i < count - 1; i++)
        {
            if (opr.exponent < 0) res = res.Pow(0.1);
            else res = res.Pow(10.0);
        }        
        return res;
    }


    //任意底数的幂运算 a^_T
    static BigNumeric Vpow(const _T e, const _T exp)
    {
        BigNumeric res(1.0);
        BigNumeric bgnmant(e);
        int chk = bgnmant.exponent == 0 ? exp : CheckRange(exp * bgnmant.exponent);
        if (chk <= 300) return BigNumeric(std::pow(e, exp));
        int count = exp / 300;
        BigNumeric bgn(std::pow(bgnmant.mantissa, 300));
        for (size_t i = 0; i < count; i++)
        {
            res = res * bgn;
            res.exponent = CheckRange(res.exponent + (_T)bgnmant.exponent * 300);
        }
        BigNumeric bgn1(std::pow(bgnmant.mantissa, exp - count * 300));
        res = res * bgn1;
        res.exponent = CheckRange(res.exponent + (_T)bgnmant.exponent * (exp - count * 300));
        return res;
    }

private:
    //控制指数的值范围
    static int CheckRange(const _T num)
    {
        if (num >= (_T)2147483647) return 2147483647;
        else if( num <= (_T)-2147483648) return -2147483648;
        return num;
    }
};

#endif // !__BigNumeric_h

4.测试

//测试代码
#include "BigNumeric.hpp"


int main() 
{    
    // 乘法测试
    /*BigNumeric<double> bignum(2.13e303);
    BigNumeric<double> bignum1(5e303);
    BigNumeric<double> res = bignum * bignum1;
    res = bignum * 2;
    res = 2 * bignum;
    double t = res.value();*/

    // 除法测试
    /*BigNumeric<double> bignum(2.13e303);
    BigNumeric<double> bignum1(4.26e303);  
    BigNumeric<double> res = bignum1 / bignum;
    res = 240.0 / bignum;
    res = bignum / 2.13e302;*/

    //加法测试
    /*BigNumeric<double> bignum(987654321);
    BigNumeric<double> bignum1(9);
    BigNumeric<double> res(bignum + bignum1);
    res = bignum + 9;
    res = 9 + bignum;*/


    //减法测试
    /*BigNumeric<double> bignum(987654321);
    BigNumeric<double> bignum1(9);
    BigNumeric<double> res(bignum - bignum1);
    res = bignum - 9;
    res = 987654321 - bignum;*/


    //幂测试
    /*BigNumeric<double> bignum(2.15 * 1.0e303);
    BigNumeric<double> TMP = (bignum*54321).Sqrt();
    BigNumeric<double> TMP1 = (bignum * 54321).Pow(0.5);
        
    bignum  = BigNumeric<double>(9.15);
    BigNumeric<double> TMP2 = bignum.Pow(400);
    double t = TMP2.value();*/

    //阶乘测试
    BigNumeric<double> bignum = BigNumeric<double>::Factorial(300);
    BigNumeric<double> bignum1 = BigNumeric<double>::Epow(1000);
    bignum = bignum1 / bignum;
    double t = bignum.value();

    return 0;
}

结果：6.4369310844548986e-181，尾数部分有效位数为 12，指数部分完全准确。