double pow_di(ap, bp) double *ap; long int *bp; { double pow, x; long int n; pow = 1; x = *ap; n = *bp; if(n != 0) { if(n < 0) { if(x == 0) { return(pow); } n = -n; x = 1/x; } for( ; ; ) { if(n & 01) pow *= x; if(n >>= 1) x *= x; else break; } } return(pow); }