🚀 go-pugleaf

In article <1993Aug19.181352.22494@ast.saic.com>, hilljr@jupiter.saic.com 
  (Jay R. Hill) writes:
|> In article <CBws90.K8B@cid.aes.doe.ca>, ojastej@wfh7546.ice.ncr.ca 
  (James Ojaste) writes:
|> |> In article <1993Aug16.194800.301@cs.uah.edu>, vchagant@ 
  (Venkata Chaganti) writes:
|> |> |> I would like to know is there any public domain routines or library
|> |> |> to do Complex arthimetic in C language.
[...]
|> |> but in C++ it shouldn't be too hard - just overload 
[...]
|> I found a book on C++ by Scott Robert Ladd helpful, it has a PC disk with
|> a C++ complex class.
[...]
|> Now I object to using exp and log (complex) to do a square root!!!  I
|> replace this with
[9 lines of code deleted]
|> but I suspect this is inaccurate for some range of arg(c), near 0 and pi.  
|> ...  Do any of you good folk have an improvement for this?
 
Sorry about following up my own post, making my own thread, but I just 
derived a form which is more accurate than the C++ library on the SUN and 
is likely better than most!  For example, I compare sqrt(-1.0 + 0.0000001i):

               real part           imag part            magnitude
C++ lib:  0.000000049421560621 1.011704190886311361 1.011704190886312471
New form: 0.000000050000000000 1.000000000000001110 1.000000000000002220

The new form is much better!  Code follows:------snippity snip------------
// C++  testsqrt.cxx
// CC -O -o testsqrt testsqrt.cxx -lcomplex -lm
#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <stdlib.h>
void main(int argc, char **argv)
{
  complex c,s; double a,x,y;
  if (argc > 1) x=atof(argv[1]); else x=1;
  if (argc > 2) y=atof(argv[2]); else y=0;
  c=complex(x,y);
  s=sqrt(c);   // default library 
  printf("%.18lf %.18lf %.18lf\n",real(s),imag(s),abs(s));
// new form, optimizations of this are machine/compiler dependent.
  a=sqrt(x*x+y*y);
  if(x>0)
    s=sqrt((a+x)*0.5)*complex(1.0,y/(a+x));
  else
    if(y>=0)
      s=sqrt((a-x)*0.5)*complex(y/(a-x),1.0);
    else
      s=-sqrt((a-x)*0.5)*complex(y/(a-x),1.0);
  printf("%.18lf %.18lf %.18lf\n",real(s),imag(s),abs(s));
}
/* test cases to try:
./testsqrt -1.0 0.0
./testsqrt -0.5 0.5
./testsqrt -0.5 -0.5
./testsqrt 0.5 0.5
./testsqrt 0.5 -0.5
./testsqrt 1.0 0.0000001
./testsqrt -1.0 0.0000001
./testsqrt -1.0 -0.0000001
*/

Warmly,
Jay 'Not tonight honey, I'm on the edge of an hyperbolic component' Hill
-- 
{ hilljr@jupiter.saic.com } begin writeln(3*ln(640320)/sqrt(163):17:15) end.
void main(){double sqrt(), y=1/sqrt(2.), a=.5, m=1,z; int n=0;
for(;m*=2,z=sqrt(1-y*y),y=(1-z)/(1+z),a=a*(1+y)*(1+y)-m*y,n<4;n++);
printf("%17.15lf\n",1/a);}