#!/bin/sh
# This is a shar archive.
# The rest of this file is a shell script which will extract:
#
# 2_2.cmp 2_2_sqrt.c 2_2_tst.c 2_2a1.c 2_2a2.c 2_2b.cmp 2_2b1.c 2_2b2.c 2_2c1.c 2_2c2.c 2_2d1.c 2_2d2.c 2_2e1.c 2_2e2.c 2_2f1.c 2_2g1.c 2_2g2.c 2_2h1.c 2_2i1.c 2_2i2.c 2_2j1.c 2_2j2.c 2_2k1.c 2_2k2.c lsqrt.c makefile rsqrt.c
#
# To extract the files from this shell archive file simply
# create a directory for this file, move the archive file
# to it and enter the command
#
#	sh filename
# 
# The files will be extracted automatically.
# Note: Do not use csh.
#
# Archive created: Mon Jul 30 22:55:42 EDT 1990
#
echo x - 2_2.cmp
sed 's/^X//' > 2_2.cmp << '!EOF!'
c=( 1, -3)
sqrt(c)=( 1.44262, -1.03978)	( 1.44262, -1.03978)
c=( 1, -3)	( 1, -3)

c=( -3, 1)
sqrt(c)=( 0.284849, 1.75532)	( 0.284849, 1.75532)
c=( -3, 1)	( -3, 1)

c=( 3, 0)
sqrt(c)=( 1.73205, 0)	( 1.73205, 0)
c=( 3, 0)	( 3, 0)

c=( 0, 3)
sqrt(c)=( 1.22474, 1.22474)	( 1.22474, 1.22474)
c=( 0, 3)	( 4.44089e-16, 3)

c=( -3, 0)
sqrt(c)=( 0, 1.73205)	( 1.06057e-16, 1.73205)
c=( -3, 0)	( -3, 3.67392e-16)

c=( 0, -3)
sqrt(c)=( 1.22474, -1.22474)	( 1.22474, -1.22474)
c=( 0, -3)	( 4.44089e-16, -3)

!EOF!
ls -l 2_2.cmp
echo x - 2_2_sqrt.c
sed 's/^X//' > 2_2_sqrt.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
#include <complex.h>
#include <math.h>
// determine sqrt of x + yi
complex sqrt(complex z)
{
    double x = real(z);
    double y = imag(z);

    // the easy ones: y == 0
    if (y == 0.0)
	if (x < 0.0)
	    return complex(0.0, sqrt(-x));

	else
	    return complex(sqrt(x), 0.0);

    // almost as easy: x == 0
    if (x == 0.0)
	if (y < 0.0)
	    {
	    double x = sqrt(-y / 2);
	    return complex(x, -x);
	    }

	else
	    {
	    double x = sqrt(y / 2);
	    return complex(x, x);
	    }

    // convert to polar and take the root
    //
    //        2    2  1/2
    // r = ( x  + y )
    //
    // theta = O- = arc tan (y / x)
    //
    //  1/2    1/2
    // z    = r   (cos O-/2 + i sin O-/2)
    double root_r = sqrt(sqrt(x * x + y * y));
    double half_t = atan2(y, x) / 2.0;
    return complex(root_r * cos(half_t),
		   root_r * sin(half_t));
}
!EOF!
ls -l 2_2_sqrt.c
echo x - 2_2_tst.c
sed 's/^X//' > 2_2_tst.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
#include <stream.h>
#include <complex.h>

complex lsqrt(complex);

void tst(complex c)
{
    cout << "c=" << c << "\n";
    complex d = sqrt(c);
    complex dl = lsqrt(c);
    cout << "sqrt(c)=" << d << "\t" << dl << "\n";
    complex e = d * d;
    complex el = dl * dl;
    cout << "c=" << e << "\t" << el << "\n\n";
}

main()
{
    tst(complex(1, -3));
    tst(complex(-3, 1));
    tst(complex(3, 0));
    tst(complex(0, 3));
    tst(complex(-3, 0));
    tst(complex(0, -3));
    return 0;
}
!EOF!
ls -l 2_2_tst.c
echo x - 2_2a1.c
sed 's/^X//' > 2_2a1.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
char ch;
!EOF!
ls -l 2_2a1.c
echo x - 2_2a2.c
sed 's/^X//' > 2_2a2.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
extern char ch;
!EOF!
ls -l 2_2a2.c
echo x - 2_2b.cmp
sed 's/^X//' > 2_2b.cmp << '!EOF!'
c=( 1, -3)
sqrt(c)=( 1.44262, -1.03978)	( 1.44262, -1.03978)
c=( 1, -3)	( 1, -3)

c=( -3, 1)
sqrt(c)=( 0.284849, 1.75532)	( 0.284849, 1.75532)
c=( -3, 1)	( -3, 1)

c=( 3, 0)
sqrt(c)=( 1.73205, 0)	( 1.73205, 0)
c=( 3, 0)	( 3, 0)

c=( 0, 3)
sqrt(c)=( 1.22474, 1.22474)	( 1.22474, 1.22474)
c=( 4.44089e-16, 3)	( 4.44089e-16, 3)

c=( -3, 0)
sqrt(c)=( 0, 1.73205)	( 1.06057e-16, 1.73205)
c=( -3, 0)	( -3, 3.67392e-16)

c=( 0, -3)
sqrt(c)=( 1.22474, -1.22474)	( 1.22474, -1.22474)
c=( 4.44089e-16, -3)	( 4.44089e-16, -3)

!EOF!
ls -l 2_2b.cmp
echo x - 2_2b1.c
sed 's/^X//' > 2_2b1.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
int count = 1;
!EOF!
ls -l 2_2b1.c
echo x - 2_2b2.c
sed 's/^X//' > 2_2b2.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
extern int count;
!EOF!
ls -l 2_2b2.c
echo x - 2_2c1.c
sed 's/^X//' > 2_2c1.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
char *name = "Bjarne";
!EOF!
ls -l 2_2c1.c
echo x - 2_2c2.c
sed 's/^X//' > 2_2c2.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
extern char *name;
!EOF!
ls -l 2_2c2.c
echo x - 2_2d1.c
sed 's/^X//' > 2_2d1.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
struct complex { float re, im; };
!EOF!
ls -l 2_2d1.c
echo x - 2_2d2.c
sed 's/^X//' > 2_2d2.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
struct complex;
complex *x;
!EOF!
ls -l 2_2d2.c
echo x - 2_2e1.c
sed 's/^X//' > 2_2e1.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
complex cvar;
!EOF!
ls -l 2_2e1.c
echo x - 2_2e2.c
sed 's/^X//' > 2_2e2.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
extern complex cvar;
!EOF!
ls -l 2_2e2.c
echo x - 2_2f1.c
sed 's/^X//' > 2_2f1.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
extern complex sqrt(complex);
!EOF!
ls -l 2_2f1.c
echo x - 2_2g1.c
sed 's/^X//' > 2_2g1.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
extern int error_number;
!EOF!
ls -l 2_2g1.c
echo x - 2_2g2.c
sed 's/^X//' > 2_2g2.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
int error_number;
!EOF!
ls -l 2_2g2.c
echo x - 2_2h1.c
sed 's/^X//' > 2_2h1.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
typedef complex point;
!EOF!
ls -l 2_2h1.c
echo x - 2_2i1.c
sed 's/^X//' > 2_2i1.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
float real(complex *p) { return p->re; }
!EOF!
ls -l 2_2i1.c
echo x - 2_2i2.c
sed 's/^X//' > 2_2i2.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
extern float real(complex *p);
!EOF!
ls -l 2_2i2.c
echo x - 2_2j1.c
sed 's/^X//' > 2_2j1.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
const double pi = 3.1415926535897932385;
!EOF!
ls -l 2_2j1.c
echo x - 2_2j2.c
sed 's/^X//' > 2_2j2.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
extern const double pi = 3.1415926535897932385;
!EOF!
ls -l 2_2j2.c
echo x - 2_2k1.c
sed 's/^X//' > 2_2k1.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
struct user;
!EOF!
ls -l 2_2k1.c
echo x - 2_2k2.c
sed 's/^X//' > 2_2k2.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
struct user { char *name; int uid, gid; };
!EOF!
ls -l 2_2k2.c
echo x - lsqrt.c
sed 's/^X//' > lsqrt.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
#include <complex.h>
complex lsqrt(complex z)
{
    double x = real(z);
    double y = imag(z);
    // convert to polar and take root
    //       2    2  1/2
    // r = (x  + y )
    //
    // theta = O = arc tan (y / x)
    //
    //  1/2    1/2
    // z    = r   (cos O/2 + i sin O/2)
    double root_r = sqrt(sqrt(x * x +
			      y * y));
    double half_t = atan2(y, x) / 2.0;
    return complex(root_r * cos(half_t),
		   root_r * sin(half_t));
}
!EOF!
ls -l lsqrt.c
echo x - makefile
sed 's/^X//' > makefile << '!EOF!'
CC= CC -I. -I../../CC

all: 2_2_tst 2_2_tstb

2_2_tst: 2_2_sqrt.o 2_2_tst.o lsqrt.o
	$(CC) 2_2_sqrt.o 2_2_tst.o lsqrt.o -o 2_2_tst -lm

2_2_tstb: rsqrt.o 2_2_tst.o lsqrt.o
	$(CC) 2_2_tst.o rsqrt.o lsqrt.o -o 2_2_tstb -lm

test: all 2_2.cmp 2_2b.cmp
	./2_2_tst > 2_2.out
	cmp 2_2.out 2_2.cmp
	./2_2_tstb > 2_2b.out
	cmp 2_2b.out 2_2b.cmp
	echo tests done
!EOF!
ls -l makefile
echo x - rsqrt.c
sed 's/^X//' > rsqrt.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
/*ident	"@(#)cfront:lib/complex/sqrt.c	1.4" */
# include "complex.h"


#define	SQRT_DANGER	1e17
# define PERIL(t)	(t > SQRT_DANGER || (t < 1/SQRT_DANGER && t != 0))

/*
 *
 * 7-25-83, note from Leonie Rose -
 * Stu Feldman says that the peril tests for the following function
 * are "acceptable" for now, but certain things like
 * sqrt(1e10 + 1e-30*i) will cause floating exceptions.
 *
 */

complex
sqrt(complex z)
{
	complex  answer(1,4.5);
	double  r_old,  partial;

/*
	Check for possible overflow, and fixup if necessary.
*/

	double x = abs(z.re);
	double y = abs(z.im);

	if (x > y && PERIL(x)) {
     		z.im /= x;
    		z.re /= x;  /* z.re is replaced by 1 or -1 */
     		partial = sqrt(x);
     	}
	else if PERIL(y) { 
     		z.im /= y;  /* roles of z.re, z.im reversed from previous */
     		z.re /= y;
     		partial = sqrt(y);
     	}
	else partial = 1;

/*
	Main computation:
	Use half angle formulas to compute angular part of the square root.
	The sign of sin_old is the same as that for sin_new, which means that the
	upper half plane gets mapped to the first quadrant, and
	the lower half plane to the fourth quandrant.
*/


	if (r_old = sqrt(z.re*z.re + z.im*z.im)) {
    		double r_new = partial * sqrt(r_old);

    		double cos_old = z.re/r_old;
    		double sin_old = z.im/r_old;
    		double cos_new = sqrt( (1 + cos_old)/2 );
   		double sin_new = (cos_new == 0)? 1 : sin_old/(2*cos_new);

    		answer.re = r_new * cos_new;
    		answer.im = r_new * sin_new;
    	}

	return answer;
}
!EOF!
ls -l rsqrt.c
# The following exit is to ensure that extra garbage 
# after the end of the shar file will be ignored.
exit 0