1*47940Sbostic /*-
2*47940Sbostic * Copyright (c) 1980 The Regents of the University of California.
3*47940Sbostic * All rights reserved.
422939Skre *
5*47940Sbostic * %sccs.include.proprietary.c%
610516Sdlw */
710516Sdlw
8*47940Sbostic #ifndef lint
9*47940Sbostic static char sccsid[] = "@(#)pow_zi.c 5.3 (Berkeley) 04/12/91";
10*47940Sbostic #endif /* not lint */
11*47940Sbostic
1210516Sdlw #include "complex"
1310516Sdlw
1433391Sbostic #define Z_MULEQ(A,B) \
1533391Sbostic t = (A).dreal * (B).dreal - (A).dimag * (B).dimag,\
1633391Sbostic (A).dimag = (A).dreal * (B).dimag + (A).dimag * (B).dreal,\
1733391Sbostic (A).dreal = t /* A *= B */
1833391Sbostic
1933391Sbostic void
pow_zi(p,a,b)2010516Sdlw pow_zi(p, a, b) /* p = a**b */
2133391Sbostic dcomplex *p, *a;
2233391Sbostic long int *b;
2310516Sdlw {
2433391Sbostic register long n = *b;
2533391Sbostic double t;
2633391Sbostic dcomplex x;
2710516Sdlw
2833391Sbostic x = *a;
2933391Sbostic p->dreal = (double)1, p->dimag = (double)0;
3033391Sbostic if (!n)
3133391Sbostic return;
3233391Sbostic if (n < 0) {
3333391Sbostic z_div(&x, p, a);
3433391Sbostic n = -n;
3510516Sdlw }
3633391Sbostic while (!(n&1)) {
3733391Sbostic Z_MULEQ(x, x);
3833391Sbostic n >>= 1;
3910516Sdlw }
4033391Sbostic for (*p = x; --n > 0; Z_MULEQ(*p, x))
4133391Sbostic while (!(n&1)) {
4233391Sbostic Z_MULEQ(x, x);
4333391Sbostic n >>= 1;
4410516Sdlw }
4510516Sdlw }
46