packages feed

egison-3.7.11: lib/math/algebra/root.egi

;;;;;
;;;;;
;;;;; Algebra
;;;;;
;;;;;

;;;
;;; Root
;;;

(define $rt
  (lambda [$n $x]
    (if (integer? n)
      (match x math-expr
        {[,0 0]
         [?monomial? (rt-monomial n x)]
         [<div <poly $xs> <poly $ys>>
          (let {[$xd (reduce gcd xs)]
                [$yd (reduce gcd ys)]}
            (let {[$d (rt-monomial n (/ xd yd))]}
              (*' d
                 (rt'' n (*' (/' (sum' (map (/' $ xd) xs)) (sum' (map (/' $ yd) ys)))))
                 )))]
         [_ (rt'' n x)]})
      (rt'' n x))))

(define $rt-monomial
  (lambda [$n $x]
    (/ (rt-term n (* (numerator x)
                     (** (denominator x) (- n 1))))
       (denominator x))))

(define $rt-term
  (lambda [$n $x]
    (match x term-expr
      {[<term $a _>
        (if (lt? a 0)
          (*' (rtm1 n) (rt-positive-term n (* -1 x)))
          (rt-positive-term n x))]})))

(define $rt-positive-term
  (lambda [$n $x]
    (match [n x] [math-expr math-expr]
      {[[,3 (* $a ,i $r)] (* -1 i (rt 3 (*' a r)))]
       [[_ (* $a (,sqrt $b) $r)] (*' (rt (* n 2) (*' (**' a 2) b)) (rt n r))]
       [[_ (* $a (,rt $n' $b) $r)] (*' (rt (* n n') (*' (**' a n') b)) (rt n r))]
       [[_ _] (rt-positive-term1 n x)]
       })))

(define $rt-positive-term1
  (lambda [$n $x]
    (letrec {[$f (lambda [$xs]
                   (match xs (assoc-multiset math-expr)
                     {[<nil> [1 1]]
                      [<ncons $p $k $rs>
                       (let {[[$a $b] (f rs)]}
                         [(*' (**' p (quotient k n)) a) (*' (**' p (remainder k n)) b)])]}))]
             [$g (lambda [$n $x]
                   (let {[$d (match x term-expr
                               {[<term $m $xs> (gcd n (reduce gcd (map 2#%2 {@(to-assoc (p-f m)) @xs})))]})]}
                     (rt'' (/ n d) (rt d x))))]}
      (match x term-expr
        {[<term $m $xs>
          (match (f {@(to-assoc (p-f (abs m))) @xs}) [integer integer]
            {[[$a ,1] a]
             [[$a $b] (*' a (g n b))]})]}))))

(define $rt''
  (lambda [$n $x]
    (match [n x] [integer integer]
      {[[,2 _] (`sqrt x)]
       [[_ _] (`rt n x)]})))

(define $rtm1
  (lambda [$n]
    (match n integer
      {[,1 -1]
       [,2 i]
       [?odd? -1]
       [_ undefined]})))

(define $sqrt
  (lambda [$x]
    (if (scalar? x)
      (let {[$m (numerator x)]
            [$n (denominator x)]}
        (/ (rt 2 (* m n)) n))
      (b.sqrt x))))

(define $rt-of-unity rtu)

(define $rtu
  (lambda [$n]
    (rtu' n)))

(define $rtu'
  (lambda [$n]
    (if (integer? n)
      (match n integer
        {[,1 1]
         [,2 -1]
         [,3 w]
         [,4 i]
         [_ (`rtu n)]
         })
      (`rtu n))))