src/rat3c.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancements.                   ;;;;;
   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 (macsyma-module rat3c)
  14
  15 ;;      THIS IS THE NEW RATIONAL FUNCTION PACKAGE PART 3.
  16 ;;      IT INCLUDES THE GCD ROUTINES AND THEIR SUPPORTING FUNCTIONS
  17
  18 (defun cgcd (a b)
  19   (cond (modulus 1)
  20         ((and $keepfloat (or (floatp a) (floatp b))) 1)
  21         (t (gcd a b))))
  22
  23 (defun pquotientchk (a b)
  24   (if (equal b 1) a (pquotient a b)))
  25
  26 ;; divides polynomial x by polynomial y
  27 ;; avoids error "quotient by polynomial of higher degree"
  28 ;;  (returns nil in this case)
  29 (defun pquotientchk-safe (x y)
  30   (ignore-rat-err (pquotientchk x y)))
  31
  32 (defun ptimeschk (a b)
  33   (cond ((equal a 1) b)
  34         ((equal b 1) a)
  35         (t (ptimes a b))))
  36
  37 (defun pfloatp (x)
  38   (catch 'float (if (pcoefp x) (floatp x) (pfloatp1 x))))
  39
  40 (defun pfloatp1 (x)
  41   (mapc #'(lambda (q) (cond ((pcoefp q) (when (floatp q) (throw 'float t)))
  42                             ((pfloatp1 q))))
  43         (cdr x))
  44   nil)
  45
  46 (defun pgcd (x y)
  47   (setq x (car (pgcda x y nil)))
  48   (cond ((pminusp x) (pminus x))
  49         (modulus (monize x))
  50         (t x)))
  51
  52 (defun plcm (x y)
  53   (setq x (pgcdcofacts x y))
  54   (ptimes (car x) (ptimes (cadr x) (caddr x))))
  55
  56 (defun plcmcofacts (x y)
  57   (setq x (pgcdcofacts x y))
  58   (list (ptimes (car x) (ptimes (cadr x) (caddr x)))
  59         (caddr x) (cadr x)))
  60
  61 ; returns list (gcd xx yy alg)
  62 ; where x * y = gcd^2 * xx * yy / alg^2
  63 ; and alg is non-nil only when $algebraic is true
  64 (defun pgcdcofacts (x y)
  65   (let ((a (pgcda x y t)))
  66     (cond ((cdr a) a)
  67           ((equal (setq a (car a)) 1) (list 1 x y))
  68           ((and $algebraic (not (pcoefp a)))
  69            (cons a (prog2 (setq x (rquotient x a)
  70                                 y (rquotient y a)
  71                                 a (pgcdcofacts (cdr x) (cdr y)))
  72                        (list (ptimes (car x) (caddr a))
  73                              (ptimes (car y) (cadr a))
  74                              (ptimes (cadr a) (cdr y))))))
  75           ((eq a x) (list x 1 (pquotient y x)))
  76           ((eq a y) (list a (pquotient x y) 1))
  77           (t (list a (pquotient x a) (pquotient y a))))))
  78
  79 (defun pgcda (x y cofac? &aux a c)
  80   (cond ((not $gcd) (list 1 x y))
  81         ((and $keepfloat (or (pfloatp x) (pfloatp y)))
  82          (cond ((or (pcoefp x) (pcoefp y)
  83                     (pcoefp (setq a (car (ptermcont x))))
  84                     (pcoefp (setq a (pgcd a (car (ptermcont y))))))
  85                 (list 1 x y))
  86                (t (list a))))
  87         ((pcoefp x)
  88          (cond ((pcoefp y)
  89                 (cons (setq a (cgcd x y))
  90                       (and cofac?
  91                            (list (cquotient x a) ;(CQUOTIENT 0 0) = 0
  92                                  (cquotient y a)))))
  93                ((zerop x) (list y x 1))
  94                (t (list (pcontent1 (cdr y) x)))))
  95         ((pcoefp y) (cond ((zerop y) (list x 1 y))
  96                           (t (list (pcontent1 (cdr x) y)))))
  97         ((equal x y) (list x 1 1))
  98         ($ratfac (multiple-value-list (fpgcdco x y)))
  99         ((not (eq (p-var x) (p-var y)))
 100          (list (if (pointergp (p-var x) (p-var y))
 101                    (oldcontent1 (p-terms x) y)
 102                    (oldcontent1 (p-terms y) x))))
 103         ((progn (desetq (a x) (ptermcont x))
 104                 (desetq (c y) (ptermcont y))
 105                 (not (and (equal a 1) (equal c 1))))
 106          (mapcar #'ptimes (monomgcdco a c cofac?) (pgcda x y cofac?)))
 107         ((and (not $algebraic) (not modulus)
 108               (desetq (a . c) (lin-var-find (nreverse (pdegreevector x))
 109                                             (nreverse (pdegreevector y))
 110                                             (reverse genvar))))
 111          (cond ((= a 1) (linhack x y (car c) (cadr c) cofac?))
 112                (t (setq a (linhack y x a (cadr c) cofac?))
 113                   (if (cdr a) (rplacd a (nreverse (cdr a))))
 114                   a)))
 115         ((eq $gcd '$spmod) (list (zgcd x y)))
 116         ((eq $gcd '$subres) (list (oldgcd x y)))
 117         ((eq $gcd '$algebraic)
 118          (if (or (palgp x) (palgp y))
 119              (let (($gcd '$subres)) (list (oldgcd x y)))
 120              (let (($gcd '$spmod)) (list (zgcd x y)))))
 121         ((eq $gcd '$ez) (ezgcd2 x y))
 122         ((eq $gcd '$red) (list (oldgcd x y)))
 123         ((eq $gcd '$mod) (newgcd x y modulus))
 124         ((not (member $gcd *gcdl* :test #'eq))
 125          (merror (intl:gettext "gcd: 'gcd' variable must be one of ~M; found: ~M") *gcdl* $gcd))
 126         (t (list 1 x y))))
 127
 128 (defun monomgcdco (p q cofac?)
 129   (let ((gcd (monomgcd p q)))
 130     (cons gcd (if cofac? (list (pquotient p gcd) (pquotient q gcd)) ()))))
 131
 132 (defun monomgcd (p q)
 133   (cond ((or (pcoefp p) (pcoefp q)) 1)
 134         ((eq (p-var p) (p-var q))
 135          (make-poly (p-var p) (min (p-le p) (p-le q))
 136                     (monomgcd (p-lc p) (p-lc q))))
 137         ((pointergp (car p) (car q)) (monomgcd (p-lc p) q))
 138         (t (monomgcd p (p-lc q)))))
 139
 140 (defun linhack (pol1 pol2 nonlindeg var cofac?)
 141   (prog (coeff11 coeff12 gcdab rpol1 rpol2 gcdcd gcdcoef)
 142      (desetq (coeff11 . coeff12) (bothprodcoef (make-poly var) pol1))
 143      (setq gcdab (if (pzerop coeff12) coeff11
 144                      (pgcd coeff11 coeff12)))
 145      (cond ((equal gcdab 1)
 146             (cond ((setq coeff11 (testdivide pol2 pol1))
 147                    (return (list pol1 1 coeff11)))
 148                   (t (return (list 1 pol1 pol2))))))
 149      (setq rpol1 (pquotient pol1 gcdab))
 150      (desetq (gcdcd rpol2) (linhackcontent var pol2 nonlindeg))
 151      (cond ((equal gcdcd 1)
 152             (cond ((setq coeff12 (testdivide rpol2 rpol1))
 153                    (return (list rpol1 gcdab coeff12)))
 154                   (t (return (list 1 pol1 pol2))))))
 155      (cond (cofac? (desetq (gcdcoef coeff11 coeff12)
 156                            (pgcdcofacts gcdab gcdcd))
 157                    (cond ((setq gcdcd (testdivide* rpol2 rpol1))
 158                           (return (list (ptimes gcdcoef rpol1)
 159                                         coeff11
 160                                         (ptimes coeff12 gcdcd))))
 161                          (t (return (list gcdcoef
 162                                           (ptimes coeff11 rpol1)
 163                                           (ptimes coeff12 rpol2))))))
 164            (t (setq gcdcoef (pgcd gcdcd gcdab))
 165               (cond ((testdivide* rpol2 rpol1)
 166                      (return (list (ptimes gcdcoef rpol1))))
 167                     (t (return (list gcdcoef))))))))
 168
 169 (defun lin-var-find (a b c)
 170   (do ((varl c (cdr varl))
 171        (degl1 a (cdr degl1))
 172        (degl2 b (cdr degl2)))
 173       ((or (null degl1) (null degl2)) nil)
 174     (if (equal (min (car degl1) (car degl2)) 1)
 175         (return (list (car degl1) (car degl2) (car varl))))))
 176
 177 (defun linhackcontent (var pol nonlindeg &aux (npol pol) coef gcd)
 178   (do ((i nonlindeg (1- i)))
 179       ((= i 0) (list (setq gcd (pgcd gcd npol)) (pquotient pol gcd)))
 180     (desetq (coef . npol) (bothprodcoef (make-poly var i 1) npol))
 181     (unless (pzerop coef)
 182       (setq gcd (if (null gcd) coef (pgcd coef gcd)))
 183       (if (equal gcd 1) (return (list 1 pol))))))
 184
 185 ;;*** THIS IS THE REDUCED POLYNOMIAL REMAINDER SEQUENCE GCD (COLLINS')
 186
 187 (defun oldgcd (x y &aux u v s egcd)     ;only called from pgcda
 188   (desetq (x  u) (oldcontent x))
 189   (desetq (y  v) (oldcontent y))
 190   (setq egcd (gcd (pgcdexpon u) (pgcdexpon v)))
 191   (if (> egcd 1)
 192       (setq u (pexpon*// u egcd nil)
 193             v (pexpon*// v egcd nil)))
 194   (if (> (p-le v) (p-le u)) (rotatef u v))
 195   (setq s (case $gcd
 196             ($red (redgcd u v))
 197             ($subres (subresgcd u v))
 198             (t (merror "OLDGCD: found gcd = ~M; how did that happen?" $gcd))))
 199   ;; Check for gcd that simplifies to 0. SourceForge bugs 831445 and 1313987
 200   (unless (ignore-rat-err (rainv s))
 201     (setq s 1))
 202   (unless (equal s 1)
 203     (setq s (pexpon*// (primpart
 204                         (if $algebraic s
 205                             (pquotient s (pquotient (p-lc s)
 206                                                     (pgcd (p-lc u) (p-lc v))))))
 207                        egcd t)))
 208   (setq s (ptimeschk s (pgcd x y)))
 209   (and $algebraic (not (pcoefp (setq u (leadalgcoef s))))
 210        (not (equal u s)) (setq s (algnormal s)))
 211   (cond (modulus (monize s))
 212         ((pminusp s) (pminus s))
 213         (t s)))
 214
 215 (defun pgcdexpon (p)
 216   (if (pcoefp p) 0
 217       (do ((d (cadr p) (gcd d (car l)))
 218            (l (cdddr p) (cddr l)))
 219           ((or (null l) (= d 1)) d))))
 220
 221 (defun pexpon*// (p n *?)
 222   (if (or (pcoefp p) (= n 1)) p
 223       (do ((ans (list (car p))
 224                 (cons (cadr l)
 225                       (cons (if *? (* (car l) n)
 226                                 (truncate (car l) n))
 227                             ans)))
 228            (l (cdr p) (cddr l)))
 229           ((null l) (nreverse ans)))))
 230
 231 ;;polynomial gcd using reduced prs
 232
 233 (defun redgcd (p q &aux (d 0))
 234   (loop until (zerop (pdegree q (p-var p)))
 235          do (psetq p q
 236                    q (pquotientchk-safe (prem p q) (pexpt (p-lc p) d))
 237                    d (+ (p-le p) 1 (- (p-le q))))
 238          (if (< d 1) (return 1))
 239          finally (return (if (pzerop q) p 1))))
 240
 241 ;;computes gcd's using subresultant prs
 242 ;;ACM Transactions On Mathematical Software Sept. 1978
 243
 244 (defun subresgcd (p q)
 245   (loop for g = 1 then (p-lc p)
 246          for h = 1 then (pquotientchk-safe (pexpt g d) h^1-d)
 247          for d = (- (p-le p) (p-le q))
 248          for h^1-d = 1 then (if (< d 1)
 249                                 (return 1)
 250                               (pexpt h (1- d)))
 251          do (psetq p q
 252                    q (pquotientchk-safe (prem p q) (ptimes g (ptimes h h^1-d))))
 253          if (zerop (pdegree q (p-var p))) return (if (pzerop q) p 1)))
 254
 255 ;;*** THIS COMPUTES PSEUDO REMAINDERS
 256
 257 (defun psquorem1 (u v quop)
 258   (prog (k (m 0) lcu lcv quo lc)
 259      (declare (special lcu lcv))
 260      (setq lcv (pt-lc v))
 261      (setq k (- (pt-le u) (pt-le v)))
 262      (cond ((minusp k) (return (list 1 '(0 0) u))))
 263      (if quop (setq lc (pexpt (pt-lc v) (1+ k))))
 264      a     (setq lcu (pminus (pt-lc u)))
 265      (if quop (setq quo (cons (ptimes (pt-lc u) (pexpt (pt-lc v) k))
 266                               (cons k quo))))
 267      (cond ((null (setq u (pgcd2 (pt-red u) (pt-red v) k)))
 268             (return (list lc (nreverse quo) '(0 0))))
 269            ((minusp (setq m (- (pt-le u) (pt-le v))))
 270             (setq u (cond ((zerop k) u)
 271                           (t (pctimes1 (pexpt lcv k) u))))
 272             (return (list lc (nreverse quo) u)))
 273            ((> (1- k) m)
 274             (setq u (pctimes1 (pexpt lcv (- (1- k) m)) u))))
 275      (setq k m)
 276      (go a)))
 277
 278 (defun prem (p q)
 279   (cond ((pcoefp p)
 280          (if (pcoefp q)
 281              (if (or modulus (floatp p) (floatp q))
 282                  0
 283                  (rem p q))
 284              p))
 285         ((pcoefp q) (pzero))
 286         (t (psimp (p-var p) (pgcd1 (p-terms p) (p-terms q))))))
 287
 288 (defun pgcd1 (u v) (caddr (psquorem1 u v nil)))
 289
 290 (defun pgcd2 (u v k &aux (i 0))
 291   (declare (special lcu lcv) (fixnum i))
 292   (cond ((null u) (pcetimes1 v k lcu))
 293         ((null v) (pctimes1 lcv u))
 294         ((zerop (setq i (+ (pt-le u) (- k) (- (car v)))))
 295          (pcoefadd (pt-le u) (pplus (ptimes lcv (pt-lc u))
 296                                     (ptimes lcu (pt-lc v)))
 297                    (pgcd2 (pt-red u) (pt-red v) k)))
 298         ((minusp i)
 299          (list* (+ (pt-le v) k) (ptimes lcu (pt-lc v)) (pgcd2 u (pt-red v) k)))
 300         (t (list* (pt-le u) (ptimes lcv (pt-lc u)) (pgcd2 (pt-red u) v k)))))
 301
 302 ;;;*** OLDCONTENT REMOVES ALL BUT MAIN VARIABLE AND PUTS THAT IN CONTENT
 303 ;;;***  OLDCONTENT OF 3*A*X IS 3*A (WITH MAINVAR=X)
 304
 305 (defun rcontent (p)                     ;RETURNS RAT-FORMS
 306   (let ((q (oldcontenta p)))
 307     (list (cons q 1) (cond ($algebraic (rquotient p q))
 308                            (t (cons (pquotient p q) 1))))))
 309
 310 (defun oldcontenta (x)
 311   (cond ((pcoefp x) x)
 312         (t (setq x (contsort (cdr x)))
 313            (oldcontent2 (cdr x) (car x)))))
 314
 315 (defun oldcontent (x)
 316   (cond ((pcoefp x) (list x 1))
 317         ((null (p-red x))
 318          (list (p-lc x) (make-poly (p-var x) (p-le x) 1)))
 319         (t (let ((u (contsort (cdr x))) v)
 320              (setq u (oldcontent2 (cdr u) (car u))
 321                    v (cond ($algebraic (car (rquotient x u)))
 322                            (t (pcquotient x u))))
 323              (cond ((pminusp v) (list (pminus u) (pminus v)))
 324                    (t (list u v)))))))
 325
 326 (defun oldcontent1 (x gcd)
 327   (cond ((equal gcd 1) 1)
 328         ((null x) gcd)
 329         (t (oldcontent2 (contsort x) gcd))))
 330
 331 (defun oldcontent2 (x gcd)
 332   (do ((x x (cdr x))
 333        (gcd gcd (pgcd (car x) gcd)))
 334       ((or (null x) (equal gcd 1)) gcd)))
 335
 336 (defun contsort (x)
 337   (setq x (coefl x))
 338   (cond ((member 1 x) '(1))
 339         ((null (cdr x)) x)
 340         (t (sort x #'contodr))))
 341
 342 (defun coefl (x)
 343   (do ((x x (cddr x))
 344        (ans nil (cons (cadr x) ans)))
 345       ((null x) ans)))
 346
 347 (defun contodr (a b)
 348   (cond ((pcoefp a) t)
 349         ((pcoefp b) nil)
 350         ((eq (car a) (car b)) (not (> (cadr a) (cadr b))))
 351         (t (pointergp (car b)(car a)))))
 352
 353 ;;;*** PCONTENT COMPUTES INTEGER CONTENT
 354 ;;;*** PCONTENT OF 3*A*X IS 3 IF MODULUS = NIL  1 OTHERWISE
 355
 356 (defun pcontent (x)
 357   (cond ((pcoefp x) (list x 1))
 358         (t (let ((u (pcontentz x)))
 359              (if (equal u 1) (list 1 x)
 360                  (list u (pcquotient x u)))))))
 361
 362 (defun pcontent1 (x gcd)
 363   (do ((x x (cddr x))
 364        (gcd gcd (cgcd gcd (pcontentz (cadr x)))))
 365       ((or (null x) (equal gcd 1)) gcd)))
 366
 367 (defun pcontentz (p)
 368   (cond ((pcoefp p) p)
 369         (t (pcontent1 (p-red p) (pcontentz (p-lc p))))))
 370
 371 (defun ucontent (p)                     ;CONTENT OF UNIV. POLY
 372   (cond ((pcoefp p) (abs p))
 373         (t (setq p (mapcar #'abs (coefl (cdr p))))
 374            (let ((m (reduce #'min p)))
 375              (oldcontent2 (delete m p :test #'equal) m)))))
 376
 377 ;;***   PGCDU CORRESPONDS TO BROWN'S ALGORITHM U
 378
 379 ;;;PGCDU IS NOT NOW IN RAT;UFACT >
 380
 381 (defun pgcdu (p q)
 382   (do () ((pzerop q) (monize p))
 383     (psetq p q q (pmodrem p q))))
 384
 385 (defun pmodrem (x y)
 386   (cond ((null modulus)
 387          (merror "PMODREM: null modulus; how did that happen?"))
 388         ((pacoefp y) (if (pzerop y) x 0))
 389         ((pacoefp x) x)
 390         ((eq (p-var x) (p-var y))
 391          (psimp (car x) (pgcdu1 (p-terms x) (p-terms y) nil)))
 392         (t (merror "PMODREM: I can't handle this; x = ~M, y = ~M" x y))))
 393
 394 (defun pmodquo (u v &aux quo)
 395   (declare (special quo))
 396   (cond ((null modulus)
 397          (merror "PMODQUO: null modulus; how did that happen?"))
 398         ((pcoefp v) (cons (ptimes (crecip v) u) 0))
 399         ((alg v) (cons (ptimes (painvmod v) u) 0))
 400         ((pacoefp u) (cons 0 u))
 401         ((not (eq (p-var u) (p-var v)))
 402          (merror "PMODQUO: arguments have different variables; how did that happen?"))
 403         (t (xcons (psimp (car u) (pgcdu1 (cdr u) (cdr v) t))
 404                   (psimp (car u) quo)))))
 405
 406
 407 (defun pgcdu1 (u v pquo*)
 408   (let ((invv (painvmod (pt-lc v))) (k 0) q*)
 409     (declare (special k quo q*))
 410     (loop until (minusp (setq k (- (pt-le u) (pt-le v))))
 411            do (setq q* (ptimes invv (pt-lc u)))
 412            if pquo* do (setq quo (nconc quo (list k q*)))
 413            when (ptzerop (setq u (ptpt-subtract-powered-product
 414                                   (pt-red u) (pt-red v) q* k)))
 415            return (ptzero)
 416            finally (return u))))
 417
 418 (defun newprime (p)
 419   (do ((pl *bigprimes* (cdr pl)))
 420       ((null pl)
 421        (setq p (next-prime (1- p) -1))
 422        (setq *bigprimes* (nconc *bigprimes* (list p)))
 423        p)
 424     (when (< (car pl) p)
 425       (return (car pl)))))
 426
 427 (defun leadcoefficient (p)
 428   (if (pcoefp p) p (leadcoefficient (caddr p))))
 429
 430 (defun maxcoefficient (p)
 431   (if (pcoefp p) (abs p) (maxcoef1 (cdr p))))
 432
 433 (defun maxcoef1 (p)
 434   (if (null p) 0 (max (maxcoefficient (cadr p)) (maxcoef1 (cddr p)))))
 435
 436 (defun maxnorm (poly)
 437   (if (null poly) 0 (max (norm (cadr poly)) (maxnorm (cddr poly)))))
 438
 439 (defun norm (poly)
 440   (cond ((null poly) 0)
 441         ((pcoefp poly) (abs poly))
 442         (t (+ (norm (caddr poly)) (norm1 (cdddr poly)) )) ))
 443
 444 (defun norm1 (poly)
 445   (if (null poly) 0 (+ (norm (cadr poly)) (norm1 (cddr poly)) )) )
 446
 447 (defun pdegree (p var)
 448   (cond ((pcoefp p) 0)
 449         ((eq var (p-var p)) (p-le p))
 450         ((pointergp var (p-var p)) 0)
 451         (t (do ((l (p-red p) (pt-red l))
 452                 (e (pdegree (p-lc p) var) (max e (pdegree (pt-lc l) var))))
 453                ((null l) e)))))
 454
 455 (defun poly-in-var (p v)
 456   (cond ((or (pcoefp p) (pointergp v (p-var p))) (list 0 p))
 457         ((eq (p-var p) v) (p-terms p))
 458         ((loop with ans
 459                 for (exp coef) on (p-terms p) by #'cddr
 460                 do (setq ans (ptptplus ans
 461                                        (everysubst2 (poly-in-var coef v)
 462                                                     (list (p-var p) exp 1))))
 463                 finally (return ans)))))
 464
 465 (defun univar (x)
 466   (or (null x) (and (pcoefp (pt-lc x)) (univar (pt-red x)))))
 467
 468 ;;**THE CHINESE REMAINDER ALGORITHM IS A SPECIAL CASE OF LAGRANGE INTERPOLATION
 469
 470 (defun lagrange3 (u uk p qk)
 471   (set-modulus p)
 472   (setq uk (pdifference uk (pmod u)))
 473   (cond ((pzerop uk) (setq modulus nil) u)
 474         (t (setq uk (pctimes (crecip (cmod qk)) uk))
 475            (setq modulus nil)
 476            (pplus u (pctimes qk uk)))))
 477
 478
 479 (defun lagrange33 (u uk qk xk)
 480   (declare (special xv))
 481   (setq uk (pdifference uk (pcsubst u xk xv)))
 482   (cond ((pzerop uk) u)
 483         (t (pplus u (ptimes
 484                      (pctimes (crecip (pcsubst qk xk xv)) uk)
 485                      qk)))))
 486
 487
 488 ;;;*************************************************************
 489
 490 ;;      THIS IS THE END OF THE NEW RATIONAL FUNCTION PACKAGE PART 3.
 491 ;;      IT INCLUDES THE GCD ROUTINES AND THEIR SUPPORTING FUNCTIONS