Menu

#4014 defint() produces incorrect output

None
open
nobody
None
5
2022-08-09
2022-08-08
No
assume(sin(phi) > 0);
f: 1/(w^4 - 2*cos(2*phi) * w^2 + 1);
defint(f, w, -inf, inf);

produces incorrect output.
presumably this would be correct

gold:  %pi/(2*sin(phi));

walking through the code it turns out the variabe "*semirat*" which is used in "(defun polelist ...)" is responsible for this behaviour. It is set in defint.lisp and will cause polelist to return unwanted poles outside the given "region".

If I disable this in residue.lisp

@@ -44,30 +44,31 @@
 ;; Finally, the fourth part is NIL, unless *semirat* is T.
 (defun polelist (d region region1)
   (prog (roots $breakup r rr ss r1 s pole wflag cf)
      (setq wflag t)
      (setq leadcoef (polyinx d var 'leadcoef))
      (setq roots (solvecase d))
      (if (eq roots 'failure) (return ()))
      ;; Loop over all the roots.  SOLVECASE returns the roots in the form
      ;; ((x = r1) mult1
      ;;  (x = r2) mult2
      ;;  ...)

    loop1
      (cond ((null roots)
        (cond ((and *semirat*
+                        nil
            (> (+ (length s) (length r))
               (+ (length ss) (length rr))))
           ;; Return CF, repeated roots (*semirat*), simple
           ;; roots (*semirat*), roots in region 1.
           (return (list cf rr ss r1)))
          (t
           ;; Return CF, repeated roots, simple roots, roots in region 1.
           (return (list cf r s r1)))))
       (t
        ;; Set POLE to the actual root and set D to the
        ;; multiplicity of the root.
        (setq pole (caddar roots))
        (setq d (cadr roots))
        (cond (leadcoef
           ;; Is it possible for LEADCOEF to be NIL ever?

then this

assume(sin(phi) > 0);
f: 1/(w^4 - 2*cos(2*phi) * w^2 + 1);
defint(f, w, -inf, inf);
trigsimp(demoivre(factor(exponentialize(%)));

produces the expected result.

Discussion

  • Robert Dodier

    Robert Dodier - 2022-08-09

    Robert, thanks very much for investigating. After making that change, what result do you get from run_testsuite() ? Are any of the results different? If no existing correct results are made incorrect, then that's a good sign.

     
    • Robert Larice

      Robert Larice - 2022-08-09

      Hello,
      I forgot to mention, there is no change of importance in test-suite.log.
      I stumbled on comments from rtoy concerning this semirat variable, saying he didn't know what it is good for. Yet I obviously can't simply remove it, my "nil" is just a way to point at some point. I traced through zmtorat etc etc, and guessed the code tried do integrate by residues. Then I spotted incorrect residues and fell into "polelist". skimming through it one cant avoid guessing semirat causing to return exact the opposite of what would have been returned without it. And voila, disabling it the correct result is produced. From there to a proper bug fix is propably quite a bit more difficult.
      I attach a file in the style of the other "test.mac" files, which I finally used to check the output.

       

Log in to post a comment.