xref: /csrg-svn/usr.bin/f77/libF77/pow_di.c (revision 10512)
1*10512Sdlw /*
2*10512Sdlw  *	"@(#)pow_di.c	1.1"
3*10512Sdlw  */
4*10512Sdlw 
5*10512Sdlw double pow_di(ap, bp)
6*10512Sdlw double *ap;
7*10512Sdlw long int *bp;
8*10512Sdlw {
9*10512Sdlw double pow, x;
10*10512Sdlw long int n;
11*10512Sdlw 
12*10512Sdlw pow = 1;
13*10512Sdlw x = *ap;
14*10512Sdlw n = *bp;
15*10512Sdlw 
16*10512Sdlw if(n != 0)
17*10512Sdlw 	{
18*10512Sdlw 	if(n < 0)
19*10512Sdlw 		{
20*10512Sdlw 		if(x == 0)
21*10512Sdlw 			{
22*10512Sdlw 			return(pow);
23*10512Sdlw 			}
24*10512Sdlw 		n = -n;
25*10512Sdlw 		x = 1/x;
26*10512Sdlw 		}
27*10512Sdlw 	for( ; ; )
28*10512Sdlw 		{
29*10512Sdlw 		if(n & 01)
30*10512Sdlw 			pow *= x;
31*10512Sdlw 		if(n >>= 1)
32*10512Sdlw 			x *= x;
33*10512Sdlw 		else
34*10512Sdlw 			break;
35*10512Sdlw 		}
36*10512Sdlw 	}
37*10512Sdlw return(pow);
38*10512Sdlw }
39