PI的高精度计算

;;;----------------------------------------------------;
;;;highflybird    2008.4.20 Haikou                     ;
;;;采用LISP计算高精度Pi的数值，可计算到小数点后16000位 ;
;;;测试运行部分                                        ;
;;;----------------------------------------------------;
(defun c:test (/ digits i n Result t1 t2 time str FILE KEY PATH)
  (setq digits 10000)                                   ;每次输出位数
  (initget 7)                                           ;非零非负非空
  (setq n (getint "\n\n请输入精度<不超过16000>:"))      ;输入精度
  (if (> n 16000)                               
    (alert "输入超过16000!")      
    (progn  
      (setq t1     (getvar "TDUSRTIMER"))               ;开始计时
      (setq Result (CalPi digits n))                    ;计算Pi值
      (setq t2     (getvar "TDUSRTIMER"))               ;计时结束
      (setq time   (* 86400 (- t2 t1)))
      (princ "\n计算Pi共用时(秒):")                           
      (princ time)                                      ;和所耗时间
      (initget 1 "Yes No")
      (setq key (getkword "\n要打印到文件吗[是(Yes)/否(No)]："))
      (if (and (= key "Yes")
               (setq path (getfiled "输入文件名:" "c:/" "txt" 1))
          )
        (progn
          (setq file (open path "w"))
          (princ "3\n.141" file)                        ;这个是为了打印需要
          (setq i 2)
          (foreach n (cdr Result)
            (setq str (itoa n))
            (while (< (strlen str) 4)
              (setq str (strcat "0" str))
            )
            (princ str file)
            (if (= (rem i 25) 0)
              (princ (strcat "    --" (itoa (* i 4)) "\n") file)
            )
            (setq i (1+ i))
          )
          (close file)
          (startapp "Notepad" path)
        )
      )
    )
  )    
  (princ)                                               ;退出
)

;;;----------------------------------------------------;
;;;highflybird    2008.4.20 Haikou                     ;
;;;采用LISP计算高精度Pi的数值，可计算到小数点后16000位 ;
;;;函数求值部分                                        ;
;;;----------------------------------------------------;
(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)                                           ;把表项反转
)  
 

;;;----------------------------------------------------;
;;;highflybird    2008.4.20 Haikou                     ;
;;;采用LISP计算高精度e的数值，可计算到小数点后16000位  ;
;;;函数求值部分                                        ;
;;;----------------------------------------------------;
(defun CalE (a c / b d e f p s x)
  (setq b 0 e 0 P "")
  (repeat (1+ c)
    (setq f (cons a f))
  )
  (while (> c 0)
    (setq d 0)
    (setq b c)
    (setq x nil)
    (while (> b 0)    
      (setq d (+ d (* (car f) a)))
      (setq f (cdr f))
      (setq x (cons (rem d b) x))
      (setq d (/ d b))
      (setq b (1- b))
    )
   
    (setq f (reverse x))
    (repeat 14
      (setq f (cdr f))
    )
    (setq c (- c 14))
    (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))
  )
  p
)

;|以下是C的源代码:                                      
--------------------------------------------------------
for(;b-c;)                                              
f[b++]=a;                                               
for(;d=0,g=c;c-=14,printf("%.4d ",e+d/a),e=d%a)         
for(b=c;d+=f[b]*a,f[b]=d%b,d/=b,--b;);                  
--------------------------------------------------------
|;;