integration of atan2



Comments are welcomed; my second waffle:

 (1) delete the atan2 --> atan conversion in integrator (defined in sin.lisp)
 (2) append $atan2 to the list of functions risch attempts to integrate (append to rischform)
 (3) revised trigin1 (defined in risch.lisp)

Before and after:

Nounform :(

  (%i1) integrate(exp(cos(x)) * cos(sin(x)),x,0,2*%pi);
  (%o1) integrate(%e^cos(x)*cos(sin(x)),x,0,2*%pi)

Bogus (integrand is positive):

 (%i2) expand(integrate(atan2(1,x),x,-1,1));
 (%o2) -%pi

Invalid for x < 0:

 (%i3) integrate(atan2(x,x),x);
 (%o3) (%pi*x)/4

Experimental code (see below)

 (%i4) load("F:/larry.lisp")$
 (%i5) load("F:/sin.lisp")$

Correct, (I think) and not a nounform:

 (%i6) integrate(exp(cos(x)) * cos(sin(x)),x,0,2*%pi);
 (%o6) 2*%pi

Correct:

 (%i7) expand(integrate(atan2(1,x),x,-1,1));
 (%o7) %pi

OK on reals:

 (%i8) integrate(atan2(x,x),x);
 (%o8) x*atan2(x,x)

 (%i9) diff(%,x);
 (%o9) atan2(x,x)

--Barton

defun rischform (l)
  (cond ((or (atom l) (alike1 intvar l) (freeof intvar l)) nil)
    ((polylogp l)
     (if (and (integerp (car (subfunsubs l)))
          (signp g (car (subfunsubs l))))
         (rischform (car (subfunargs l)))
         (setq operator t)))
    ((atom (caar l))
     (case (caar l)
       ((%sin %cos %tan %cot %sec %csc)
        (setq trigint t $exponentialize t)
        (rischform (cadr l)))
       ((%asin %acos %atan %acot %asec %acsc $atan2) ; <--- append $atan2
        (setq trigint t $logarc t)
        (rischform (cadr l)))
       ((%sinh %cosh %tanh %coth %sech %csch)
        (setq hypertrigint t $exponentialize t)
        (rischform (cadr l)))
       ((%asinh %acosh %atanh %acoth %asech %acsch)
        (setq hypertrigint t $logarc t)
        (rischform (cadr l)))
       ((mtimes mplus mexpt rat %erf %log)
        (mapc #'rischform (cdr l)))
       (t (setq operator (caar l)))))
    (t (setq operator (caar l)))))


(defun attempt-analytic-continuation (e)
  (let (($opsubst t))
    ;; atan2(cos(x),sin(x)) --> x
    (setq e ($substitute #'(lambda (a b)
                 (if (and (consp a) (consp (car a)) (eq (mop a) '%cos)
                      (consp b) (consp (car b)) (eq (mop b) '%sin)
                      (alike (second a) (second b)))
                 (second a) (take '($atan2) a b))) '$atan2 e))
    ;; atan(tan(x)) --> x
    (setq e (let (($triginverses '$all) ($exponentialize nil) ($demoivre t)) ($trigreduce e))) ; trigreduce does sin(x)/cos(x) -> tan(x).

    ;; Convert incomplete_gamma(0,x) to exponential_ei functions; and convert
    ;; exponential_ei(z) --> %gamma + (log(z) - log(1/z))/2 + z * hypergeometric([1,1],[2,2],z).
    (let (($gamma_expand t) (dosimp t) ($exponentialize t) ($logexpand '$all))
      (setq e (sratsimp (resimplify e))) ; gotta resimplify to convert from incomplete_gamma to expintegral_ei.
      (setq e ($substitute #'(lambda (z)
                   (add '%gamma
                    (div (sub (take '(%log) z) (take '(%log) (div 1 z))) 2)
                    (mul z (take '($hypergeometric) (list '(mlist) 1 1) (list '(mlist) 2 2) z)))) '%expintegral_ei e))
      (sratsimp e))))


(defun trigin1 (*exp var)
  (let ((rischp var) (rp-polylogp t) ($logarc nil) ($exponentialize nil)
    (result (hypertrigint1 *exp var nil)))
    ;; When result has a integral nounform, set *in-risch-p* to true and resimplify.
    ;; Presumably setting *in-risch-p*to true avoids an infinite loop.
    (if (isinop result '%integrate)
    (let ((*in-risch-p* t)) (setq result (resimplify result))))

    (setq result (let ((generate-atan2 t)) ($rectform result)))
    (sratsimp (attempt-analytic-continuation (sratsimp result)))))

Replaces:


(defun trigin1 (*exp var)

  (let ((yyy (hypertrigint1 *exp var nil)))

    (setq yyy (div ($expand ($num yyy))

                   ($expand ($denom yyy))))

    (let ((rischp var) (rp-polylogp t) $logarc $exponentialize result)

      (setq result (sratsimp (if (and (freeof '$%i *exp) (freeof '$li yyy))

                                 ($realpart yyy)

                                 ($rectform yyy))))

      ;; The result can contain solveable integrals. Look for this case.

      (if (isinop result '%integrate)

          ;; Found an integral. Evaluate the result again.

          ;; Set the flag *in-risch-p* to make sure that we do not call

          ;; rischint again from the integrator. This avoids endless loops.

          (let ((*in-risch-p* t))

            (meval (list '($ev) result '$nouns)))

          result))))




________________________________