PI的高精度计算

本程序采用LISP计算高精度Pi的数值,可计算到小数点后16000位。

  1. ;;;----------------------------------------------------;
  2. ;;;highflybird    2008.4.20 Haikou                     ;
  3. ;;;采用LISP计算高精度Pi的数值,可计算到小数点后16000位 ;
  4. ;;;测试运行部分                                        ;
  5. ;;;----------------------------------------------------;
  6. (defun c:test (/ digits i n Result t1 t2 time str FILE KEY PATH)
  7.   (setq digits 10000)                                   ;每次输出位数
  8.   (initget 7)                                           ;非零非负非空
  9.   (setq n (getint "nn请输入精度<不超过16000>:"))      ;输入精度
  10.   (if (> n 16000)                               
  11.     (alert "输入超过16000!")                                            
  12.     (progn  
  13.       (setq t1     (getvar "TDUSRTIMER"))               ;开始计时
  14.       (setq Result (CalPi digits n))                    ;计算Pi值
  15.       (setq t2     (getvar "TDUSRTIMER"))               ;计时结束
  16.       (setq time   (* 86400 (- t2 t1)))
  17.       (princ "n计算Pi共用时(秒):")                           
  18.       (princ time)                                      ;和所耗时间
  19.       (initget 1 "Yes No")
  20.       (setq key (getkword "n要打印到文件吗[是(Yes)/否(No)]:"))
  21.       (if (and (= key "Yes")
  22.                (setq path (getfiled "输入文件名:" "c:/" "txt" 1))
  23.           )
  24.         (progn
  25.           (setq file (open path "w"))
  26.           (princ "3n.141" file)                        ;这个是为了打印需要
  27.           (setq i 2)
  28.           (foreach n (cdr Result)
  29.             (setq str (itoa n))
  30.             (while (< (strlen str) 4)
  31.               (setq str (strcat "0" str))
  32.             )
  33.             (princ str file)
  34.             (if (= (rem i 25) 0)
  35.               (princ (strcat "    --" (itoa (* i 4)) "n") file)
  36.             )
  37.             (setq i (1+ i))
  38.           )
  39.           (close file)
  40.           (startapp "Notepad" path)
  41.         )
  42.       )
  43.     )
  44.   )    
  45.   (princ)                                               ;退出
  46. )
  47.  
  48. ;;;----------------------------------------------------;
  49. ;;;highflybird    2008.4.20 Haikou                     ;
  50. ;;;采用LISP计算高精度Pi的数值,可计算到小数点后16000位 ;
  51. ;;;函数求值部分                                        ;
  52. ;;;----------------------------------------------------;
  53. (defun CalPi (digits n / b c d e f g h r s x)
  54.   (setq c (/ (1+ n) (/ (log 2) (log 10))))              ;需要迭代的次数
  55.   (setq c (fix c))                                      ;转化为整数
  56.   (setq e 0 r nil)                                      ;存储结果的字符串赋空值
  57.   (setq h (/ digits 5))                                 ;从小数后算起
  58.   (repeat c                                     
  59.     (setq f (cons h f))                                 ;初始余数为10000 * 2 / 10
  60.   )
  61.   (repeat (1+ (/ n 4))                                  ;重复1+ 800/4 = 201次
  62.     (setq d 0)                                          ;每次末位小数为0
  63.     (setq g (+ c c))                                    ;分母。因为每次循环都输出了4位,所以在后面运算时乘以了a,所以这里得 -2
  64.     (setq b c)                                          ;分子
  65.     (setq x nil)
  66.     (while (> b 0)
  67.       ;;根据公式,乘以分子
  68.       (setq d (* d b))                                  
  69.       (setq b (1- b))
  70.       (setq d (+ d (* (car f) digits)))                 ;因为每次外循环都输出了4位
  71.       ;;根据公式,除以分母
  72.       (setq f (cdr f))
  73.       (setq g (1- g))
  74.       (setq x (cons (rem d g) x))                       ;带分数的 分子部分
  75.       (setq d (/ d g))                                  ;带分数的 整数部分
  76.       (setq g (1- g))
  77.     )
  78.     (setq f (reverse x))
  79.     (repeat 13                                          
  80.       (setq f (cdr f))
  81.     )
  82.     (setq s (+ e (/ d digits)))                         ;printf("%.4d", e+d/a);
  83.     (setq r (cons s r))                                 ;算出的每一项,注意表的每项如果不足4位要加零补全
  84.     (setq e (rem d digits))                             ;e = d % a;
  85.     (setq c (- c 13))                                   ;精度固定为800位,每输出4位后,相当于精度需求降低了4位,故每次可少算13项
  86.   )
  87.   (reverse r)                                           ;把表项反转
  88. )  
  89.  
  90.  
  91. ;;;----------------------------------------------------;
  92. ;;;highflybird    2008.4.20 Haikou                     ;
  93. ;;;采用LISP计算高精度e的数值,可计算到小数点后16000位  ;
  94. ;;;函数求值部分                                        ;
  95. ;;;----------------------------------------------------;
  96. (defun CalE (a c / b d e f p s x)
  97.   (setq b 0 e 0 P "")
  98.   (repeat (1+ c)
  99.     (setq f (cons a f))
  100.   )
  101.   (while (> c 0)
  102.     (setq d 0)
  103.     (setq b c)
  104.     (setq x nil)
  105.     (while (> b 0)    
  106.       (setq d (+ d (* (car f) a)))
  107.       (setq f (cdr f))
  108.       (setq x (cons (rem d b) x))
  109.       (setq d (/ d b))
  110.       (setq b (1- b))
  111.     )
  112.  
  113.     (setq f (reverse x))
  114.     (repeat 14
  115.       (setq f (cdr f))
  116.     )
  117.     (setq c (- c 14))
  118.     (setq s (+ e (/ d a)))
  119.     (setq s (itoa s))
  120.     (while (< (strlen s) 4)
  121.       (setq s (strcat "0" s))
  122.     )
  123.     (setq P (strcat p s))
  124.     (setq e (rem d a))
  125.   )  
  126.   p
  127. )
  128.  
  129. ;|以下是C的源代码:                                      
  130. --------------------------------------------------------
  131. for(;b-c;)                                              
  132. f[b++]=a;                                               
  133. for(;d=0,g=c;c-=14,printf("%.4d ",e+d/a),e=d%a)         
  134. for(b=c;d+=f[b]*a,f[b]=d%b,d/=b,--b;);                  
  135. --------------------------------------------------------
  136. |;;

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注