采用连分式求解pell方程

aimisiyou

;;; 求解形如X^2-D*Y^2=1的pell方程的最小正整数解(D为非平方数)当D=9781时X为156位数 Y为154位数

(defun sr(a)
   (setq n (strlen a) i 1 lst nil)
   (while (<= i n)
        (setq lst (cons (atoi (substr a i 1)) lst))
        (setq i (+ 1 i))
   )
   (reverse lst)
)
(defun add (lst1 lst2)
  (setq k1 (length lst1))
  (setq k2 (length lst2))
  (if (> k1 k2)
          (repeat (- k1 k2)  (setq lst2 (cons 0 lst2)) )
          (repeat (- k2 k1)  (setq lst1 (cons 0 lst1)) )
   )
  (setq addlst (mapcar '(lambda (x y) (+ x y)) lst1 lst2 ))
  addlst
)
(defun lnum (lst n n0)
   (setq lst (mapcar '(lambda (x) (* x n))  lst))
   (setq zlst (reverse lst))
   (repeat n0 (setq zlst (cons 0 zlst) ) )
(reverse zlst)
)  
(defun xc(lst1 lst2)
   (setq n1 (length lst1))
   (setq n2 (length lst2))
   (setq j 1 fflst '(0))
   (if (> n1 n2)
          (while (<= j n2)
              (setq fflst (add (lnum lst1 (nth (- j 1) lst2) (- n2 j) ) fflst))
              (setq j (+ 1 j))
           )
          (while (<= j n1)
              (setq fflst (add (lnum lst2 (nth (- j 1) lst1) (- n1 j) ) fflst))
              (setq j (+ 1 j))
           )
    )
fflst
)
(defun jinw (blst)
    (while (not (apply 'and (mapcar '(lambda (x) (<= x 9)) blst )))
         (setq hlst (mapcar '(lambda (x) (rem x 10)) blst))
         (setq qlst (mapcar '(lambda (x) (/ x 10)) blst))
         (setq qlst (append qlst (list 0)))
         (setq blst (add hlst qlst))
    )
  blst
)
(defun cheng (str1 str2)
    (setq lst1 (sr str1))
    (setq lst2 (sr str2))
    (setq tenlst (jinw (xc lst1 lst2)))
    (setq str (apply 'strcat (mapcar '(lambda (x) (itoa x)) tenlst ) ) )
    (if (= (car tenlst) 0)
      (setq str (vl-string-left-trim "0" str))
    )
str
)
(defun sadd (str1 str2)
  (setq lst1 (sr str1))
  (setq lst2 (sr str2))
  (setq addlst (add lst1 lst2))
  (setq tenlst (jinw addlst))
  (setq str (apply 'strcat (mapcar '(lambda (x) (itoa x)) tenlst ) ) )
  (if (= (car tenlst) 0)
      (setq str (vl-string-left-trim "0" str))
   )
str
)

(defun lf (n)
    (setq flst nil)
    (defun hj (lst)
       (if (/= (setq k (gcd (gcd (car lst) (cadr lst)) (gcd (cadr lst) (caddr lst)))) 1)
          (setq lst (list (/ (car lst) k) (/ (cadr lst) k) (/ (caddr lst) k)))
      )
      lst
    )
   (defun qz (lst)
        (hj lst)
        (setq va (fix (/ (+ (car lst) (* (cadr lst) (sqrt n))) (caddr lst))) )
         va
    )
   (defun hx (lst)
         (setq j (qz lst))
         (setq flst (cons j flst))
         (setq lst (list
                (- (* j (caddr lst) (caddr lst)) (* (car lst) (caddr lst)) )
                (* (cadr lst) (caddr lst) )
                (- (* n (cadr lst) (cadr lst)) (* (- (car lst) (* j (caddr lst))) (- (car lst) (* j (caddr lst)))) )
               )
          )
         (setq lst (hj lst))
          lst
     )
   (setq i (fix (sqrt n)))
   (setq flst (cons i flst))
   (setq lst (list i 1 (- n (* i i))))
   (while (not (equal lst (list i 1 1)))
       (setq lst (hx lst))
   )
  (reverse flst)
)
(defun  ff (lst)
  (if (= (rem (length lst) 2) 1)  
       (setq lst (append lst (list (* 2 (car lst))) (cdr lst)))
  )
  lst     
)

(defun pell (lst)
   (setq pelst (list (car lst) "1"))
   (foreach e (cdr lst)
          (setq pelst (list   (sadd (cadr pelst) (cheng (car pelst) e) )   (car pelst)    )    )

  )
pelst

)
(defun c:pfun ()
   (setq n (getint "请输入非完全平方数D"))
   (if (= (fix (sqrt n)) (sqrt n))
     (alert "你输入的是完全平方数！请重新输入！")
     (progn
        (setq lst (mapcar 'itoa (reverse (ff (lf n)))))
        (pell lst)
      )
   )
)

newer

给力，可以贴写背景材料吗？什么应用要解这样的方程？

/db_自贡黄明儒_

C语言中有许多经典程序，转成lisp就靠你们通多种语言的了

aimisiyou

newer 发表于 2014-11-23 20:52
给力，可以贴写背景材料吗？什么应用要解这样的方程？

背景材料我没想过，只是觉得此连分式算法比单纯while、for循环要高效多，感兴趣而已就写成lisp代码。

aimisiyou

/db_自贡黄明儒_ 发表于 2014-11-23 20:55
C语言中有许多经典程序，转成lisp就靠你们通多种语言的了

语言是实现手段而已，算法才是精髓！

newer

aimisiyou 发表于 2014-11-23 21:00
背景材料我没想过，只是觉得此连分式算法比单纯while、for循环要高效多，感兴趣而已就写成lisp代码。

什么情况下要解这样的方程呢？

aimisiyou

newer 发表于 2014-11-23 21:10
什么情况下要解这样的方程呢？

扯远点的话，估计就是数论方面的事情了。