home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numth.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  11.6 KB  |  372 lines
  1. ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
  2. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  3. ;;;     The data in this file contains enhancments.                    ;;;;;
  4. ;;;                                                                    ;;;;;
  5. ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
  6. ;;;     All rights reserved                                            ;;;;;
  7. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  8. ;;;     (c) Copyright 1981 Massachusetts Institute of Technology         ;;;
  9. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10.  
  11. (in-package "MAXIMA")
  12. ;;;   *****************************************************************
  13. ;;;   ***** NUMTH ******* VARIOUS NUMBER THEORY FUNCTIONS *************
  14. ;;;   *****************************************************************
  15.  
  16. (macsyma-module numth)
  17.  
  18. (declare-top(special primechannel $intfaclim))
  19.  
  20. (load-macsyma-macros rzmac)
  21.  
  22. (comment PRIME number generator)
  23.  
  24. (defmvar $maxprime 489318.)
  25.  
  26. (or (boundp 'primechannel) (setq primechannel nil))
  27.  
  28. ;#+ITS
  29. ;(defun open-primechannel nil
  30. ;      (setq primechannel (open '|mc:maxdmp;ptable >| '(in fixnum))
  31. ;        $maxprime (f1- (car (syscall 1 'fillen primechannel)))))
  32.  
  33. ;#+LISPM
  34. ;(defun open-primechannel nil
  35. ;       (setq primechannel
  36. ;         (open "mc:maxdmp;ptable >" '(:read :fixnum :byte-size 9))))
  37.                         
  38. #-(or cl NIL)
  39. (defun prime (n)
  40.        (cond ((or (< n 1) (> n $maxprime))
  41.           nil)
  42.          ((input-word n))))
  43.  
  44.  
  45. #+cl
  46. (defmacro byte-in (file) `(read-byte ,file))
  47. #+obsolete
  48. (defun input-word (n)
  49.        (funcall primechannel ':set-pointer (f* 4 (f1- n)))
  50.        (dpb (byte-in primechannel) 3311
  51.         (dpb (byte-in primechannel) 2211
  52.          (dpb (byte-in primechannel) 1111 (byte-in primechannel)))))
  53.  
  54. (defmfun $prime (n)
  55.   #+(or cl NIL) (MERROR "PRIME doesn't work yet in NIL.")
  56.   #-(or cl NIL) (prog2 (open-primechannel)
  57.            (if (eq (ml-typep n) 'fixnum)
  58.            (or (prime n) (list '($prime) n))
  59.            (list '($prime) n))
  60.            (close primechannel)))
  61.  
  62. (comment Sum of divisors and Totient functions)
  63.  
  64. (DEFMFUN $divsum n
  65.        (or (< n 3)
  66.        (merror "To many arguments to DIVSUM"))
  67.        ((lambda ($intfaclim k n) 
  68.         (cond ((and (integerp k) (integerp n))
  69.                (setq n (abs n))
  70.                (cond ((equal k 0)
  71.                   (cond ((= n 1) 1)
  72.                     ((= n 0) 1)
  73.                     (t (do ((l (cfactorw n) (cddr l))
  74.                         (a 1 (times a (f1+ (cadr l)))))
  75.                        ((null l) a)))))
  76.                  (t (divsum (cfactorw n) k))))
  77.               ((list '($divsum) n k))))
  78.     nil
  79.     (cond ((= n 1) 1)
  80.           ((arg 2)))
  81.     (arg 1)))
  82.  
  83. (defun divsum (l k)
  84.        (do ((l l (cddr l))
  85.         (ans 1) (temp))
  86.        ((null l) ans)
  87.        (cond ((equal (car l) 1))
  88.          ((setq temp (expt (car l) k)
  89.             ans (times
  90.                  (*quo (sub1 (expt temp (add1 (cadr l))))
  91.                    (sub1 temp))
  92.                  ans))))))
  93.  
  94. (defmfun $totient (n)
  95.   (cond ((integerp n)
  96.      (setq n (abs n))
  97.      (cond ((lessp n 1) 0)
  98.            ((equal n 1) 1)
  99.            (t (do ((factors (let ($intfaclim) (cfactorw n))
  100.                 (cddr factors))
  101.                (total 1 (times total
  102.                        (sub1 (car factors))
  103.                        (expt (car factors)
  104.                          (sub1 (cadr factors))))))
  105.               ((null factors) total)))))
  106.     (t (list '($totient) n))))
  107.  
  108. (comment JACOBI symbol and Gaussian factoring)
  109.  
  110. (declare-top(special *incl modulus $intfaclim))
  111.  
  112. (setq *incl (list 2 4))
  113.  
  114. (and (nconc *incl *incl) 'noprint)
  115.  
  116. (defun rtzerl2 (n)
  117.        (cond ((zerop n) 0)
  118.          (t (do ((n n (quotient n 4)))
  119.             ((not (zerop (haipart n -2))) n)))))
  120.  
  121. (defun imodp (p)
  122.        (cond ((not (equal (remainder p 4) 1)) nil)
  123.          ((equal (remainder p 8) 5) (imodp1 2 p))
  124.          ((equal (remainder p 24) 17) (imodp1 3 p))    ;p=2(mod 3)
  125.          (t (do ((i 5 (plus i (car j)))    ;p=1(mod 24)
  126.              (j *incl (cdr j)))
  127.             ((equal (jacobi i p) -1) (imodp1 i p))))))
  128.  
  129. (defun imodp1 (i modulus)
  130.        (abs (cexpt i (quotient (sub1 modulus) 4) )))
  131.  
  132. (DEFMFUN $jacobi (p q)
  133.        (cond ((null (and (integerp p) (integerp q)))
  134.           (list '($jacobi) p q))
  135.          ((zerop q) (merror "Zero denominator?"))
  136.          ((minusp q) ($jacobi p (minus q)))
  137.          ((and (evenp (setq q (rtzerl2 q)))
  138.            (setq q (quotient q 2))
  139.            (evenp p)) 0)
  140.          ((equal q 1) 1)
  141.          ((minusp (setq p (remainder p q)))
  142.           (jacobi (rtzerl2 (plus p q)) q))
  143.          (t (jacobi (rtzerl2 p) q))))
  144.  
  145. (defun jacobi (p q)
  146.        (do ((r1 p (rtzerl2 (remainder r2 r1)))
  147.         (r2 q r1)
  148.         (bit2 (haipart q -2))
  149.         (odd 0 (boole boole-xor odd (boole  boole-and bit2
  150.                        (setq bit2 (haipart r1 -2))))))
  151.        ((zerop r1) 0)
  152.        (cond ((evenp r1) (setq r1 (quotient r1 2))
  153.                  (setq odd (boole boole-xor odd
  154.                           (lsh (^ (haipart r2 -4) 2) -2)))))
  155.        (and (equal r1 1) (return (expt -1 (boole  boole-and 1 (lsh odd -1)))))))
  156.  
  157. (defun psumsq (p)
  158.        ((lambda (x)
  159.         (cond ((equal p 2) (list 1 1))
  160.               ((null x) nil)
  161.               (t (psumsq1 p x))))
  162.     (imodp p)))
  163.  
  164. (defun psumsq1 (p x)
  165.        (do ((sp ($isqrt p))
  166.         (r1 p r2)
  167.         (r2 x (remainder r1 r2)))
  168.        ((not (greaterp r1 sp)) (list r1 r2))))
  169.  
  170.  
  171. (defun gctimes (a b c d)
  172.        (list (difference (times a c)
  173.              (times b d))
  174.          (plus (times a d)
  175.            (times b c))))
  176.  
  177.  
  178.  
  179. (declare-top(*lexpr $rat))
  180.  
  181.  
  182. ;(DEFMFUN $gcfactor (n)
  183. ;       (setq n (cdr ($totaldisrep ($bothcoef ($rat n '$%i) '$%i))))
  184. ;       (cond ((and (integerp (car n)) (integerp (cadr n)))
  185. ;          
  186. ;          (setq n (map2c #'(lambda (term exp)
  187. ;                 (cond ((= exp 1) (gcdisp term))
  188. ;                       (t (list '(mexpt) (gcdisp term) exp))))
  189. ;                 (gcfactor (cadr n) (car n))))
  190. ;          (cond ((null (cdr n)) (car n))
  191. ;            (t (cons '(mtimes simp) (nreverse n)))))
  192. ;         (t (gcdisp (nreverse n)))))
  193.  
  194.  
  195. (DEFMFUN $gcfactor (n)
  196.        (setq n (cdr ($totaldisrep ($bothcoef ($rat n '$%i) '$%i))))
  197.        (cond ((and (integerp (car n)) (integerp (cadr n)))
  198.           (setq n (sloop for (term exp) on (gcfactor (cadr n) (car n))
  199.                 collecting
  200.                  (cond ((= exp 1) (gcdisp term))
  201.                        (t (list '(mexpt) (gcdisp term) exp)))))
  202.           (cond ((null (cdr n)) (car n))
  203.             (t (cons '(mtimes simp) (nreverse n)))))
  204.          (t (gcdisp (nreverse n)))))
  205.  
  206. (defun gcdisp (term)
  207.        (cond ((atom term) term)
  208.          ((let ((rp (car term))
  209.             (ip (cadr term)))
  210.            (setq ip (cond ((equal ip 1) '$%i)
  211.                   (t (list '(mtimes) ip '$%i))))
  212.            (cond ((equal rp 0) ip)
  213.              (t (list '(mplus) rp ip)))))))
  214. ;
  215. ;(defun gcfactor (a b &aux tem) 
  216. ;       (prog (gl cd dc econt p e1 e2 ans plis nl $intfaclim)
  217. ;       (setq e1 0
  218. ;         e2 0
  219. ;         econt 0
  220. ;         gl (gcd a b)
  221. ;         a (quotient a gl)
  222. ;         b (quotient b gl)
  223. ;         nl (cfactorw (plus (times a a) (times b b)))
  224. ;         gl (cfactorw gl))
  225. ;       (and (equal 1 (car gl)) (setq gl nil))
  226. ;       (and (equal 1 (car nl)) (setq nl nil))
  227. ;loop   (show e1 e2 ans gl nl)
  228. ;       (cond ((null gl)
  229. ;          (cond ((null nl) (go ret))
  230. ;            ((setq p (car nl)))))
  231. ;         ((null nl) (setq p (car gl)))
  232. ;         (t (setq p (max (car gl) (car nl)))))
  233. ;       (setq cd (psumsq p))
  234. ;       (show (list p cd ))
  235. ;       (cond ((null cd)
  236. ;          (setq plis (cons p (cons (cadr gl) plis)))
  237. ;          (setq gl (cddr gl)) (go loop))
  238. ;         ((equal p (car nl))
  239. ;          (cond ((zerop (remainder (setq tem (plus (times a (car cd))    ;gcremainder
  240. ;                         (times b (cadr cd))))
  241. ;                       p))    ;remainder(real((a+bi)cd~),p)   z~ is complex conjugate
  242. ;             (setq e1 (cadr nl)) (setq dc cd))
  243. ;            (t (setq e2 (cadr nl))
  244. ;               (setq dc (reverse cd))))
  245. ;          (show tem dc)
  246. ;          (setq dc (gcexpt dc (cadr nl))    ;
  247. ;            dc (gctimes a b (car dc) (minus (cadr dc)))
  248. ;            a (quotient (car dc) p)
  249. ;            b (quotient (cadr dc) p)
  250. ;            nl (cddr nl))))
  251. ;       (cond ((equal p (car gl))
  252. ;          (setq econt (plus econt (cadr gl)))
  253. ;          (cond ((equal p 2)
  254. ;             (setq e1 (f+ e1 (f* 2 (cadr gl)))))
  255. ;            (t (setq e1 (f+ e1 (cadr gl))
  256. ;                 e2 (f+ e2 (cadr gl)))))
  257. ;          (setq gl (cddr gl))))
  258. ;       (show a b e1 e2 dc cd)
  259. ;       (and (not (zerop e1))
  260. ;        (setq ans (cons cd (cons e1 ans)))
  261. ;        (setq e1 0))
  262. ;       (and (not (zerop e2))
  263. ;        (setq ans (cons (reverse cd) (cons e2 ans)))
  264. ;        (setq e2 0)) 
  265. ;       (go loop) 
  266. ;ret    (show 'ret ans)
  267. ;       (setq cd (gcexpt (list 0 -1)
  268. ;            (remainder econt 4)))
  269. ;       (setq a (gctimes a b (car cd) (cadr cd)))
  270. ;       (cond ((or (equal (car a) -1) (equal (cadr a) -1))
  271. ;          (setq plis (cons -1 (cons 1 plis)))))
  272. ;       (cond ((equal (car a) 0)
  273. ;          (setq ans (cons '(0 1) (cons 1 ans)))))
  274. ;       (return (nconc plis ans))))
  275.  
  276. ;(defun test-gcfactor (n &aux facts numb prod orig-facts prod1 tem dif)
  277. ;  (setq orig-facts  (sloop for i below (f+ 2 (random n))
  278. ;             do (setq tem (list (random 10) (random 12)))
  279. ;             when (not (equal tem (list 0 0)))
  280. ;             collecting tem
  281. ;             collecting (f1+ (random 5))))
  282. ;;  (show orig-facts)
  283. ;  (setq prod (multiply-gcfactors orig-facts))
  284. ;;  (show prod)
  285. ;  (setq numb (add* (car prod) (mul* '$%i (second prod))))
  286. ;  (displa numb)
  287. ;  (setq facts  ($gcfactor numb))
  288. ;  (displa facts)
  289. ;  (setq prod1 ($ratsimp facts))
  290. ;  (assert (equal 0 (setq dif ($ratsimp (sub* numb prod1)))))
  291. ;;  (show dif)
  292. ;  dif)
  293.     
  294. (defun gcfactor (a b &aux tem) 
  295.        (prog (gl cd dc econt p e1 e2 ans plis nl $intfaclim )
  296.        (setq e1 0
  297.          e2 0
  298.          econt 0
  299.          gl (gcd a b)
  300.          a (quotient a gl)
  301.          b (quotient b gl)
  302.          nl (cfactorw (plus (times a a) (times b b)))
  303.          gl (cfactorw gl))
  304.        (and (equal 1 (car gl)) (setq gl nil))
  305.        (and (equal 1 (car nl)) (setq nl nil))
  306. loop
  307.        (cond ((null gl)
  308.           (cond ((null nl) (go ret))
  309.             ((setq p (car nl)))))
  310.          ((null nl) (setq p (car gl)))
  311.          (t (setq p (max (car gl) (car nl)))))
  312.        (setq cd (psumsq p))
  313.        (cond ((null cd)
  314.           (setq plis (cons p (cons (cadr gl) plis)))
  315.           (setq gl (cddr gl)) (go loop))
  316.          ((equal p (car nl))
  317.           (cond ((zerop (remainder (setq tem (plus (times a (car cd))    ;gcremainder
  318.                          (times b (cadr cd))))
  319.                        p))    ;remainder(real((a+bi)cd~),p)   z~ is complex conjugate
  320.              (setq e1 (cadr nl)) (setq dc cd))
  321.             (t (setq e2 (cadr nl))
  322.                (setq dc (reverse cd))))
  323.           (setq dc (gcexpt dc (cadr nl))    ;
  324.             dc (gctimes a b (car dc) (minus (cadr dc)))
  325.             a (quotient (car dc) p)
  326.             b (quotient (cadr dc) p)
  327.             nl (cddr nl))))
  328.        (cond ((equal p (car gl))
  329.           (setq econt (plus econt (cadr gl)))
  330.           (cond ((equal p 2)
  331.              (setq e1 (f+ e1 (f* 2 (cadr gl)))))
  332.             (t (setq e1 (f+ e1 (cadr gl))
  333.                  e2 (f+ e2 (cadr gl)))))
  334.           (setq gl (cddr gl))))
  335.        (and (not (zerop e1))
  336.         (setq ans (cons cd (cons e1 ans)))
  337.         (setq e1 0))
  338.        (and (not (zerop e2))
  339.         (setq ans (cons (reverse cd) (cons e2 ans)))
  340.         (setq e2 0))
  341.        (go loop) 
  342. ret    (setq cd (gcexpt (list 0 -1)
  343.             (remainder econt 4)))
  344.        (setq a (gctimes a b (car cd) (cadr cd)))
  345.        ;;a hasn't been divided by p yet..
  346.        (setq a (mapcar 'signum a)) 
  347.        #+cl (assert (or (zerop (car a))(zerop (second a))))
  348.        (cond ((or (equal (car a) -1) (equal (cadr a) -1))
  349.           (setq plis (cons -1 (cons 1 plis)))))
  350.        (cond ((equal (car a) 0)
  351.           (setq ans (cons '(0 1) (cons 1 ans)))))
  352.        (setq ans (nconc plis ans))
  353.        (return ans)))
  354.  
  355. (defun multiply-gcfactors (lis)
  356.   (sloop for (term exp) on (cddr lis) by 'cddr
  357.     with answ = (cond ((numberp (car lis))(list (pexpt (car lis) (second lis)) 0))
  358.               (t(gcexpt (car lis) (second lis))))
  359.     when (numberp term)
  360.     do (setq answ (list (times (first answ) term) (times (second answ) term)))
  361.     (show answ)
  362.     else
  363.     do (setq answ (apply 'gctimes (append answ (gcexpt term exp))))
  364.     finally (return answ)))
  365.  
  366. (defun gcexpt (a n)
  367.        (cond ((zerop n) '(1 0))
  368.          ((equal n 1) a)
  369.          (t (gctime1 a (gcexpt a (f1- n))))))
  370.  
  371. (defun gctime1 (a b) (gctimes (car a) (cadr a) (car b) (cadr b)))
  372.