src/nalgfa.lisp

   1
   2 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   3 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   4 ;;;     The data in this file contains enhancements.                   ;;;;;
   5 ;;;                                                                    ;;;;;
   6 ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
   7 ;;;     All rights reserved                                            ;;;;;
   8 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   9 ;;;     (c) Copyright 1980 Massachusetts Institute of Technology         ;;;
  10 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  11
  12 (in-package :maxima)
  13
  14 (macsyma-module nalgfa)
  15
  16 (declare-top (special vlist *nosplitf *algvar *denom *num *ans))
  17
  18 (load-macsyma-macros rzmac ratmac)
  19
  20 (defun new-alg ()
  21   (newsym ($gensym "alg")))
  22
  23
  24 (defun psqfrp (p var)
  25   (zerop (pdegree (pgcd p (pderivative p var)) var)))
  26
  27 (defun pgsubst (val var p) ;;generalized psubst substitutes any
  28   (cond ((pcoefp p) p) ;;expression for any var in p
  29         ((eq var (car p))
  30          (cond ((pzerop val)
  31                 (ptterm (cdr p) 0))
  32                ((do ((ld (cadr p) (car a))
  33                      (a (cdddr p) (cddr a))
  34                      (ans (caddr p)
  35                           (pplus
  36                            (ptimes ans (pexpt val (- ld (car a))))
  37                            (cadr a))))
  38                     ((null a) (ptimes ans (pexpt val ld)))))))
  39         ((pointergp var (car p)) p)
  40         ((do ((a (cdddr p) (cddr a))
  41               (ans (ptimes (list (car p) (cadr p) 1)
  42                            (pgsubst val var (caddr p)))
  43                    (pplus ans
  44                           (ptimes (list (car p) (car a) 1)
  45                                   (pgsubst val var (cadr a))))))
  46              ((null a) ans)))))
  47
  48 (defun pvsubst (nvar ovar p)
  49   (cond ((or (pcoefp p) (pointergp ovar (car p))) p)
  50         ((eq ovar (car p))
  51          (cons nvar (cdr p)))
  52         (t (pgsubst (make-poly nvar) ovar p))))
  53
  54 (defun ordervar (var l)
  55   (let ((mvar (lmainvar l)))
  56     (cond ((null mvar) l)
  57           ((null (pointergp mvar var)) (cons var l))
  58           ((let ((newvar (gensym)))
  59              (setq genvar (append genvar (list newvar)))
  60              (setf (symbol-plist newvar) (symbol-plist var))
  61              (setf (symbol-value newvar) (1+ (symbol-value mvar)))
  62              (cons newvar (mapcar #'(lambda (p) (pvsubst newvar var p)) l)))))))
  63
  64 (defun lmainvar (l) ;;main var of list of poly's
  65   (do ((l l (cdr l))
  66        (v))
  67       ((null l) v)
  68     (cond ((pcoefp (car l)))
  69           ((null v) (setq v (caar l)))
  70           ((pointergp (caar l) v)
  71            (setq v (caar l))))))
  72
  73 (defun presult (p1 p2 var) ;;change call in algsys?
  74   (let ((genvar genvar))
  75     (setq var (ordervar var (list p1 p2))
  76           p1 (cadr var)
  77           p2 (caddr var)
  78           var (car var))
  79     (cond ((zerop (pdegree p1 var))
  80            (cond ((zerop (pdegree p2 var)) 1)
  81                  ((pexpt p1 (cadr p2)))))
  82           ((zerop (pdegree p2 var))
  83            (pexpt p2 (cadr p1)))
  84           ((resultant p1 p2)))))
  85
  86 (defun pcoefvec (p)
  87   (cond ((pcoefp p) (list p))
  88         ((do ((l)
  89               (i (cadr p) (1- i))
  90               (p (cdr p)))
  91              ((signp l i) (nreverse l))
  92            (push (cond ((and p (= (car p) i))
  93                         (prog1 (cadr p) (setq p (cddr p))))
  94                        ( 0 ))
  95                  l)))))
  96
  97 (defun algtrace1 (bvec tvec)
  98   (do ((i (- (length tvec) (length bvec)) (1- i)))
  99       ((zerop i) (algtrace* bvec tvec))
 100     (setq bvec (cons 0 bvec))))
 101
 102 (defun algtrace* (bvec tvec)
 103   (do ((b bvec (cdr b))
 104        (tr (car (last bvec))
 105            (pplus tr (car (last b)))))
 106       ((null (cdr b)) tr)
 107     (or (pzerop (car b))
 108         (do ((l (cdr b) (cdr l))
 109              (tv tvec (cdr tv)))
 110             ((null l))
 111           (rplaca l (pdifference (car l) (ptimes (car b) (car tv))))))))
 112
 113 (defun algtrace (f p)
 114   (let* ((r (cadr (pdivide (car f) p))))
 115     (if (or (constantp (car r))                  ; r constant
 116             (not (eq (caar r) (car p))))          ; r constant in main var of p
 117         (ratreduce (pctimes (cadr p) (car r)) (cdr r)) ; r*deg(p)
 118         (ratreduce (algtrace1 (pcoefvec (car r)) (cdr (pcoefvec p)))
 119                     (cdr r)))))
 120
 121
 122 (defmfun $algtrace (f p var)
 123   (let ((varlist (list var))
 124         (genvar nil))
 125     (rdis* (algtrace (rform f) (car (rform p))))))
 126
 127
 128 (defun good-form (l) ;;bad -> good
 129   (do ((l l (cddr l))
 130        (ans))
 131       ((null l) (nreverse ans))
 132     (push (cons (cadr l) (car l)) ans)))
 133
 134 (defun bad-form (l) ;;good -> bad
 135   (mapcar #'(lambda (q) (list (cdr q) (car q))) l))
 136
 137 (defmfun $algfac (a1 &optional (a2 nil a2?) (a3 nil a3?))
 138   (if a3?
 139       ($pfactoralg a1 a2 a3)
 140       (let ((varlist))
 141         (cond (a2?
 142                (newvar a2)
 143                (if (alike1 a2 (car varlist))
 144                    ($pfactoralg a1 nil a2)
 145                    ($pfactoralg a1 a2 (car (last varlist)))))
 146               (t
 147                (newvar a1)
 148                  (setq varlist (mapcan #'(lambda (q) (if (algpget q) (list q) nil)) varlist))
 149                  (cond ((= (length varlist) 1)
 150                         ($pfactoralg a1 nil (car varlist)))
 151                        ((> (length varlist) 1)
 152                         ;; MEANING OF NEXT MESSAGE IS UNCLEAR
 153                         (merror (intl:gettext "algfac: too many algebraics.")))
 154                        (t
 155                         ;; MEANING OF NEXT MESSAGE IS UNCLEAR
 156                         (merror (intl:gettext "algfac: no algebraics.")))))))))
 157
 158 (defmfun $pfactoralg (f p alg)
 159   (let ((varlist (list alg))
 160         (genvar) (vlist) (tellratlist) ($ratfac)
 161         ($gcd '$algebraic)
 162         ($algebraic) ($ratalgdenom t)
 163         (*denom 1) (*num 1) (*ans))
 164     (cond ((and (null p) (radfunp alg t)) (newvar (cadr alg)))
 165           (t (newvar p)))
 166     (newvar1 f)
 167     (cond ((null vlist) (merror (intl:gettext "pfactoralg: attempt to factor a constant."))))
 168     (setq varlist (nconc varlist (sortgreat vlist)))
 169     (cond (p (setq p (cadr (ratrep* p)))
 170              (push (cons alg (mapcar #'pdis (cdr p)))
 171                    tellratlist))
 172           (t (setq p (algpget alg))
 173              (setq p (pdifference
 174                       (pexpt (cadr (ratrep* alg)) (car p))
 175                       (cadr p)))))
 176     (setq $algebraic t)
 177     (setq f (cadr (ratrep* f)))
 178     (setq f (pfactoralg1 f p 0))
 179     (cons '(mtimes)
 180           (cons (rdis (ratreduce *num *denom))
 181                 (mapcar 'pdis f)))))
 182
 183 (defun nalgfac (p mp)
 184   (let ((*num 1) (*denom 1) (*ans) (algfac*) ($nalgfac)
 185         ($algebraic t) ($gcd '$algebraic))
 186     (setq p (pfactoralg1 p mp 0))
 187     (setq adn* (* adn* *denom))
 188     (cond ((equal *num 1) p)
 189           (t (cons *num p)))))
 190
 191 (setq *nosplitf t)
 192
 193 (defun pfactoralg1 (f p ind)
 194   (cond ((pcoefp f) (setq *num (ptimeschk f *num)) *ans)
 195         ((= (cadr f) 1) (setq f (pshift f (car p) ind))
 196          (push (algnormal f) *ans)
 197          (setq f (rquotient f (car *ans))
 198                *denom (ptimeschk (cdr f) *denom)
 199                *num (ptimeschk (car f) *num))
 200          *ans)
 201         ((equal (cdr f) (cdr p))
 202          (push (pdifference (make-poly (car f)) (make-poly (car p))) *ans)
 203          (setq f (rquotient f (car *ans))
 204                *denom (ptimeschk (cdr f) *denom))
 205          (pfactoralg1 (car f) p ind))
 206         ((zerop (pdegree f (car p)))
 207          (mapc #'(lambda (q)
 208                    (if (pcoefp q) nil
 209                      (pfactoralg1 (pshift q (car p) -1) p (1+ ind))))
 210                (let (($algebraic nil)
 211                      ($gcd (car *gcdl*)))
 212                  (pfactor1 f)))
 213          *ans)
 214         (t (do ((l (let (($algebraic nil)
 215                          ($gcd (car *gcdl*)))
 216                      (pfactor (algnorm f p)))
 217                    (cddr l))
 218                 (polys)
 219                 (temp)
 220                 (alg (car p)))
 221                ((null l)
 222                 (setq *num (ptimeschk f *num))
 223                 (mapc #'(lambda (q) (pfactoralg1
 224                                      (pshift q alg -1) p (1+ ind)))
 225                       polys)
 226                 *ans)
 227                (cond ((pcoefp (car l)) nil)
 228                      (t (setq temp (cond ((null (cddr l)) f)
 229                                          (t (pgcd f (car l)))))
 230                         (cond ((pcoefp temp) nil)
 231                               ((= (cadr temp) 1)
 232                                (setq temp (algnormal temp))
 233                                (push (pshift temp alg ind) *ans))
 234                               ((= (cadr l) 1)
 235                                (setq temp (algnormal temp))
 236                                (push (pshift temp alg ind) *ans)
 237                                (or *nosplitf
 238                                    (setq *nosplitf
 239                                          (list (car l) temp ind))))
 240                               (t (push temp polys)))
 241                         (setq f (rquotient f temp)
 242                               *denom (ptimeschk (cdr f) *denom)
 243                               f (car f)))) ))))
 244
 245 (defun pshift (f alg c)
 246   (if (= c 0) f
 247       (pgsubst (pplus (make-poly (car f)) (pctimes c (make-poly alg)))
 248                (car f) f)))
 249
 250
 251
 252 (defmfun $splitfield (p var)
 253   (let ((varlist)
 254         (genvar)
 255         (genpairs)
 256         (*algvar)
 257         ($gcd '$algebraic))
 258     (newsym var)
 259     (setq *algvar (caar (new-alg)))
 260     (setq p (psplit-field (cadr (ratf p))))
 261     (cons
 262      '(mlist)
 263      (cons (pdis* (car p))
 264            (mapcar 'rdis* (cdr p))))))
 265
 266 (defun psplit-field (p)                 ;modresult?
 267   (let ((l (mapcar #'(lambda (q) (psplit-field1 (cdr q)))
 268                    (good-form (pfactor p)))) ;don't normalize lcfs?
 269         ($algebraic t))
 270     (if (null (cdr l)) (car l)
 271         (do ((l l (cdr l))
 272              (prim) (zeroes) (temp))
 273             ((null l) (cons (or prim '|$splits in q|) zeroes))
 274           (cond ((eq (caar l) 'linear)
 275                  (setq zeroes (cons (cdar l) zeroes)))
 276                 ((null prim)
 277                  (setq prim (caar l)
 278                        zeroes (nconc (cdar l) zeroes)))
 279                 ((setq temp
 280                        (primelmt (cons (car p) (cdr prim))
 281                                  (cons (car p) (cdaar l))
 282                                  *algvar)
 283                        zeroes
 284                        (nconc
 285                         (mapcar
 286                          #'(lambda (q)
 287                              (ratgsubst (cadddr temp) (caaar l) q))
 288                          (cdar l))
 289                         (mapcar
 290                          #'(lambda (q)
 291                              (ratgsubst (caddr temp) (car prim) q))
 292                          zeroes))
 293                        prim (car temp))))))))
 294
 295 (defun plsolve (p)
 296   (ratreduce (ptterm (cdr p) 0)
 297              (pminus (ptterm (cdr p) 1))))
 298
 299
 300 (defun psplit-field1 (p)
 301   ;;returns minimal poly and list of roots
 302   ;;p must be square free
 303   (*bind* ((minpoly (cons *algvar (cdr p)))
 304            (zeroes) ($algebraic t)
 305            ($ratalgdenom t))
 306           (if (equal (cadr p) 1) (return (cons 'linear (plsolve p))))
 307           (do ((polys (list p) )
 308                (nminpoly)
 309                (*nosplitf nil nil)
 310                (alpha (cons (make-poly (car minpoly)) 1)))
 311               ((null polys)
 312                (cons minpoly zeroes))
 313             (push alpha zeroes)
 314             (putprop (car minpoly) (cdr minpoly) 'tellrat)
 315             (rplaca polys
 316                     (car
 317                      (rquotient (pctimes (cdr alpha) (car polys))
 318                                 (pdifference
 319                                  (pctimes (cdr alpha) (pget (caar polys)))
 320                                  (car alpha)))))
 321             (setq polys
 322                   (mapcan
 323                    #'(lambda (q)
 324                        (cond ((equal (cadr q) 1) ;;linear factor
 325                               (push (plsolve q) zeroes)
 326                               nil) ;;flush linear factor
 327                              ((list q))))
 328                    (mapcan #'(lambda (q)
 329                                (let ((*ans) (*num 1) (*denom 1))
 330                                  (nreverse (pfactoralg1 q minpoly 0))))
 331                            polys)))
 332             (when *nosplitf
 333               (setq nminpoly (car *nosplitf)
 334                     *nosplitf (cdr *nosplitf))
 335               (putprop *algvar (cdr nminpoly) 'tellrat)
 336               (let ((beta
 337                      (plsolve (pgcd (cons (caar *nosplitf) (cdr minpoly))
 338                                     (exchangevar (car *nosplitf) *algvar)))))
 339                 (setq alpha (ratplus (cons (make-poly *algvar) 1)
 340                                      (rattimes (cons (- (cadr *nosplitf)) 1)
 341                                                beta t)))
 342                 (setq zeroes
 343                       (mapcar
 344                        #'(lambda (q) (ratgsubst beta (car minpoly) q))
 345                        zeroes))
 346                 (setq polys
 347                       (mapcar
 348                        #'(lambda (q) (car (rgsubst beta (car minpoly) q)))
 349                        polys))
 350                 (setq minpoly
 351                       (cons *algvar (cdr nminpoly))))))))
 352
 353 (defun exchangevar (poly var)
 354   (let ((newvar (gensym))
 355         (ovar (car poly)))
 356     (setf (symbol-value newvar) (1+ (eval ovar)))
 357     (pvsubst ovar newvar
 358              (pvsubst var ovar
 359                       (pvsubst newvar var poly)))))
 360
 361 (defun rgsubst (val var p) ;;generalized psubst substitutes any
 362   (cond ((pcoefp p)
 363          (cons p 1)) ;;expression for any var in p
 364         ((eq var (car p))
 365          (cond ((pzerop val)
 366                 (cons (ptterm (cdr p) 0) 1))
 367                ((do ((ld (cadr p) (car a))
 368                      (a (cdddr p) (cddr a))
 369                      (ans (cons (caddr p) 1)
 370                           (ratplus
 371                            (rattimes ans
 372                                      (ratexpt val
 373                                               (- ld (car a)))
 374                                      t)
 375                            (cons (cadr a) 1))))
 376                     ((null a) (rattimes ans (ratexpt val ld) t))))))
 377         ((pointergp var (car p)) (cons p 1))
 378         (t (let ((newsym (gensym)))
 379              (setf (symbol-value newsym) (1+ (symbol-value (car p))))
 380              (rgsubst val newsym (pvsubst newsym var p))))))
 381
 382 (defun ratgsubst (val var rat)
 383   (ratquotient (rgsubst val var (car rat))
 384                (rgsubst val var (cdr rat))))
 385
 386 (defun algnorm (f p)
 387   (presult f p (car p)))
 388
 389 (defmfun $algnorm (r p var)
 390   (let ((varlist (list var))
 391         (genvar))
 392     (setq r (ratf r)
 393           p (cadr (ratf p)))
 394     (rdis* (cons (algnorm (cadr r) p)
 395                  (algnorm (cddr r) p)))))
 396
 397 (defun sqfrnorm (f p fvar) ;;f must be sqfr, p is minpoly, fvar # pvar
 398   (*bind* ((pvar (car p)))
 399           (setq f (cdr (ordervar pvar (list f p))) ;;new main var will be car of p
 400                 p (cadr f) f (car f))  ;make mainvar of f = mainvar(p)
 401           (do ((i 0 (1+ i))
 402                (dif (pdifference (make-poly fvar) (make-poly (car p))))
 403                (f f (pgsubst dif fvar f))
 404                (res))
 405               ((and (eq (car f) (car p))
 406                     (setq res (primpart (algnorm f p)))
 407                     (psqfrp res fvar))
 408                (list res
 409                      (*bind* (($algebraic t))           ;;;modified f
 410                              (putprop pvar (cdr p) 'tellrat)
 411                              (pvsubst pvar (car p) f))
 412                      (car p)
 413                      p
 414                      i)))))
 415
 416 (defun primelmt (a b gvar &aux ($algebraic nil))
 417   ;;a is a poly with coeff's in k(b)
 418   ;;gvar is new variable
 419   (let ((norm (sqfrnorm (cons gvar (cdr a)) b gvar))
 420         (alpha) (beta) ($ratalgdenom t))
 421     (rplaca norm (primpart (car norm)))
 422     (putprop gvar (cdar norm) 'tellrat)
 423     (setq $algebraic t
 424           beta (subresgcd (cadddr norm)
 425                      (pvsubst (caddr norm)
 426                               (car b)
 427                               (cadr norm))))
 428     (setq beta (plsolve beta)
 429           alpha (ratplus (cons (make-poly gvar) 1)
 430                          (rattimes (cons (- (cadddr (cdr norm))) 1)
 431                                    beta t)))
 432     (list (car norm) ;;minimal poly
 433           (pplus (make-poly (car a)) ;;new prim elm in old guys
 434                  (list (car b) 1  (car (last norm))))
 435           alpha beta))) ;;in terms of gamma
 436
 437 (defun $primelmt (f_b p_a c)
 438   ;;p_a(a) is an irreducible polynomial in K defining an extension
 439   ;;K[a] of degree n_a. Then f_b(b) is a polynomial in K[a] of degree
 440   ;;n_f , which defines a new extension K[a,b]. The output is a
 441   ;;polynomial of degree n_a*n_b, and the expression of its variable
 442   ;;in terms of a root a of p_a and a root b of f_b.
 443   ;;One assumes that p_a only depends on one variable a and f_b of two
 444   ;;different variables b and a. So this works only for algebraic
 445   ;;numbers, not algebraic curves.
 446
 447
 448
 449   (let* ((vla (newvar p_a))  ; ensure varlist = ($a $b)
 450          (vlf (newvar f_b))
 451          (varlist (cons (car vla) (remove (car vla) vlf)))
 452          (genvar)  ; start with a clean space of vars
 453          (genpairs)
 454          (p (cadr (ratf p_a)))          ;p_a in rat form
 455          (f (cadr (ratf f_b)))          ;f_b in rat form mainvar b
 456          (algvar (caar (newsym c)))  ;new primitive element
 457          (prim  (primelmt f p algvar)))
 458     (list '(mlist) (pdis (car prim)) (pdis (cadr prim)))))
 459           ;(rdis (caddr prim)) (rdis (cadddr prim))))) ;debug alpha beta
 460
 461
 462
 463 ;; discriminant of a basis
 464
 465 (defmfun $bdiscr (&rest args)
 466   (let ((varlist) (genvar))
 467     (xcons (bdiscr (mapcar #'rform (butlast args))
 468                    (car (rform (car (last args)))))
 469            (list 'mrat 'simp varlist genvar))))
 470
 471 (defun bdiscr (l minp)
 472   (det (mapcar #'(lambda (q1)
 473                    (mapcar #'(lambda (q2)
 474                                (algtrace (cons (ptimes (car q1)
 475                                                  (car q2)) 1)
 476                                          minp)) l)) l)))