search

rtoy / maxima

A Clone of Maxima's repo
Other
0 stars 0 forks source link

defint() produces incorrect output #3920

Open rtoy opened 1 month ago

rtoy commented 1 month ago

Imported from SourceForge on 2024-07-08 21:17:15 Created by rlar on 2022-08-08 17:31:36 Original: https://sourceforge.net/p/maxima/bugs/4014

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.

rtoy commented 1 month ago

Imported from SourceForge on 2024-07-08 21:17:16 Created by robert_dodier on 2022-08-09 15:15:35 Original: https://sourceforge.net/p/maxima/bugs/4014/#d33c

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.

rtoy commented 1 month ago

Imported from SourceForge on 2024-07-08 21:17:20 Created by rlar on 2022-08-09 19:18:34 Original: https://sourceforge.net/p/maxima/bugs/4014/#d33c/8ea7

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.

Attachments: