This is the mail archive of the
mailing list for the GCC project.
Re: [fortran, patch] Expand powi to optimal multiplications
I think that file should have libF77's copyright, as I remarked here:
While comparing the code to libF77's I saw this: is there a reason you
say x = pow/x instead of x = 1/x in the case of a negative exponent?
Paul Brook wrote:
RCS file: libgfortran/m4/pow.m4
diff -N libgfortran/m4/pow.m4
*** /dev/null 1 Jan 1970 00:00:00 -0000
--- b/libgfortran/m4/pow.m4 18 May 2004 21:51:12 -0000
*** 0 ****
--- 1,68 ----
+ `/* Support routines for the intrinsic power (**) operator.
+ Copyright 2004 Free Software Foundation, Inc.
+ Contributed by Paul Brook
+ This file is part of the GNU Fortran 95 runtime library (libgfor).
+ Libgfor is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+ Ligbfor 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 Lesser General Public License for more details.
+ You should have received a copy of the GNU Lesser General Public
+ License along with libgfor; see the file COPYING.LIB. If not,
+ write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+ Boston, MA 02111-1307, USA. */
+ #include "config.h"
+ #include "libgfortran.h"'
+ /* Uuse 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. */
+ `prefix(pow_'rtype_code`_'atype_code) (rtype_name a, atype_name b)
+ rtype_name pow, x;
+ atype_name n, u;
+ n = b;
+ x = a;
+ pow = 1;
+ if (n != 0)
+ if (n < 0)
+ if (x == 1)
+ return 1;
+ if (x == -1)
+ return (n & 1) ? -1 : 1;
+ return (x == 0) ? 1 / x : 0;
+ n = -n;
+ x = pow / x;
+ u = n;
+ for (;;)
+ if (u & 1)
+ pow *= x;
+ u >>= 1;
+ if (u)
+ x *= x;
+ return pow;