#include "f2c.h" double pow_di(ap, bp) doublereal *ap; integer *bp; { double pow, x; integer 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); }