C++ BigFloat类代码示例

int main(void)
{
    // 

#ifdef T_addBigInt
    string num1("123");
    string num2("4567");

    string sum = addBigInt(num1, num2);

    if (sum == "")
    {
        cerr << "error in addBigInt" << endl;
        return -1;
    }

    cout << num1 << " + " << num2 << " = " << sum << endl;
#endif 

#ifdef T_mulOneBit
    BigFloat bigNum;
    cout << "enter a big float: ";
    cin >> bigNum;
    cout << "digits = " << bigNum.digits << endl;
    cout << "exp = " << bigNum.exp << endl;

    string resu(bigNum.mulOneBit('0'));

    if ("" == resu)
    {
        cerr << "mulOneBit error"  << endl;
        return -1;
    }
    cout << "result is " << resu << endl;
#endif

#ifdef T_opMul

    BigFloat num1(100);
    BigFloat num2;
    cout << "Enter a number : ";
    cin >> num2;
    BigFloat prdt = num1 * num2;

    cout << "result is : \n" << prdt << endl;
#endif

    return 0;
}

int BigFloat::ucmp(const BigFloat &x) const{
    //  Compare function that ignores the sign.
    //  This is needed to determine which direction subtractions will go.

    //  Magnitude
    int64_t magA = exp + L;
    int64_t magB = x.exp + x.L;
    if (magA > magB)
        return 1;
    if (magA < magB)
        return -1;

    //  Compare
    int64_t mag = magA;
    while (mag >= exp || mag >= x.exp){
        uint32_t wordA = word_at(mag);
        uint32_t wordB = x.word_at(mag);
        if (wordA < wordB)
            return -1;
        if (wordA > wordB)
            return 1;
        mag--;
    }
    return 0;
}

BigFloat BigFloat::usub(const BigFloat &x,size_t p) const{
    //  Perform subtraction ignoring the sign of the two operands.

    //  "this" must be greater than or equal to x. Otherwise, the behavior
    //  is undefined.

    //  Magnitude
    int64_t magA = exp + L;
    int64_t magB = x.exp + x.L;
    int64_t top = std::max(magA,magB);
    int64_t bot = std::min(exp,x.exp);

    //  Truncate precision
    int64_t TL = top - bot;

    if (p == 0){
        //  Default value. No trunction.
        p = (size_t)TL;
    }else{
        //  Increase precision
        p += YCL_BIGFLOAT_EXTRA_PRECISION;
    }

    if (TL > (int64_t)p){
        bot = top - p;
        TL = p;
    }

    //  Compute basic fields.
    BigFloat z;
    z.sign  = sign;
    z.exp   = bot;
    z.L     = (uint32_t)TL;

    //  Allocate mantissa
    z.T = std::unique_ptr<uint32_t[]>(new uint32_t[z.L]);

    //  Subtract
    int32_t carry = 0;
    for (size_t c = 0; bot < top; bot++, c++){
        int32_t word = (int32_t)word_at(bot) - (int32_t)x.word_at(bot) - carry;
        carry = 0;
        if (word < 0){
            word += 1000000000;
            carry = 1;
        }
        z.T[c] = word;
    }

    //  Strip leading zeros
    while (z.L > 0 && z.T[z.L - 1] == 0)
        z.L--;
    if (z.L == 0){
        z.exp = 0;
        z.sign = true;
        z.T.reset();
    }

    return z;
}

inline BigFloat Epsilon<BigFloat>( const BigFloat& alpha )
{
    // NOTE: This 'precision' is the number of bits in the mantissa
    auto p = alpha.Precision();

    // epsilon = b^{-(p-1)} / 2 = 2^{-p}
    BigFloat one(1,p);
    return one >> p;
}

inline BigFloat SafeMin( const BigFloat& alpha )
{
    // TODO: Decide how to only recompute this when the precision changes
    const BigFloat one = BigFloat(1,alpha.Precision());
    const BigFloat eps = Epsilon(alpha);
    const BigFloat minVal = Min(alpha);
    const BigFloat invMax = one/Max(alpha);
    return ( invMax>minVal ? invMax*(one+eps) : minVal );
}

BigFloat invsqrt(uint32_t x,size_t p,int tds){
    //  Compute inverse square root using Newton's Method.

    //            (  r0^2 * x - 1  )
    //  r1 = r0 - (----------------) * r0
    //            (       2        )

    if (x == 0)
        throw "Divide by Zero";

    //  End of recursion. Generate starting point.
    if (p == 0){
        double val = 1. / sqrt((double)x);

        int64_t exponent = 0;

        //  Scale
        while (val < 1000000000.){
            val *= 1000000000.;
            exponent--;
        }

        //  Rebuild a BigFloat.
        uint64_t val64 = (uint64_t)val;

        BigFloat out;
        out.sign = true;

        out.T = std::unique_ptr<uint32_t[]>(new uint32_t[2]);
        out.T[0] = (uint32_t)(val64 % 1000000000);
        out.T[1] = (uint32_t)(val64 / 1000000000);
        out.L = 2;
        out.exp = exponent;

        return out;
    }

    //  Half the precision
    size_t s = p / 2 + 1;
    if (p == 1) s = 0;
    if (p == 2) s = 1;

    //  Recurse at half the precision
    BigFloat T = invsqrt(x,s,tds);

    BigFloat temp = T.mul(T,p);     //  r0^2
    temp = temp.mul(x,p,tds);       //  r0^2 * x
    temp = temp.sub(BigFloat(1),p); //  r0^2 * x - 1
    temp = temp.mul(500000000);         //  (r0^2 * x - 1) / 2
    temp.exp--;
    temp = temp.mul(T,p,tds);       //  (r0^2 * x - 1) / 2 * r0
    return T.sub(temp,p);           //  r0 - (r0^2 * x - 1) / 2 * r0
}

BigFloat BigFloat::sub(const BigFloat &x,size_t p) const{
    //  Subtraction

    //  The target precision is p.
    //  If (p = 0), then no truncation is done. The entire operation is done
    //  at maximum precision with no data loss.

    //  Different sign. Add.
    if (sign != x.sign)
        return uadd(x,p);

    //  this > x
    if (ucmp(x) > 0)
        return usub(x,p);

    //  this < x
    BigFloat z = x.usub(*this,p);
    z.negate();
    return z;
}

int main()
{
	BigFloat *sum;
	sum = new BigFloat(0.0);
	
	for ( int i = 1; i <= 100 ; i ++)
	{
		BigFloat index(1.0);
		BigFloat temp((double)i);
		if ( 0 == index.divide(&index, &temp))
			sum->add(sum,&index);
		else
		{
			printf("error");
			exit(1);
		}
	}
	sum->printAll();
	sum->print_fraction_portion_in_decimal(1000,NULL);
}

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
//  Pi
void Pi_BSR(BigFloat &P,BigFloat &Q,BigFloat &R,uint32_t a,uint32_t b,size_t p){
    //  Binary Splitting recursion for the Chudnovsky Formula.

    if (b - a == 1){
        //  P = (13591409 + 545140134 b)(2b-1)(6b-5)(6b-1) (-1)^b
        P = BigFloat(b).mul(545140134);
        P = P.add(BigFloat(13591409));
        P = P.mul(2*b - 1);
        P = P.mul(6*b - 5);
        P = P.mul(6*b - 1);
        if (b % 2 == 1)
            P.negate();

        //  Q = 10939058860032000 * b^3
        Q = BigFloat(b);
        Q = Q.mul(Q).mul(Q).mul(26726400).mul(409297880);

        //  R = (2b-1)(6b-5)(6b-1)
        R = BigFloat(2*b - 1);
        R = R.mul(6*b - 5);
        R = R.mul(6*b - 1);

        return;
    }

    uint32_t m = (a + b) / 2;

    BigFloat P0,Q0,R0,P1,Q1,R1;
    Pi_BSR(P0,Q0,R0,a,m,p);
    Pi_BSR(P1,Q1,R1,m,b,p);

    P = P0.mul(Q1,p).add(P1.mul(R0,p),p);
    Q = Q0.mul(Q1,p);
    R = R0.mul(R1,p);
}

bool GaussElimination::changeValuesInRowByBigFloat(vector<vector<BigFloat>>& matrix, int indexRow)
{
	if( matrix[indexRow][indexRow].compare( &BigFloat(0.0) ) != 0 ) return true;

	// wenn i.zeile zahl in spalte i = 0, tausche mit zeile darunter!!, wo zahl in i.spalte != 0 -> falls nicht da, gs nicht lösbar
	for(unsigned int rowCounter = indexRow+1; rowCounter < row; ++rowCounter)
	{
		BigFloat value = matrix[rowCounter][indexRow];
		if(value.compare( &BigFloat(0.0) ) != 0 )
		{
			// tausche reihen
			for(unsigned int i = 0; i <= col; ++i) //zusätzliche spalte angefügt durch vector
			{
				BigFloat currowValue = matrix[indexRow][i];
				BigFloat changeValue = matrix[rowCounter][i];
				matrix[indexRow][i] = changeValue;
				matrix[rowCounter][i] = currowValue;
			}
			return true;
		}
	}

	return false;
}

	// Zoom vers la cible du clic
	static void clicGauche(SDL_Event& event, Affichage* disp)
	{

		// Calcul de la position du clic dans le plan complexe
		float x = disp->center.x + disp->scale*(((float)event.motion.x) / WIDTH - 0.5f);
		float y = disp->center.y + disp->scale*(((float)event.motion.y) / HEIGHT - 0.5f);

		// MAJ de la position du nouveau centre dans le plan complexe
		if (INTERACTIVE) {
			disp->center.x = x + (disp->center.x - x) * ZOOM_FACTOR;
			disp->center.y = y + (disp->center.y - y) * ZOOM_FACTOR;
			updateBigCenter(event, true);
		}
		

		// MAJ de l'echelle
		disp->scale *= ZOOM_FACTOR;
		//BigFloat zoomFactor(true, 0, 0x80000000, 0, 0);
		BigFloat temp, temp2;
		BigFloat::mult(ZOOM_FACTOR, *bigScale, temp);
		bigScale->reset();
		temp2.copy(*bigScale);
		BigFloat::add(temp, temp2, *bigScale);
		
		//Nouveau calcul de la fractale avec chrono
		disp->start = chrono::system_clock::now();

		if (GPU && BIG_FLOAT_SIZE == 0)
			affichageGPU(disp);
		else if (GPU)
			computeBigMandelGPU(disp, xCenter->pos, xCenter->decimals, yCenter->pos, yCenter->decimals, bigScale->decimals);
		else
			if (BIG_FLOAT_SIZE == 0)
				computeMandel(disp->pixels, disp->center, disp->scale);
			else {
				computeBigMandel(disp->pixels, *xCenter, *yCenter, *bigScale);
			}

		disp->end = chrono::system_clock::now();
		disp->duration = disp->end - disp->start;
		cout << "Frame computing time : " << disp->duration.count() << endl;
		//cout << "Frame computing scale : " << disp->scale << endl;
		
		// Affichage de la fractale
		disp->dessin();
	}

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
//  Pi
void Pi_BSR(BigFloat &P,BigFloat &Q,BigFloat &R,uint32_t a,uint32_t b,size_t p,int tds = 1){
    //  Binary Splitting recursion for the Chudnovsky Formula.

    if (b - a == 1){
        //  P = (13591409 + 545140134 b)(2b-1)(6b-5)(6b-1) (-1)^b
        P = BigFloat(b).mul(545140134);
        P = P.add(BigFloat(13591409));
        P = P.mul(2*b - 1);
        P = P.mul(6*b - 5);
        P = P.mul(6*b - 1);
        if (b % 2 == 1)
            P.negate();

        //  Q = 10939058860032000 * b^3
        Q = BigFloat(b);
        Q = Q.mul(Q).mul(Q).mul(26726400).mul(409297880);

        //  R = (2b-1)(6b-5)(6b-1)
        R = BigFloat(2*b - 1);
        R = R.mul(6*b - 5);
        R = R.mul(6*b - 1);

        return;
    }

    uint32_t m = (a + b) / 2;

    BigFloat P0,Q0,R0,P1,Q1,R1;

    if (b - a < 1000 || tds < 2){
        //  No more threads.
        Pi_BSR(P0,Q0,R0,a,m,p);
        Pi_BSR(P1,Q1,R1,m,b,p);
    }else{
        //  Run sub-recursions in parallel.
        int tds0 = tds / 2;
        int tds1 = tds - tds0;
#pragma omp parallel num_threads(2)
        {
            int tid = omp_get_thread_num();
            if (tid == 0){
                Pi_BSR(P0,Q0,R0,a,m,p,tds0);
            }
            if (tid != 0 || omp_get_num_threads() < 2){
                Pi_BSR(P1,Q1,R1,m,b,p,tds1);
            }
        }
    }

    P = P0.mul(Q1,p,tds).add(P1.mul(R0,p,tds),p);
    Q = Q0.mul(Q1,p,tds);
    R = R0.mul(R1,p,tds);
}

BigFloat BigFloat::add(const BigFloat &x,size_t p) const{
    //  Addition

    //  The target precision is p.
    //  If (p = 0), then no truncation is done. The entire operation is done
    //  at maximum precision with no data loss.

    //  Same sign. Add.
    if (sign == x.sign)
        return uadd(x,p);

    //  this > x
    if (ucmp(x) > 0)
        return usub(x,p);

    //  this < x
    return x.usub(*this,p);
}

void ComputeFloatSession<wtype>::write_digits(const BigFloat<wtype>& x, const std::string& name, const std::string& algorithm){
    std::string dec_name, hex_name;
    if (algorithm.empty()){
        dec_name = name + " - Dec.txt";
        hex_name = name + " - Hex.txt";
    }else{
        dec_name = name + " - Dec - " + algorithm + ".txt";
        hex_name = name + " - Hex - " + algorithm + ".txt";
    }

    //  Special Case: If the number is close to one, don't count the digits in
    //  front of the decimal place. This keeps the behavior consistent with
    //  y-cruncher and other Pi programs that start counting digits after the
    //  decimal place.
    upL_t dec_digits = m_dec;
    upL_t hex_digits = m_hex;
    if (x.get_mag() == 1){
        wtype top_word = x[0];
        dec_digits += std::to_string(top_word).size();
        while (top_word != 0){
            top_word >>= 4;
            hex_digits++;
        }
    }

BigFloat BigFloat::rcp(size_t p,int tds) const{
    //  Compute reciprocal using Newton's Method.

    //  r1 = r0 - (r0 * x - 1) * r0

    if (L == 0)
        throw "Divide by Zero";

    //  Collect operand
    int64_t Aexp = exp;
    size_t AL = L;
    uint32_t *AT = T.get();

    //  End of recursion. Generate starting point.
    if (p == 0){
        //  Truncate precision to 3.
        p = 3;
        if (AL > p){
            size_t chop = AL - p;
            AL = p;
            Aexp += chop;
            AT += chop;
        }

        //  Convert number to floating-point.
        double val = AT[0];
        if (AL >= 2)
            val += AT[1] * 1000000000.;
        if (AL >= 3)
            val += AT[2] * 1000000000000000000.;

        //  Compute reciprocal.
        val = 1. / val;
        Aexp = -Aexp;

        //  Scale
        while (val < 1000000000.){
            val *= 1000000000.;
            Aexp--;
        }

        //  Rebuild a BigFloat.
        uint64_t val64 = (uint64_t)val;

        BigFloat out;
        out.sign = sign;

        out.T = std::unique_ptr<uint32_t[]>(new uint32_t[2]);
        out.T[0] = (uint32_t)(val64 % 1000000000);
        out.T[1] = (uint32_t)(val64 / 1000000000);
        out.L = 2;
        out.exp = Aexp;

        return out;
    }

    //  Half the precision
    size_t s = p / 2 + 1;
    if (p == 1) s = 0;
    if (p == 2) s = 1;

    //  Recurse at half the precision
    BigFloat T = rcp(s,tds);

    //  r1 = r0 - (r0 * x - 1) * r0
    return T.sub(this->mul(T,p,tds).sub(BigFloat(1),p).mul(T,p,tds),p);
}