[风之影][Lisp大挑战第一季]圆周率

cabinsummer

jicqj 发表于 2011-11-5 19:57 大学时候用c语言玩过圆周率 ， 但是不能和大师比啊

思路是什么？可以在lisp里尝试一下

czb203

看眼界了,希望多来点这种有挑战性的,可以看到更多的大师浮出水面!

cabinsummer

遗憾，计算4000位时竟然出错了

cabinsummer

cabinsummer 发表于 2011-11-11 19:48
遗憾，计算4000位时竟然出错了

找到原因了，不是算法问题，是输出字符串函数的问题。LISP只能处理的字符串最多只能有4096个字符

jicqj

cabinsummer 发表于 2011-11-12 09:45
找到原因了，不是算法问题，是输出字符串函数的问题。LISP只能处理的字符串最多只能有4096个字符

哇  原来如此

VBALISPER

小弟不才,搞个最简单的,还是在朋友们的帮助下完成的.先贴上,慢慢搞一个莱布尼兹,再搞马青...
(defun C:uu();密率π=355/113,π=16arctan1/5-4arctan1/239马青公式
(setq fz 355)
(setq fm 113)
(setq n 0)
(setq lpai nil)
(repeat 30       
        (setq shang (fix (/ fz fm)));求商
          (setq fz (rem fz fm));求余数
        (setq lpai(append lpai (list shang)));第一个参数,每一个除法
        (setq fz (* fz 10))
        (setq n (+ n 1))
);repeat
lpai
)
结果:7位小数精度3. 1 4 1 5 9 2 9 2 0 3 5 3 9 8 2 3 0 0 8 8 4 9 5 5 7 5 2 2 1
是表的形式,小数点是我手加上去的.这两步以后再完善.

qjchen

问题很有趣，前些日子编了一下，始终比Highflybird兄的外星人代码算法慢许多，采用的是Machin算法，大数计算。一直想改进后再发的，不过近日开始忙碌了，就先到此为止，先发下，以后再慢慢改进吧。
希望cabinsummer多出些有趣的题目来做 ：）

1. ;;; By qjchen@gmail.com ，采用Machin算法，不过由于采用了大数计算，速度随n增大急剧下降，
   
2. ;;; 应该是O(n2)的算法，以后再改进吧，大概只能算800位内时间还能接受，更大的就不要了
   
3. ;;; 远比不上Highflybird兄的外星人算法
   
4. ;;; 2011.11.16
   
5. (defun q:math:big+(a b n / f i res x)
   
6.   (setq i (1- n) f 0)
   
7.   (repeat n
   
8.     (if (> (setq x (+ (nth i a) (nth i b) f)) 10000) (setq x (- x 10000) f 1) (setq f 0))
   
9.     (setq res (cons x res) i (1- i))
   
10.   )
    
11.   res
    
12. )
    
13. 
    
14. (defun q:math:big-(a b n / f i res x)
    
15.   (setq i (1- n) f 0)
    
16.   (repeat n
    
17.     (if (< (setq x (- (nth i a) (nth i b) f)) 0) (setq x (+ x 10000) f 1) (setq f 0))
    
18.     (setq res (cons x res) i (1- i))
    
19.   )
    
20.   res
    
21. )
    
22. 
    
23. (defun q:math:round(x)
    
24. (setq y (fix x))
    
25. (if (> (- x y) 0.5) (1+ y) y)
    
26. )
    
27. 
    
28. (defun q:math:bigdiv(a b n / c f i res x)
    
29.   (setq i 0 f 0)
    
30.   (repeat n
    
31.     (setq c (+ (nth i a) (* f 10000))
    
32.           x (fix (/ (+ c 0.1) b))
    
33.           res (cons x res)
    
34.           i (1+ i)
    
35.           f (rem c b))
    
36.   )
    
37.   (reverse res)
    
38. )
    
39. 
    
40. (defun q:math:bigreci(f b n / c i res x)
    
41.   (setq res (cons (fix (/ f b)) res)
    
42.         f (rem f b)
    
43.         i 1)
    
44.   (while (< i n)
    
45.     (setq c (* f 10000)
    
46.           x (fix (/ (+ c 0.1) b))
    
47.           res (cons x res) i (1+ i) f (rem c b))
    
48.   )
    
49.   (reverse res)
    
50. )
    
51. 
    
52. (defun q:math:inin (v n / res)
    
53.   (repeat (1- n)
    
54.     (setq res (cons 0 res))
    
55.   )
    
56.   (cons v res)
    
57. )
    
58. 
    
59. (defun q:pi:deal1(a b / i n1 n2 res)
    
60.   (setq n1 (length a) n2 (length b) res b i 0)
    
61.   (if (< n2 n1)
    
62.      (repeat (- n1 n2)
    
63.        (setq res (append res (list (nth (+ n2 i) a))) i (1+ i))
    
64.      )
    
65.   )
    
66.   res
    
67. )
    
68. 
    
69. (defun q:math:arctg(k v zf N / i n2)
    
70.   (setq i (q:math:round (/ (* N 23.1) (/ (log (* k k)) (log 10)))) nn i)
    
71.   (while (>= i 0)
    
72.     (setq n2 (1+ (q:math:round (/ (* (- nn i) N) nn))))
    
73.     (if (> n2 N) (setq n2 N))
    
74.     (setq c (q:math:bigreci v (1+ (* 2 i)) n2)
    
75.           b1 (q:math:bigdiv b (* k k) n2)
    
76.           b (q:pi:deal1 b b1)
    
77.           b2 (q:math:big- c b n2)
    
78.           b (q:pi:deal1 b b2)
    
79.           i (1- i))
    
80.   )
    
81.   (setq b (q:math:bigdiv b k N))  
    
82.   (if (> zf 0) (setq a (q:math:big+  a b N)) (setq a (q:math:big- a b N)) )
    
83.   
    
84. )
    
85. 
    
86. (defun c:test( / n x x1 x2)
    
87.   (setq N (getint "\n请输入所需的PI的位数（4的倍数）:") N (+ 2 (fix (/ N 4)))
    
88.         a (q:math:inin 0 N)
    
89.         b (q:math:inin 0 N)
    
90.         x1 (q:math:arctg 5 16 1 N)
    
91.         x2 (q:math:arctg 239 4 -1 N))
    
92.   (princ "\n")
    
93.   (foreach x (reverse (cdr (reverse a))) (princ x))
    
94.   (princ)
    
95. )
    
96. 
    
97. 
    
98. (princ "\n By qjchen@gmail.com, 算Pi，命令test")
    
99. (princ)
    
100. 
     

;;; By <a href="mailto:qjchen@gmail.com">qjchen@gmail.com</a> ，采用Machin算法，不过由于采用了大数计算，速度随n增大急剧下降，<br /> ;;; 应该是O(n2)的算法，以后再改进吧，大概只能算800位内时间还能接受，更大的就不要了<br /> ;;; 远比不上Highflybird兄的外星人算法<br /> ;;; 2011.11.16<br /> (defun q:math:big+(a b n / f i res x)<br />   (setq i (1- n) f 0)<br />   (repeat n<br />     (if (> (setq x (+ (nth i a) (nth i b) f)) 10000) (setq x (- x 10000) f 1) (setq f 0))<br />     (setq res (cons x res) i (1- i))<br />   )<br />   res<br /> )<br /> <br /> (defun q:math:big-(a b n / f i res x)<br />   (setq i (1- n) f 0)<br />   (repeat n<br />     (if (< (setq x (- (nth i a) (nth i b) f)) 0) (setq x (+ x 10000) f 1) (setq f 0))<br />     (setq res (cons x res) i (1- i))<br />   )<br />   res<br /> )<br /> <br /> (defun q:math:round(x)<br /> (setq y (fix x))<br /> (if (> (- x y) 0.5) (1+ y) y)<br /> )<br /> <br /> (defun q:math:bigdiv(a b n / c f i res x)<br />   (setq i 0 f 0)<br />   (repeat n<br />     (setq c (+ (nth i a) (* f 10000))<br />           x (fix (/ (+ c 0.1) b))<br />           res (cons x res) <br />           i (1+ i) <br />           f (rem c b))<br />   )<br />   (reverse res)<br /> )<br /> <br /> (defun q:math:bigreci(f b n / c i res x)<br />   (setq res (cons (fix (/ f b)) res)<br />         f (rem f b)<br />         i 1)<br />   (while (< i n)<br />     (setq c (* f 10000)<br />           x (fix (/ (+ c 0.1) b))<br />           res (cons x res) i (1+ i) f (rem c b))<br />   )<br />   (reverse res)<br /> )<br /> <br /> (defun q:math:inin (v n / res)<br />   (repeat (1- n)<br />     (setq res (cons 0 res))<br />   )<br />   (cons v res)<br /> )<br /> <br /> (defun q:pi:deal1(a b / i n1 n2 res)<br />   (setq n1 (length a) n2 (length b) res b i 0)<br />   (if (< n2 n1)<br />      (repeat (- n1 n2)<br />        (setq res (append res (list (nth (+ n2 i) a))) i (1+ i))<br />      )<br />   )<br />   res<br /> )<br /> <br /> (defun q:math:arctg(k v zf N / i n2)<br />   (setq i (q:math:round (/ (* N 23.1) (/ (log (* k k)) (log 10)))) nn i)<br />   (while (>= i 0)<br />     (setq n2 (1+ (q:math:round (/ (* (- nn i) N) nn))))<br />     (if (> n2 N) (setq n2 N))<br />     (setq c (q:math:bigreci v (1+ (* 2 i)) n2)<br />           b1 (q:math:bigdiv b (* k k) n2)<br />           b (q:pi:deal1 b b1)<br />           b2 (q:math:big- c b n2)<br />           b (q:pi:deal1 b b2)<br />           i (1- i))<br />   )<br />   (setq b (q:math:bigdiv b k N))  <br />   (if (> zf 0) (setq a (q:math:big+  a b N)) (setq a (q:math:big- a b N)) )<br />   <br /> )<br /> <br /> (defun c:test( / n x x1 x2)<br />   (setq N (getint "\n请输入所需的PI的位数（4的倍数）:") N (+ 2 (fix (/ N 4)))<br />         a (q:math:inin 0 N) <br />         b (q:math:inin 0 N)<br />         x1 (q:math:arctg 5 16 1 N)<br />         x2 (q:math:arctg 239 4 -1 N))<br />   (princ "\n")<br />   (foreach x (reverse (cdr (reverse a))) (princ x))<br />   (princ)<br /> )<br /> <br /> <br /> (princ "\n By <a href="mailto:qjchen@gmail.com">qjchen@gmail.com</a>, 算Pi，命令test")<br /> (princ)<br /> <br />

cabinsummer

qjchen 发表于 2011-11-16 23:38
问题很有趣，前些日子编了一下，始终比Highflybird兄的外星人代码算法慢许多，采用的是Machin算法，大数计算 ...

你的方法和我的方法大同小异，如果没有对稀疏数据的优化，速度是不可能提很快的。同样的算法，在VB上运行比LISP上快10倍以上。

cabinsummer

qjchen 发表于 2011-11-16 23:38
问题很有趣，前些日子编了一下，始终比Highflybird兄的外星人代码算法慢许多，采用的是Machin算法，大数计算 ...

qjchen兄可以去挑战一下第二季，光路计算
http://bbs.mjtd.com/thread-90226-1-1.html

zqb05

学习了
支持