/****************      Euler's Formula       ************************
Author:  Steve Meardon  meardon@dis1.cide.mx
                        smeardon@econ.duke.edu
Date: 8 Feb 98
This code is provided gratis, without any warrantees or guarantees, for
public use.

The problem: in a system of non-linear equations I was attempting
to solve using NLSYS, some of the endogenous variables were taken to the
power of a parameter that I had set to a non-integer value.  The only
solutions of interest to me -- perhaps the only solutions that existed --
were positive solutions for those endogenous variables.  Nevertheless, the
algorithm used by NLSYS tries out both positive and negative values on its
way towards a solution -- and when it tries taking a negative value to a
non-integer power, the program terminates with an error message.

Solving this problem turns out to require a bit more than 'telling' Gauss
(by way of complex(.,.)) that it will be dealing with a complex number if a
negative number is raised to a non-integer power.  Although that will make
Gauss in general accept the operation, NLSYS still doesn't want to deal with
complex numbers.

I figured if I could get NLSYS to find a solution to the system:
y = -x - 3
y = x^(2.8)

then I could get it to deal adequately enough with negative numbers to
non-integer powers.  I can graph this system in Mathcad; since the second
equation is an 'almost cube' equation, then it ought to be almost symmetric
in the first and third quadrants.  Thus, if I'm looking for the solution
that Mathcad is plotting graphically, I should hope to find something close
to the negatives of the solutions to y=x+3, y=x^(2.8).  That is, something
around (-1.2, -1.7).

It took me a while to understand why Mathcad can plot the solution in the
space of REAL numbers, yet Gauss can't find a solution at all.  When
plugging negative x's into the almost cube equation, Mathcad was
transforming the resulting complex numbers into reals.  That, then, is
exactly what one wants to do in order to get NLSYS to work in this case.
This transformation is new to me; those who are more familiar will have
other sources, but I found the formula (Euler's Formula) on p.235 of:
Ross, Shepley. _Differential Equations_. Waltham, MA: Blaisdell Publishing,
1964.

The following program is designed along the lines suggested by Michael
Loretan.  The key is to define a procedure (here,pow(base,expo)) that
first tells Gauss to understand it is dealing with a
complex number if the base is negative and the exponent is non-integer;
second, that converts the complex number to a real.
This program does so, and returns the expected solution.
***********************************************************************/
  library nlsys;
  #include nlsys.ext;
  nlset;

  proc pow(base,expo);
    local s1,is1,rs1,s2,rs2;

    if ( (base<0) and (round(expo) ne expo) );
      s1=(complex(base,0)^expo);
      is1=(s1-real(s1));
      rs1=real(s1);
      s2=rs1+ln(cos(is1/complex(0,1))+complex(0,1)*sin(is1/complex(0,1)));
      rs2=real(s2);
      retp(rs2);

    else;
      retp(base^expo);
    endif;
  endp;

  proc f(v);
    local eqns,x,y;
    x=v[1];y=v[2];

    eqns = y - pow(x, 2.8)|
           y - (-x - 3);

    retp(eqns);
  endp;

  x0={-3,-3};
  _nlalgr=1;
  _nlchpf=0;

  {v,fvp,jc,tcode}=nlprt(nlsys(&f,x0));
  print v;