圆周率小数后800位！

Highflybird

下面给出一个程序，纯lisp的，能把圆周率精确到小数后800位。稍加修改就可以到10000位，当然远没有C语言快。
不过程序运行命令 test,估计10秒中左右就能把这个数打印出来。
[pcode=lisp,true]
;;;highflybird    2008.4.20 Haikou
(defun c:test (/ a c p t1 t2 time)
  (setq a 10000 c 2800)
  (setq t1 (getvar "TDUSRTIMER"))
  (setq p  (CalPi a c))
  (setq t2 (getvar "TDUSRTIMER"))
  (setq time (* 86400 (- t2 t1)))
  (setq p (strcat "3." (substr p 2)))
  (princ p)
  (princ "\n共耗时<秒>:")
  (princ time)
  (princ)
)
;;;计算函数
(defun CalPi (a c / b d e f g h p s x)
  (setq b 0 e 0 P "")
  (setq h (/ a 5))
  (repeat (1+ c)
    (setq f (cons h f))
  )
  (while (> c 0)
    (setq d 0)
    (setq g (+ c c))
    (setq b c)
    (setq x nil)
    (while (> b 0)
      (setq d (* d b))
      (setq b (1- b))
      (setq d (+ d (* (car f) a)))
      (setq f (cdr f))
      (setq g (1- g))
      (setq x (cons (rem d g) x))
      (setq d (/ d g))
      (setq g (1- g))
    )
    (setq f (reverse x))
    (repeat 14
      (setq f (cdr f))
    )
    (setq s (+ e (/ d a)))
    (setq s (itoa s))
    (while (< (strlen s) 4)
      (setq s (strcat "0" s))
    )
    (setq P (strcat p s))
    (setq e (rem d a))
    (setq c (- c 14))
  )
  p
)  
;;;以后给出完全注释版的
[/pcode]

今天重新优化了这个程序，使得速度提高了不少，而且编译成vlx文件，速度更是大步提升。
你若对此有趣，不妨对编译前后做一下比较。
在我的机器上，运行它不到0.5秒。
lisp文件：

请点击此处下载

查看状态：需购买或无权限

您的用户组是：游客

文件名称:Pi.lsp 
下载次数:15  文件大小:1.09 KB 
下载权限: 不限 以上  [免费赚D豆]

vlx文件：

请点击此处下载

查看状态：需购买或无权限

您的用户组是：游客

文件名称:PI.rar 
下载次数:1  文件大小:1.05 KB 
下载权限: 不限 以上  [免费赚D豆]

原理介绍参见这篇帖子：http://blog.csdn.net/zyl910/article/details/1368387

详细注视的lisp文件在这里：

请点击此处下载

查看状态：需购买或无权限

您的用户组是：游客

文件名称:Pi高精度.lsp 
下载次数:17  文件大小:3.54 KB 
下载权限: 不限 以上  [免费赚D豆]

部分源代码：
[pcode=lisp,true];;;高精度计算Pi函数
(defun CalPi (digits n / b c d e f g h r s x)
  (setq c (/ (1+ n) (/ (log 2) (log 10))))      ;需要迭代的次数
  (setq c (fix c))                          ;转化为整数
  (setq  e 0 r nil)                                 ;存储结果的字符串赋空值
  (setq h (/ digits 5))                           ;从小数后算起
  (repeat c                                    
    (setq f (cons h f))                           ;初始余数为10000 * 2 / 10
  )
  (repeat (1+ (/ n 4))                            ;重复1+ 800/4 = 201次
    (setq d 0)                                    ;每次末位小数为0
    (setq g (+ c c))                              ;分母。因为每次循环都输出了4位，所以在后面运算时乘以了a，所以这里得 -2
    (setq b c)            ;分子
    (setq x nil)
    (while (> b 0)
      ;;根据公式，乘以分子
      (setq d (* d b))         
      (setq b (1- b))
      (setq d (+ d (* (car f) digits)))      ;因为每次外循环都输出了4位
      ;;根据公式，除以分母
      (setq f (cdr f))
      (setq g (1- g))
      (setq x (cons (rem d g) x))      ;带分数的 分子部分
      (setq d (/ d g))          ;带分数的 整数部分
      (setq g (1- g))
    )
    (setq f (reverse x))
    (repeat 13            
      (setq f (cdr f))
    )
    (setq s (+ e (/ d digits)))        ;printf("%.4d", e+d/a);
    (setq r (cons s r))          ;算出的每一项，注意表的每项如果不足4位要加零补全
    (setq e (rem d digits))        ;e = d % a;
    (setq c (- c 13))          ;因为精度固定为800位，每输出4位后，相当于精度需求降低了4位，所以每次可以少算13项
  )
  (reverse r)            ;把表项反转
)  
[/pcode]

yangliucq

虽然看不懂，但是肯定是好贴！

sachindkini

dear sir,

thx for sharing