[Axiom-developer] Complex argument

[Waldek Hebisch]

Current version of complex argument (in gaussian.spad.pamphlet) is
responsible for multiple bugs (15, 47, 184, 293, 314 and a few wrong
answers in mapleok.output).  The patch below uses more complicated
but correct formula for complex argument.

The existing formula assumed that x has positive real part.  But
for integration it is typical to have constant imaginary part
and real part which changes sign. So the existing formula produced
spurious jump discontinuities.

With the patch Axiom gives now correct answers for most such problems.
Unfortunatly, correct answer is more complicated and ATM Axiom is
unable to simplify it.  Because of this some test outputs (especially
in mapleok.output) got bigger -- I am still checking if all new
results correct.

In principle we could try to determine when the old formula works
correctly -- we should be able to determine in which a constant
lives and for many functions it is also possible to prove that
the function has constant sign.  OTOH such machinery belongs to
other packages.  In fact we have various sign finding utilities
which work on expressions, but:
- sign finding may be pretty expensive
- it is not clear if such utilities are applicable to all
  domains allowed as arguments to complex (which probably
  means that we need to add a maze of conditionals choosing
  correct version of 'argument') 

Also, currently integrator converts exponentials and logarithms
to trigonometric and arc trigonometric functions by taking
real part of the answer -- this looks wrong for me.

Anyway, the patch follows:

diff -u wh-sandbox2.bb/src/algebra/gaussian.spad.pamphlet 
wh-sandbox2/src/algebra/gaussian.spad.pamphlet
--- wh-sandbox2.bb/src/algebra/gaussian.spad.pamphlet   2006-12-03 
04:23:11.000000000 +0100
+++ wh-sandbox2/src/algebra/gaussian.spad.pamphlet      2007-01-03 
18:01:06.000000000 +0100
@@ -367,10 +367,18 @@
            argument x == atan2loc(imag x, real x)
 
          else
-           -- Not ordered so dictate two quadrants
-           argument x ==
-             zero? real x => pi()$R * half
-             atan(imag(x) * recip(real x)::R)
+           if R has RadicalCategory then
+             argument x ==
+               n1 := sqrt(norm(x))
+               x1 := real(x) + n1
+               (2::R)*atan(imag(x) * recip(x1)::R)
+
+           else
+             -- Emulate sqrt using exp and log
+             argument x ==
+               n1 := exp(half*log(norm(x)))
+               x1 := real(x) + n1
+               (2::R)*atan(imag(x) * recip(x1)::R)
 
          pi()  == pi()$R :: %
 


-- 
                              Waldek Hebisch
address@hidden