share/contrib/ratpow.lisp

   1 ;;; ratpow version 0.3 4 June 2015
   2 ;;;
   3 ;;; Copyright 2013-2015 Stavros Macrakis
   4 ;;; Released under the LGPL (http://www.gnu.org/copyleft/lesser.html)
   5
   6 ;;; These functions operate on the *numerator* of a CRE
   7 ;;; to find the exponents in the denominator, use ratXXpow(ratdenom(...))
   8
   9 ;;; ratp_hipow -- highest power of the main variable in a CRE polynomial
  10 ;;; ratp_lopow -- lowest power of the main variable in a CRE polynomial
  11 ;;; ratp_coeffs -- list of powers and coefficients of the main variable
  12 ;;; ratp_dense_coeffs -- list of coefficients of the main variable, highest first
  13 ;;; ratp_dense_coeffs_lo -- list of coefficients of the main variable, lowest first
  14
  15 ;; To get a list of vars in a CRE, use showratvars
  16
  17 ;;; ratp_hipow( expr, var ) => highest power of var in ratnumer(expr)
  18 ;;;  ratp_hipow( x^(5/2) + x^2 , x) => 2
  19 ;;;  ratp_hipow( x^(5/2) + x^2 , sqrt(x)) => 5
  20 (defun $ratp_hipow (e v) (ratp_pow e v 'hipow))
  21
  22 ;;; ratp_lopow( expr, var ) => lowest power of var in ratnumer(expr)
  23 ;;;  ratp_lopow( x^(5/2) + x^3 , x) => 0
  24 ;;     CRE is {x}^3 + {x^(1/2)}^5 * {x}^0, where {} denotes kernels
  25 (defun $ratp_lopow (e v) (ratp_pow e v 'lopow))
  26
  27 ;;; ratp_coeffs( expr, var ) => list of powers/coefficients in ratnumer(expr)
  28 ;;; returned coefficients are in CRE form except for numbers
  29 ;;;  ratp_coeffs( 4*x^3 + x + sqrt(x), x) => [[3,4],[1,1],[0,sqrt(x)]]
  30 (defun $ratp_coeffs (e v) (ratp_pow e v 'sparse_coeffs))
  31
  32 ;;; ratp_dense_coeffs( expr, var ) => list of coefficients in ratnumer(expr), highest first
  33 ;;; returned coefficients are in CRE form except for numbers
  34 ;;;  ratp_dense_coeffs( 4*x^3 + x + sqrt(x), x) => [4,0,1,sqrt(x)]
  35 (defun $ratp_dense_coeffs (e v) (ratp_pow e v 'dense_coeffs))
  36
  37 ;;; ratp_dense_coeffs_lo( expr, var ) => list of coefficients in ratnumer(expr), lowest first
  38 ;;; returned coefficients are in CRE form except for numbers
  39 ;;;  ratp_dense_coeffs_lo( 4*x^3 + x + sqrt(x), x) => [sqrt(x),1,0,4]
  40 (defun $ratp_dense_coeffs_lo (e v) (ratp_pow e v 'dense_coeffs_low))
  41
  42
  43 (defun ratp_pow (e v func)
  44   (if (mbagp e)
  45       (simplify (cons (list (caar e))
  46                       (mapcar #'(lambda (i) (ratp_pow i v func)) (cdr e))))
  47     (let* ((e ($rat (if ($taylorp e) ($ratdisrep e) e) v)) ; change main variable if necessary
  48            (mrat (car e))
  49            (kernels (caddr mrat))
  50            (gensyms (cadddr mrat))
  51            (mainvar-gensym (nth (1- (length kernels)) gensyms))
  52            (numerat (cadr e))
  53            (res ())
  54            (lis
  55             (if (or (pcoefp numerat)
  56                     (not (eq (car numerat) mainvar-gensym)))
  57                 (list 0 numerat)
  58               (cdr numerat))))
  59       (cond ((eq func 'hipow) (car lis))
  60             ((eq func 'lopow) (car (last lis 2)))
  61             (t (setq kernels (butlast kernels 1))
  62                (setq mrat (list 'mrat 'simp kernels gensyms))
  63   ;;; Putprop needed for pdisrep, but not part-rat
  64   ;;;        (mapc #'(lambda (y z) (putprop y z 'disrep)) gensyms kernels)
  65                (cond
  66                 ((eq func 'sparse_coeffs)
  67                  (do ((lis lis (cddr lis)))
  68                      ((null lis)
  69                       (cons '(mlist simp) (nreverse res)))
  70                    (push (list '(mlist simp)
  71                                (car lis)
  72                                (part-rat mrat (cadr lis)))
  73                          res)))
  74                 ((memq func '(dense_coeffs dense_coeffs_low))
  75                  (do ((lis lis (cddr lis))
  76                       (pow (car lis) (1- pow)))
  77                      ((null lis)
  78                       (cons '(mlist simp)
  79                             (if (eq func 'dense_coeffs_low)
  80                                 res
  81                               (nreverse res))))
  82                    (while (> pow (car lis))
  83                      (push 0 res)
  84                      (setq pow (1- pow)))
  85                    (push (part-rat mrat (cadr lis)) res)))))))))
  86
  87 (defun part-rat (mrat expr)
  88   (if (pcoefp expr)
  89       expr
  90     (list* mrat expr 1)))