| /* Support routines for the intrinsic power (**) operator. |

| Copyright 2004, 2007, 2009 Free Software Foundation, Inc. |

| Contributed by Paul Brook |

| |

| This file is part of the GNU Fortran 95 runtime library (libgfortran). |

| |

| Libgfortran is free software; you can redistribute it and/or |

| modify it under the terms of the GNU General Public |

| License as published by the Free Software Foundation; either |

| version 3 of the License, or (at your option) any later version. |

| |

| Libgfortran is distributed in the hope that it will be useful, |

| but WITHOUT ANY WARRANTY; without even the implied warranty of |

| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |

| GNU General Public License for more details. |

| |

| Under Section 7 of GPL version 3, you are granted additional |

| permissions described in the GCC Runtime Library Exception, version |

| 3.1, as published by the Free Software Foundation. |

| |

| You should have received a copy of the GNU General Public License and |

| a copy of the GCC Runtime Library Exception along with this program; |

| see the files COPYING3 and COPYING.RUNTIME respectively. If not, see |

| <http://www.gnu.org/licenses/>. */ |

| |

| #include "libgfortran.h" |

| |

| |

| /* Use Binary Method to calculate the powi. This is not an optimal but |

| a simple and reasonable arithmetic. See section 4.6.3, "Evaluation of |

| Powers" of Donald E. Knuth, "Seminumerical Algorithms", Vol. 2, "The Art |

| of Computer Programming", 3rd Edition, 1998. */ |

| |

| #if defined (HAVE_GFC_COMPLEX_4) && defined (HAVE_GFC_INTEGER_16) |

| |

| GFC_COMPLEX_4 pow_c4_i16 (GFC_COMPLEX_4 a, GFC_INTEGER_16 b); |

| export_proto(pow_c4_i16); |

| |

| GFC_COMPLEX_4 |

| pow_c4_i16 (GFC_COMPLEX_4 a, GFC_INTEGER_16 b) |

| { |

| GFC_COMPLEX_4 pow, x; |

| GFC_INTEGER_16 n; |

| GFC_UINTEGER_16 u; |

| |

| n = b; |

| x = a; |

| pow = 1; |

| if (n != 0) |

| { |

| if (n < 0) |

| { |

| |

| u = -n; |

| x = pow / x; |

| } |

| else |

| { |

| u = n; |

| } |

| for (;;) |

| { |

| if (u & 1) |

| pow *= x; |

| u >>= 1; |

| if (u) |

| x *= x; |

| else |

| break; |

| } |

| } |

| return pow; |

| } |

| |

| #endif |