4.3BSD-Reno/src/lib/libF77/pow_zi.c

Compare this file to the similar file:
Show the results in this format:

/*
 * Copyright (c) 1980 Regents of the University of California.
 * All rights reserved.  The Berkeley software License Agreement
 * specifies the terms and conditions for redistribution.
 *
 *	@(#)pow_zi.c	5.2	1/24/88
 */

#include "complex"

#define	Z_MULEQ(A,B)	\
	t = (A).dreal * (B).dreal - (A).dimag * (B).dimag,\
	(A).dimag = (A).dreal * (B).dimag + (A).dimag * (B).dreal,\
	(A).dreal = t	/* A *= B */

void
pow_zi(p, a, b) 	/* p = a**b  */
	dcomplex *p, *a;
	long int *b;
{
	register long n = *b;
	double t;
	dcomplex x;

	x = *a;
	p->dreal = (double)1, p->dimag = (double)0;
	if (!n)
		return;
	if (n < 0) {
		z_div(&x, p, a);
		n = -n;
	}
	while (!(n&1)) {
		Z_MULEQ(x, x);
		n >>= 1;
	}
	for (*p = x; --n > 0; Z_MULEQ(*p, x))
		while (!(n&1)) {
			Z_MULEQ(x, x);
			n >>= 1;
		}
}