1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
#include "f2c.h" double pow_di (doublereal * ap, integer * bp) { double pow, x; integer n; unsigned long u; pow = 1; x = *ap; n = *bp; if (n != 0) { if (n < 0) { n = -n; x = 1 / x; } for (u = n;;) { if (u & 01) pow *= x; if (u >>= 1) x *= x; else break; } } return (pow); }