开心一刻，考考你的计算能力

G-Spider

7# mathe
我用matlab算了一下，最后一行数据不一样。
---------------------------------------------------------------------
。。。
3799248216530963946104908731687 718950889

sheng_jianguo

7# mathe
我用matlab算了一下，最后一行数据不一样。
---------------------------------------------------------------------
。。。
3799248216530963946104908731687 718950889
G-Spider 发表于 2011-1-13 19:30

用matlab算，怎样才能得到如此高的精度？（我也用matlab算，好像精度不够）

G-Spider

12# sheng_jianguo
>>help vpa
VPA    Variable precision arithmetic.
    R = VPA(S) numerically evaluates each element of the double matrix
    S using variable precision floating point arithmetic with D decimal
    digit accuracy, where D is the current setting of DIGITS.
    The resulting R is a SYM.
  
    VPA(S,D) uses D digits, instead of the current setting of DIGITS.
    D is an integer or the SYM representation of a number.

    It is important to avoid the evaluation of an expression using double
    precision floating point arithmetic before it is passed to VPA.
    For example,
       phi = vpa((1+sqrt(5))/2)
    first computes a 16-digit approximation to the golden ratio, then
    converts that approximation to one with d digits, where d is the current
    setting of DIGITS.  To get full precision, use unevaluated string or
    symbolic arguments,
       phi = vpa('(1+sqrt(5))/2')
    or
       s = sym('sqrt(5)')
       phi = vpa((1+s)/2);

//////////////////////////
可见matlab有良好的大数运算功能。
vpa('173746*sin(10^22)+94228*log(17.1)-78487*exp(0.42)',1000)  %注意要加上引号，不然每一项不到15位的精度，最后的结果比1000的精度略低一点，存在舍入

sheng_jianguo

13# G-Spider
谢谢！

winxos

我是打酱油的，
过来看看各位老大。

wayne

9# wayne
代码不好，在windows的MinGW下编译，即使开了-Wall也不警告，运行也很正常。
在linux下编译没问题，但运行时会出现segmentation fault .
检查了很长时间，原来是argv[1] 使然。
现在更正如下：

1.    #include <stdio.h>
   
2.    #include <stdlib.h>
   
3.    #include <mpfr.h>
   
4.    int main (int argc, char *argv[])
   
5.    {int m;
   
6.    if (argc>1) m= atoi(argv[1])-1;
   
7.            else printf("usage: %s n\n",argv[0]),exit(1);
   
8.    mp_prec_t p = 100+5*m;
   
9.    mpfr_t a, b, c, d;
   
10.    mpfr_inits2 (p, a, b, c, d, (mpfr_ptr) 0);
    
11.    mpfr_set_str (a, "1e22", 10, GMP_RNDN);
    
12.    mpfr_sin (a, a, GMP_RNDN);
    
13.    mpfr_mul_ui (a, a, 173746, GMP_RNDN);
    
14.    mpfr_set_str (b, "17.1", 10, GMP_RNDN);
    
15.    mpfr_log (b, b, GMP_RNDN);
    
16.    mpfr_mul_ui (b, b, 94228, GMP_RNDN);
    
17.    mpfr_set_str (c, "0.42", 10, GMP_RNDN);
    
18.    mpfr_exp (c, c, GMP_RNDN);
    
19.    mpfr_mul_si (c, c, -78487, GMP_RNDN);
    
20.    mpfr_add (d, a, b, GMP_RNDN);
    
21.    mpfr_add (d, d, c, GMP_RNDN);
    
22.    mpfr_printf ("d=%.*Re\n",m,d);
    
23.    mpfr_clears (a, b, c, d, NULL);
    
24.    return 0;
    
25.    }

复制代码