！壞人必需做好事！

獻給同是用C++做數值分析的研究生們！^^ 科科？

程式庫載點：http://ppt.cc/ENG9
程式載點：http://ppt.cc/ZVT-

從下面的結果可以看到一般 double 所不能做的運算，mpfr都能輕易做到。

計算

I(n) = 0積到1 {(x^n)/(x+10)} dx, n=0~30

積法一：
I(0) = ln(11/10)
I(n) = 1/n-10 * I(n-1)
normal c
x 0=    0.09531017980432493500000000000000
x 1=    0.04689820195675065100000000000000
x 2=    0.03101798043249348600000000000000
x 3=    0.02315352900839845500000000000000
x 4=    0.01846470991601545400000000000000
x 5=    0.01535290083984547400000000000000
x 6=    0.01313765826821192100000000000000
x 7=    0.01148056017502363500000000000000
x 8=    0.01019439824976364800000000000000
x 9=    0.009167128613474628800000000000000
x10=    0.008328713865253717400000000000000
x11=    0.007621952256553737900000000000000
x12=    0.0071138107677959500000000000000
x13=    0.005784969245117427300000000000000
x14=    0.01357887897739715200000000000000
x15=    -0.06912212310730485300000000000000
x16=    0.7537212310730485600000000000000
x17=    -7.478388781318720600000000000000
x18=    74.83944336874276400000000000000
x19=    -748.3418021084802400000000000000
x20=    7483.468021084802800000000000000
x21=    -74834.63259180041500000000000000
x22=    748346.3713725495600000000000000
x23=    -7483463.670247235300000000000000
x24=    74834636.74413903100000000000000
x25=    -748346367.4013903100000000000000
x26=    7483463674.052364300000000000000
x27=    -74834636740.48661800000000000000
x28=    748346367404.9019800000000000000
x29=    -7483463674048.985400000000000000
x30=    74834636740489.89100000000000000
mpfr
x 0=    0.095310179804324860043952123280765092220605365308644
x 1=    0.046898201956751399560478767192349077793946346913558
x 2=    0.03101798043248600439521232807650922206053653086442
x 3=    0.023153529008473289381210052568241112727968024689134
x 4=    0.018464709915267106187899474317588872720319753108659
x 5=    0.015352900847328938121005256824111272796802468913415
x 6=    0.013137658193377285456614098425553938698641977532519
x 7=    0.011480560923370002576716158601603470156437367531957
x 8=    0.010194390766299974232838413983965298435626324680429
x 9=    0.0091672034481113687827269712714581267548478643068237
x10=    0.0083279655188863121727302872854187324515213569317632
x11=    0.0076294357202277873636062180549035845756955215914593
x12=    0.0070389761310554596972711527842974875763781174187403
x13=    0.0065333156125223261042115490801020473131419027356735
x14=    0.0060954153033481675293130806275509554400095440718364
x15=    0.0057125136331849913735358603911571122665712259483028
x16=    0.0053748636681500862646413960884288773342877405169724
x17=    0.0050748927302638432359389802921818148924167124773344
x18=    0.0048066282529171231961657526337374066313884307822119
x19=    0.0045652964181971890909740526099943547387472711252497
x20=    0.0043470358180281090902594739000564526125272887475027
x21=    0.004148689438766528145024308618483092922346160144021
x22=    0.0039676510668801730952114592697145253219929440143356
x23=    0.0038017502007634864391897551289417033018096902914264
x24=    0.0036491646590318022747691153772496336485697637524026
x25=    0.0035083534096819772523088462275036635143023624759742
x26=    0.0033780043647186890153730761865018263954379137017964
x27=    0.0032569933898501468833062751720187730826579000190726
x28=    0.0031443518157842454526515339940979834591352855235601
x29=    0.00303924046284720064589845316246844127071611028164
x30=    0.0029409287048613268743488017086489206261722305169328

積法二．
I(30)=0
I(n-1) = 1/(10*n) - 1/10*I(n)

normal c
x30=    0.000000
x29=    0.003333333333333333500000000000000
x28=    0.003114942528735632400000000000000
x27=    0.003259934318555008700000000000000
x26=    0.003377710271848203200000000000000
x25=    0.003508382818969026300000000000000
x24=    0.003649161718103097500000000000000
x23=    0.003801750494856356900000000000000
x22=    0.003967651037470885800000000000000
x21=    0.004148689441707457600000000000000
x20=    0.004347035817734016700000000000000
x19=    0.004565296418226598300000000000000
x18=    0.004806628252914182100000000000000
x17=    0.005074892730264137500000000000000
x16=    0.005374863668150056700000000000000
x15=    0.005712513633184994700000000000000
x14=    0.006095415303348167600000000000000
x13=    0.006533315612522326700000000000000
x12=    0.007038976131055460500000000000000
x11=    0.007629435720227787500000000000000
x10=    0.008327965518886313800000000000000
x 9=    0.009167203448111369700000000000000
x 8=    0.01019439076629997400000000000000
x 7=    0.01148056092337000300000000000000
x 6=    0.01313765819337728600000000000000
x 5=    0.01535290084732893700000000000000
x 4=    0.01846470991526710800000000000000
x 3=    0.02315352900847329100000000000000
x 2=    0.03101798043248600200000000000000
x 1=    0.04689820195675140100000000000000
x 0=    0.09531017980432486500000000000000
mpfr
x30=    0
x29=    0.0033333333333333333333333333333333333333333333333333
x28=    0.0031149425287356321839080459770114942528735632183908
x27=    0.0032599343185550082101806239737274220032840722495895
x26=    0.0033777102718482028826856413063309615033752964787448
x25=    0.0035083828189690258655775897155207500035086241982794
x24=    0.0036491617181030974134422410284479249996491375801721
x23=    0.0038017504948563569253224425638218741667017529086495
x22=    0.0039676510374708860465981905262265082355037377526133
x21=    0.0041486894417074568498856354928318946309950807701932
x20=    0.0043470358177340162197733412126215724416623966848854
x19=    0.0045652964182265983780226658787378427558337603315115
x18=    0.0048066282529141822674608913068630578296797818615857
x17=    0.005074892730264137328809466424869249772587577369397
x16=    0.0053748636681500568553543474751601338462706540277662
x15=    0.0057125136331849943144645652524839866153729345972234
x14=    0.0060954153033481672352202101414182680051293732069443
x13=    0.0065333156125223261336208361287153160566299198221627
x12=    0.0070389761310554596943302240794361607020293157100914
x11=    0.0076294357202277873639003109253897172631304017623242
x10=    0.0083279655188863121727008779983701191827778689146767
x 9=    0.0091672034481113687827299122001629880817222131085323
x 8=    0.010194390766299974232838119891094812302938889800258
x 7=    0.011480560923370002576716188010890518769706111019974
x 6=    0.013137658193377285456614095484625233837315103183717
x 5=    0.015352900847328938121005257118204143282935156348295
x 4=    0.018464709915267106187899474288179585671706484365171
x 3=    0.023153529008473289381210052571182041432829351563483
x 2=    0.031017980432486004395212328076215129190050398176985
x 1=    0.046898201956751399560478767192378487080994960182301
x 0=    0.09531017980432486004395212328076215129190050398177

http://codepad.org/Glbn3yeN

#include <auto_link_mpir.hpp>
#include <auto_link_mpfr.hpp>
#include <iostream>
#include "mpreal.h"

void exp4_mpfr()
{
    using namespace mpfr;
    using namespace std;
    mpreal x = "0";
    mpreal n = "30";
    printf("x%2d=\t", 30);
    cout<<x<<endl;
    for (int i=30;i>0;i--)
    {
        x = mpreal("1")/"10"/n-x/"10";
        n -= "1";
        printf("x%2d=\t", i-1);
        cout << x << endl;
    }
}

void exp4_normal_c()
{
    int n;
    double x;
    x = 0;
    printf("x30=\t%f\n",x);
    for(n=30; n>0; n--)
    {
        x = 1.0/10.0/n - x/10.0;
        printf("x%2d=\t", n-1);
        std::cout << x << std::endl;
    }
}

void exp3_mpfr()
{
    using namespace mpfr;
    using namespace std;
    mpreal x = "11";
    mpreal n = "10";
    x = log(x/n);
    printf("x%2d=\t", 0);
    cout<<x<<endl;
    n = "1";
    for (int i=1;i<31;i++)
    {
        x = "1"/n-"10"*x;
        n += "1";
        printf("x%2d=\t", i);
        cout << x << endl;
    }
}

void exp3_normal_c()
{
    int n;
    double x;
    x = log(11.0/10.0);
    printf("x%2d=\t", 0);
    std::cout << x << std::endl;
    for(n=1; n<31; n++)
    {
        x = 1.0/n - 10*x;
        printf("x%2d=\t", n);
        std::cout << x << std::endl;
    }
}

int main(int argc, char* argv[])
{
    mpfr::mpreal::set_default_prec(1024);
    std::cout.precision(50);
    std::cout << "normal c" << std::endl;
    exp3_normal_c();
    std::cout << "mpfr" << std::endl;
    exp3_mpfr();
    std::cout << "normal c" << std::endl;
    exp4_normal_c();
    std::cout << "mpfr" << std::endl;
    exp4_mpfr();
    return 0;
}