; .EnTete "Le-Lisp (c) version 15.2" " " "Les Complexes"
; .EnPied " " "%" " "
; .Chapitre IV "Les Nombres Complexes"
;
; .Centre "*****************************************************************"
; .Centre " Ce fichier est en lecture seule hors du projet ALE de l'INRIA. "
; .Centre " Il est maintenu par ILOG SA, 2 Avenue Gallie'ni, 94250 Gentilly "
; .Centre " (c) Le-Lisp est une marque de'pose'e de l'INRIA "
; .Centre "*****************************************************************"
; .Centre "$Header: complex.ll,v 4.3 88/02/16 10:22:39 kuczynsk Rel213 $"
;==========================================================================
;
; Le←Lisp v15.2 : les complexes C.
;
;===========================================================================
; (c) Jean Vuillemin, 1986
; Institut National de Recherche en Informatique et Automatique
; B.P. 105, 78153 Le Chesnay Cedex, France. vuillemin@inria
;===========================================================================
(unless (>= (version) 15.2)
(error 'load 'erricf 'ratio))
(setq #:sys-package:colon 'R)
(add-feature 'complex)
(unless (featurep 'genr) (loadmodule 'genr))
(setq #:backup:majuscules #:system:read-case-flag
#:system:read-case-flag ())
; Un complexe est repre'sente' c par sa partie re'elle (:r c) et sa partie
; imaginaire (:i c).
(defstruct :C r i)
; Teste le type C:
(dmd :C:? (c) `(eq ':C (type-of ,c)))
; Construit un complexe [r:i]. Suppose i<>0 :
(dmd :C:x (r i) `(let ((:c (:C:make))) (:C:r :c ,r) (:C:i :c ,i) :c))
(setq #:sys-package:colon '#:R:C)
; Construit a+ib = [a:b]
(de makecomplex (a b) (if (= 0 b) a (:x a b)))
; L'imaginaire pur i*i=-1 :
(defvar i (makecomplex 0 1))
; Teste si un nombre est un complexe pur:
(de complexp (c) (and (:? c) c))
; Teste si un objet Lisp est un re'el.
(de realp (r) ; Les re'els sont forme's de
(or (rationalp r) ; Q,
(floatp r) ; des flottants,
(and (numberp r) ; et de tous les nombres qui
(send 'realp r)))) ; re'pondent au message realp.
; Les complexes ne sont pas des re'els:
(de :realp (c) ())
; .Section "Fonctions propres aux complexes"
; Acce's aux parties re'elles et imaginaires :
(de realpart (c)
(cond ((:? c) (:r c))
((numberp c) c)
(T (#:R:error 'realpart 'ERRNNA (list 'realpart c)))))
(de imagpart (c)
(cond ((:? c) (:i c))
((numberp c) 0)
(T (#:R:error 'imagpart 'ERRNNA (list 'imagpart c)))))
; Le conjugue' de a+ib est a-ib
(de conjugate (c)
(if (numberp c)
(if (:? c) (:x (:r c) (- (:i c))) c)
(#:R:error 'conjugate 'ERRNNA
(list 'conjugate c))))
; Definition du module au carre:
(de :module2 (c) (+ (* (:r c) (:r c)) (* (:i c) (:i c))))
; .Section "Arithme'tique Complexe"
(de :float (c) (:x (float (:r c)) (float (:i c))))
(de :truncate (c) (#:R:error 'truncate 'RINC (list 'fix c)))
; La comparaison est rarement bien de'finie sur les complexes:
; On renvoie le complexe forme' du re'sultat
; des comparaisons sur les 2 champs.
(de :<?> (a b)
(if (:? b)
(makecomplex (<?> (:r a) (:r b)) (<?> (:i a) (:i b)))
(:x (<?> (:r a) b) (<?> (:i a) 0))))
; Addition binaire :
(de :+ (c r)
(if (:? r)
(makecomplex (+ (:r c) (:r r)) (+ (:i c) (:i r)))
(:x (+ r (:r c)) (:i c))))
; Negation Unaire:
(de :0- (c) (:x (- (:r c)) (- (:i c))))
; Multiplication binaire :
(de :* (c r)
(if (:? r)
(makecomplex (- (* (:r c) (:r r)) (* (:i c) (:i r)))
(+ (* (:r c) (:i r)) (* (:i c) (:r r))))
(makecomplex (* r (:r c)) (* r (:i c)))))
; Inverse: 1/a+ib = a/m-ib/m avec m=a*a+b*b
(de :1/ (c)
(let ((m (:module2 c)))
(:x (/ (:r c) m) (- (/ (:i c) m)))))
; Par convention, a+ib s'e'crit [bi+a]
(de :prin (c)
(cond ((eq (precision) 'cl)
(princn #/#) (princn #/C) (princn #/()
(prin (:r c)) (princn #\SP) (prin (:i c)) (princn #/)))
(T (princn #/[)
(unless (or (= 1 (:i c)) (= 1. (:i c)))
(if (or (= -1 (:i i)) (= -1. (:i c)))
(princn #/-)
(prin (:i c))))
(princn #/i)
(unless (= 0 (:r c))
(cond ((> (:r c) 0)
(princn #/+) (prin (:r c)))
(T (princn #/-) (prin (- (:r c))))))
(princn #/])))
c)
; Pour pouvoir relire les complexes:
(dmc |[| ()
(let ((r) (j))
(setq j (with ((typecn #/i 'csep))
(selectq (peekcn)
(#/i 1)
(#/- (readcn)
(if (eq (peekcn) #/i) -1 (- (read))))
(T (read))))
r (readcn))
(selectq r
(#/i (selectq (readcn)
(#/] (setq r 0))
(#/+ (setq r (read))
(if (eq (peekcn) #/]) (readcn)
(#:R:error 'read 'RNOTC r)))
(#/- (setq r (- (read)))
(if (eq (peekcn) #/]) (readcn)
(#:R:error 'read 'RNOTC r)))
(T (#:R:error 'read 'RNOTC r))))
(#/] (setq r j j 0))
(T (#:R:error 'read 'RNOTC r)))
(makecomplex r j)))
(defsharp |C| ()
(setq n (read))
(list (:x (car n) (cadr n))))
; .Section "Conversion en coordonne'es polaires"
; De'finition de pi : soit on prend la valeur exacte si on dispose des re'els,
; soit on approxime en flottant avec une pre'cision de'pendante de la machine:
(defvar pi (if (boundp 'pi) pi (* 4 (atan 1.))))
(defvar pi/2 (/ pi 2))
; Pour un complexe c=rho*e**i*theta, le module rho est donne' par (:abs c),
; la phase theta par (:phase c) et e**i*theta par (:signum c).
(de :abs (c) (if (zerop (:r c)) (abs (:i c)) (sqrt (:module2 c))))
; La phase est 0 pour un re'el, et atan(r/j) pour c=[ji+r]
(de phase (c)
(cond ((:? c) (:arctan (:i c) (:r c)))
((numberp c) (if (< c 0) pi 0))
(T (#:R:error 'phase 'ERRNNA
(list 'phase c)))))
; i*r = r*e**i*pi/2 rc+ic=rho*e**i*theta, tan(theta)=ic/rc
(de :arctan (y x)
(selectq (<?> y 0)
(1 (selectq (<?> x 0)
(1 (atan (/ y x)))
(0 pi/2)
(-1 (+ pi (atan (/ y x))))))
(0 (if (< x 0) pi 0))
(-1 (- (:arctan (- y) x)))))
; Pour un re'el, (signum r) donne le signe -1,0,1 de r.
; Pour un complexe c=rho*e**i*theta, c'est e**i*theta, soit c/(abs c).
(de signum (c)
(cond ((:? c) (:* c (/ (:abs c))))
((numberp c) (<?> c 0))
(T (#:R:error 'signum 'ERRNNA
(list 'signum c)))))
; .Section "Logarithme et Exponentielle sur C"
; log(r*e**i*theta) = log(r)+i*theta
(de :log (c) (:x (log (:abs c)) (:arctan (:i c) (:r c))))
; e**rc+ic = e**rc*(cosic+i*sinic)
(de :exp (c) (:* (:cis (:i c)) (exp (:r c))))
(de :cis (theta) (:x (cos theta) (sin theta)))
; (cis r) = cos(r)+i*sin(r)
(de cis (r)
(if (numberp r) (:cis r)
(#:R:error 'cis 'ERRNNA (list 'cis r))))
; (sqrt c) = c**1/2 = e**(log(c)/2)
(de :sqrt (c) (:exp (/ (:log c) 2)))
; .Section "Sinus, Cosinus et Arc Tangente sur C"
; sin(c)=1/2i(e**i*c-e**-i*c)
(de :sin (c)
(* (- (exp (* i c)) (exp (* [-i] c))) [-1/2i]))
; cos(c)=1/2(e**i*c+e**-i*c)
(de :cos (c)
(/ (+ (exp (* i c)) (exp (* [-i] c))) 2))
; atan(z)=-ilog((1+iz)(sqrt1/(1+z*z)))
(de :atan (z)
(if (= i z) (/ 0)
(:* [-i] (log (/ (+ 1 (* i z)) (sqrt (+ 1 (* z z))))))))
; tan(z) = sin(z)/co
(de tan (z) (/ (sin z) (cos z)))
; .Section "Extention sur C des Fonctions Transcendantes sur R"
(setq #:sys-package:colon '#:R)
(de :log (r)
(selectq (<?> r 0)
(1 (log (float r)))
(0 (- (/ 0)))
(-1 (:C:x (log (- r)) pi))))
(de #:fix:log (r) (:log r))
(de #:float:log (r) (:log r))
(de :sqrt (r)
(selectq (<?> r 0)
(1 (sqrt (float r)))
(0 0)
(-1 (:C:x 0 (sqrt (- r))))))
(de #:fix:sqrt (r) (:sqrt r))
(de #:float:sqrt (r) (:sqrt r))
; asin(z)=-ilog(iz+sqrt(1-z*z))
(de :asin (z)
(* [-i] (log (+ (* z i) (sqrt (- 1 (* z z)))))))
(de #:fix:asin (r) (:asin r))
(de #:float:asin (r) (:asin r))
; acos(z)=pi/2-asinz
(de :acos (z) (- pi/2 (asin z)))
(de #:fix:acos (r) (:acos r))
(de #:float:acos (r) (:acos r))
; .SSection "Fonctions Hyperboliques"
; sinhx=(e**x-e**-x)/2
(de sinh (x)
(if (numberp x)
(/ (- (exp x) (exp (- x))) 2)
(#:R:error ':#:R:error 'ERRNNA
(list 'sinh x))))
; coshx=(e**x+e**-x)/2
(de cosh (x)
(if (numberp x)
(/ (+ (exp x) (exp (- x))) 2)
(#:R:error ':#:R:error 'ERRNNA
(list 'cosh x))))
; tanhx=sinhx/coshx
(de tanh (x)
(if (numberp x)
(/ (sinh x) (cosh x))
(#:R:error ':#:R:error 'ERRNNA
(list 'tanh x))))
; asinhx=log(x+sqrt(1+x*x))
(de asinh (x)
(if (numberp x)
(log (+ x (sqrt (+ 1 (* x x)))))
(#:R:error ':#:R:error 'ERRNNA
(list 'asinh x))))
; acoshx=log(x+(x+1)*sqrt(x-1/x+1))
(de acosh (x)
(if (numberp x)
(log (+ x (* (+ x 1) (sqrt (/ (- x 1) (+ x 1))))))
(#:R:error ':#:R:error 'ERRNNA
(list 'acosh x))))
; atanh(x) = 1/2*log((1+x/1-x))
(de atanh (x)
(cond
((not (numberp x))
(#:R:error ':#:R:error 'ERRNNA (list 'atanh x)))
((= 1 x) (/ 0))
((= -1 x) (- (/ 0)))
(t (/ (log (/ (+ x 1) (- 1 x))) 2))))
; .Section "Modifs a` genr.ll"
(de :+ (n m)
(if (complexp m) (:C:+ m n) (plus (float n) (float m))))
(de :* (n m)
(if (complexp m) (:C:* m n) (times (float n) (float m))))
(de :<?> (n m)
(if (complexp m) (- (:C:<?> m n)) (<?> (float n) (float m))))
(setq #:system:read-case-flag #:backup:majuscules)