This is the mail archive of the
gcc-patches@gcc.gnu.org
mailing list for the GCC project.
[fortran, patch] Expand powi to optimal multiplications
- From: Feng Wang <wf_cs at yahoo dot com>
- To: fortran <fortran at gcc dot gnu dot org>, patch <gcc-patches at gcc dot gnu dot org>
- Date: Fri, 30 Apr 2004 17:33:33 +0800 (CST)
- Subject: [fortran, patch] Expand powi to optimal multiplications
Hi, Tobi, Paul and others
This is my implementation to expand powi operator to optimal multiplications.
Since there is no echo in gcc community, we had better implement it in the FE
first. Tested on i686 and ia64.
Consider v ** n. If n is of real type, I use builtin_pow or builtin_powf. If n
is of complex type, just call the runtime library functions (implemented). If n
is an integer variable, call the library powi serial functions implemted in the
attached file powi.c. If n is a constant integer, expand the operator using
"addition chain method" and "window method" arithmetic. The expanding is
controlled by unsafe_math_optimizations and optimize_size option flags.
Feng Wang
Fortran ChangeLog:
2004-04-30 Feng Wang <fengwang@nudt.edu.cn>
* f95-lang.c (gfc_init_builtin_functions): Fix a bug of voidchain. Add
power and powerf to the built_in_decls.
* trans-decl.c (gfor_fndecl_math_pow*i): Add pow_i4i, pow_i8i, pow_r4i,
pow_r8i, pow_c4i, pow_c8i runtime library function declarations.
(gfor_fndecl_math_pow[f]): Delete.
(gfc_build_intrinsic_function_decls): Build function declarations for
gfor_fndecl_math_pow*i and delete gfor_fndecl_math_pow[f].
* trans.h (gfor_fndecl_math_pow*i): Extern the new declarations and
delete gfor_fndecl_math_pow[f].
* trans-expr.c (gfc_conv_integer_power): Delete.
(powi_table): Lookup table used for expanding powi to optimal
multiplications.
(gfc_expand_powi): New recursive function to expand power operator.
(gfc_expand_cnst_int_power): New function. Expand powi.
(gfc_conv_power_op): Use the new functions and runtime library function
decls.
libgfortran ChangeLog:
2004-04-30 Feng Wang <fengwang@nudt.edu.cn>
* runtime/powi.c: New file to implement raising a value to an integer
power.
* Makefile.am: Add it.
* Makefile.in: Regenerate.
_________________________________________________________
Do You Yahoo!?
惠普TT游戏剧，玩游戏，中大奖！
http://cn.rd.yahoo.com/mail_cn/tag/SIG=1402c0to2/**http%3A%2F%2Fhp.allyes.com%2Flaserjet%2Fgamestory%2Findex.html%3Fjumpid%3Dex_hphqapcn_MongooseLJ1010%2F201073CN407016%2FYahoodiff -c3p ../fortran/f95-lang.c ../gcc/gcc/fortran/f95-lang.c
*** ../fortran/f95-lang.c 2004-04-30 19:50:59.000000000 -0800
--- ../gcc/gcc/fortran/f95-lang.c 2004-04-30 19:56:44.000000000 -0800
*************** gfc_init_builtin_functions (void)
*** 751,757 ****
tree tmp;
tree voidchain;
! voidchain = tree_cons (NULL_TREE, void_type_node, NULL_TREE);
tmp = tree_cons (NULL_TREE, float_type_node, voidchain);
mfunc_float[0] = build_function_type (float_type_node, tmp);
--- 751,757 ----
tree tmp;
tree voidchain;
! voidchain = void_list_node;
tmp = tree_cons (NULL_TREE, float_type_node, voidchain);
mfunc_float[0] = build_function_type (float_type_node, tmp);
*************** gfc_init_builtin_functions (void)
*** 765,770 ****
--- 765,772 ----
#include "mathbuiltins.def"
+ DEFINE_MATH_BUILTIN (POW, "pow", 2)
+
/* We define there seperately as the fortran versions have different
semantics (they return an integer type) */
gfc_define_builtin ("__builtin_floor", mfunc_double[0],
diff -c3p ../fortran/trans-decl.c ../gcc/gcc/fortran/trans-decl.c
*** ../fortran/trans-decl.c 2004-04-30 19:50:59.000000000 -0800
--- ../gcc/gcc/fortran/trans-decl.c 2004-04-30 20:08:40.000000000 -0800
*************** tree gfor_fndecl_associated;
*** 93,100 ****
/* Math functions. Many other math functions are handled in
trans-intrinsic.c. */
! tree gfor_fndecl_math_powf;
! tree gfor_fndecl_math_pow;
tree gfor_fndecl_math_cpowf;
tree gfor_fndecl_math_cpow;
tree gfor_fndecl_math_cabsf;
--- 93,104 ----
/* Math functions. Many other math functions are handled in
trans-intrinsic.c. */
! tree gfor_fndecl_math_powi4i;
! tree gfor_fndecl_math_powi8i;
! tree gfor_fndecl_math_powr4i;
! tree gfor_fndecl_math_powr8i;
! tree gfor_fndecl_math_powc4i;
! tree gfor_fndecl_math_powc8i;
tree gfor_fndecl_math_cpowf;
tree gfor_fndecl_math_cpow;
tree gfor_fndecl_math_cabsf;
*************** gfc_build_intrinsic_function_decls (void
*** 1398,1411 ****
/* Power functions. */
! gfor_fndecl_math_powf =
! gfc_build_library_function_decl (get_identifier ("powf"),
gfc_real4_type_node,
! 1, gfc_real4_type_node);
! gfor_fndecl_math_pow =
! gfc_build_library_function_decl (get_identifier ("pow"),
gfc_real8_type_node,
! 1, gfc_real8_type_node);
gfor_fndecl_math_cpowf =
gfc_build_library_function_decl (get_identifier ("cpowf"),
gfc_complex4_type_node,
--- 1402,1441 ----
/* Power functions. */
! gfor_fndecl_math_powi4i =
! gfc_build_library_function_decl (get_identifier (PREFIX("pow_i4i")),
! gfc_int4_type_node,
! 2, gfc_int4_type_node,
! gfc_int4_type_node);
! gfor_fndecl_math_powi8i =
! gfc_build_library_function_decl (get_identifier (PREFIX("pow_i8i")),
! gfc_int8_type_node,
! 2, gfc_int8_type_node,
! gfc_int4_type_node);
! gfor_fndecl_math_powr4i =
! gfc_build_library_function_decl (get_identifier (PREFIX("pow_r4i")),
gfc_real4_type_node,
! 2, gfc_real4_type_node,
! gfc_int4_type_node);
!
! gfor_fndecl_math_powr8i =
! gfc_build_library_function_decl (get_identifier (PREFIX("pow_r8i")),
gfc_real8_type_node,
! 2, gfc_real8_type_node,
! gfc_int4_type_node);
!
! gfor_fndecl_math_powc4i =
! gfc_build_library_function_decl (get_identifier (PREFIX("pow_c4i")),
! gfc_complex4_type_node,
! 2, gfc_complex4_type_node,
! gfc_int4_type_node);
!
! gfor_fndecl_math_powc8i =
! gfc_build_library_function_decl (get_identifier (PREFIX("pow_c8i")),
! gfc_complex8_type_node,
! 2, gfc_complex8_type_node,
! gfc_int4_type_node);
!
gfor_fndecl_math_cpowf =
gfc_build_library_function_decl (get_identifier ("cpowf"),
gfc_complex4_type_node,
diff -c3p ../fortran/trans-expr.c ../gcc/gcc/fortran/trans-expr.c
*** ../fortran/trans-expr.c 2004-04-30 19:50:59.000000000 -0800
--- ../gcc/gcc/fortran/trans-expr.c 2004-04-30 20:04:32.000000000 -0800
*************** gfc_conv_unary_op (enum tree_code code,
*** 382,560 ****
}
! /* For power op (lhs ** rhs) We generate:
! m = lhs
! if (rhs > 0)
! count = rhs
! else if (rhs == 0)
! {
! count = 0
! m = 1
! }
! else // (rhs < 0)
! {
! count = -rhs
! m = 1 / m;
! }
! // for constant rhs we do the above at compile time
! val = m;
! for (n = 1; n < count; n++)
! val = val * m;
! */
! static void
! gfc_conv_integer_power (gfc_se * se, tree lhs, tree rhs)
{
! tree count;
! tree result;
! tree cond;
! tree neg_stmt;
! tree pos_stmt;
tree tmp;
! tree var;
! tree type;
! stmtblock_t block;
! tree exit_label;
! type = TREE_TYPE (lhs);
!
! if (INTEGER_CST_P (rhs))
{
! if (integer_zerop (rhs))
! {
! se->expr = gfc_build_const (type, integer_one_node);
! return;
! }
! /* Special cases for constant values. */
! if (TREE_INT_CST_HIGH (rhs) == -1)
! {
! /* x ** (-y) == 1 / (x ** y). */
! if (TREE_CODE (type) == INTEGER_TYPE)
! {
! se->expr = integer_zero_node;
! return;
! }
!
! tmp = gfc_build_const (type, integer_one_node);
! lhs = fold (build (RDIV_EXPR, type, tmp, lhs));
!
! rhs = fold (build1 (NEGATE_EXPR, TREE_TYPE (rhs), rhs));
! assert (INTEGER_CST_P (rhs));
! }
! else
! {
! /* TODO: really big integer powers. */
! assert (TREE_INT_CST_HIGH (rhs) == 0);
! }
!
! if (integer_onep (rhs))
! {
! se->expr = lhs;
! return;
! }
! if (TREE_INT_CST_LOW (rhs) == 2)
! {
! se->expr = build (MULT_EXPR, type, lhs, lhs);
! return;
! }
! if (TREE_INT_CST_LOW (rhs) == 3)
! {
! tmp = build (MULT_EXPR, type, lhs, lhs);
! se->expr = fold (build (MULT_EXPR, type, tmp, lhs));
! return;
! }
! /* Create the loop count variable. */
! count = gfc_create_var (TREE_TYPE (rhs), "count");
! gfc_add_modify_expr (&se->pre, count, rhs);
}
else
{
! /* Put the lhs into a temporary variable. */
! var = gfc_create_var (type, "val");
! count = gfc_create_var (TREE_TYPE (rhs), "count");
! gfc_add_modify_expr (&se->pre, var, lhs);
! lhs = var;
!
! /* Generate code for negative rhs. */
! gfc_start_block (&block);
!
! if (TREE_CODE (TREE_TYPE (lhs)) == INTEGER_TYPE)
! {
! gfc_add_modify_expr (&block, lhs, integer_zero_node);
! gfc_add_modify_expr (&block, count, integer_zero_node);
! }
! else
! {
! tmp = gfc_build_const (type, integer_one_node);
! tmp = build (RDIV_EXPR, type, tmp, lhs);
! gfc_add_modify_expr (&block, var, tmp);
! tmp = build1 (NEGATE_EXPR, TREE_TYPE (rhs), rhs);
! gfc_add_modify_expr (&block, count, tmp);
! }
! neg_stmt = gfc_finish_block (&block);
! pos_stmt = build_v (MODIFY_EXPR, count, rhs);
! /* Code for rhs == 0. */
! gfc_start_block (&block);
! gfc_add_modify_expr (&block, count, integer_zero_node);
! tmp = gfc_build_const (type, integer_one_node);
! gfc_add_modify_expr (&block, lhs, tmp);
! tmp = gfc_finish_block (&block);
! /* Select the appropriate action. */
! cond = build (EQ_EXPR, boolean_type_node, rhs, integer_zero_node);
! tmp = build_v (COND_EXPR, cond, tmp, neg_stmt);
! cond = build (GT_EXPR, boolean_type_node, rhs, integer_zero_node);
! tmp = build_v (COND_EXPR, cond, pos_stmt, tmp);
! gfc_add_expr_to_block (&se->pre, tmp);
}
! /* Create a variable for the result. */
! result = gfc_create_var (type, "pow");
! gfc_add_modify_expr (&se->pre, result, lhs);
!
! exit_label = gfc_build_label_decl (NULL_TREE);
! TREE_USED (exit_label) = 1;
!
! /* Create the loop body. */
! gfc_start_block (&block);
!
! /* First the exit condition (until count <= 1). */
! tmp = build1_v (GOTO_EXPR, exit_label);
! cond = build (LE_EXPR, TREE_TYPE (count), count, integer_one_node);
! tmp = build_v (COND_EXPR, cond, tmp, build_empty_stmt ());
! gfc_add_expr_to_block (&block, tmp);
!
! /* Multiply by the lhs. */
! tmp = build (MULT_EXPR, type, result, lhs);
! gfc_add_modify_expr (&block, result, tmp);
!
! /* Adjust the loop count. */
! tmp = build (MINUS_EXPR, TREE_TYPE (count), count, integer_one_node);
! gfc_add_modify_expr (&block, count, tmp);
!
! tmp = gfc_finish_block (&block);
!
! /* Create the the loop. */
! tmp = build_v (LOOP_EXPR, tmp);
! gfc_add_expr_to_block (&se->pre, tmp);
! /* Add the exit label. */
! tmp = build1_v (LABEL_EXPR, exit_label);
! gfc_add_expr_to_block (&se->pre, tmp);
! se->expr = result;
}
! /* Power op (**). Integer rhs has special handling. */
static void
gfc_conv_power_op (gfc_se * se, gfc_expr * expr)
--- 382,542 ----
}
+ /* Expand power operator to optimal multiplications when a value is raised
+ to an constant integer n. See section 4.6.3, "Evaluation of Powers" of
+ Donald E. Knuth, "Seminumerical Algorithms", Vol. 2, "The Art of Computer
+ Programming", 3rd Edition, 1998. */
+
+ /* We establish the "optimal power tree" lookup table with the defined size.
+ The items in the table are the exponents used to calculate the index
+ exponents. Any integer n less than the value can get an "addition chain",
+ with the first node being one. */
+ #define POWI_TABLE_SIZE 256
+
+ /* The table is from Builtins.c. */
+ static const unsigned char powi_table[POWI_TABLE_SIZE] =
+ {
+ 0, 1, 1, 2, 2, 3, 3, 4, /* 0 - 7 */
+ 4, 6, 5, 6, 6, 10, 7, 9, /* 8 - 15 */
+ 8, 16, 9, 16, 10, 12, 11, 13, /* 16 - 23 */
+ 12, 17, 13, 18, 14, 24, 15, 26, /* 24 - 31 */
+ 16, 17, 17, 19, 18, 33, 19, 26, /* 32 - 39 */
+ 20, 25, 21, 40, 22, 27, 23, 44, /* 40 - 47 */
+ 24, 32, 25, 34, 26, 29, 27, 44, /* 48 - 55 */
+ 28, 31, 29, 34, 30, 60, 31, 36, /* 56 - 63 */
+ 32, 64, 33, 34, 34, 46, 35, 37, /* 64 - 71 */
+ 36, 65, 37, 50, 38, 48, 39, 69, /* 72 - 79 */
+ 40, 49, 41, 43, 42, 51, 43, 58, /* 80 - 87 */
+ 44, 64, 45, 47, 46, 59, 47, 76, /* 88 - 95 */
+ 48, 65, 49, 66, 50, 67, 51, 66, /* 96 - 103 */
+ 52, 70, 53, 74, 54, 104, 55, 74, /* 104 - 111 */
+ 56, 64, 57, 69, 58, 78, 59, 68, /* 112 - 119 */
+ 60, 61, 61, 80, 62, 75, 63, 68, /* 120 - 127 */
+ 64, 65, 65, 128, 66, 129, 67, 90, /* 128 - 135 */
+ 68, 73, 69, 131, 70, 94, 71, 88, /* 136 - 143 */
+ 72, 128, 73, 98, 74, 132, 75, 121, /* 144 - 151 */
+ 76, 102, 77, 124, 78, 132, 79, 106, /* 152 - 159 */
+ 80, 97, 81, 160, 82, 99, 83, 134, /* 160 - 167 */
+ 84, 86, 85, 95, 86, 160, 87, 100, /* 168 - 175 */
+ 88, 113, 89, 98, 90, 107, 91, 122, /* 176 - 183 */
+ 92, 111, 93, 102, 94, 126, 95, 150, /* 184 - 191 */
+ 96, 128, 97, 130, 98, 133, 99, 195, /* 192 - 199 */
+ 100, 128, 101, 123, 102, 164, 103, 138, /* 200 - 207 */
+ 104, 145, 105, 146, 106, 109, 107, 149, /* 208 - 215 */
+ 108, 200, 109, 146, 110, 170, 111, 157, /* 216 - 223 */
+ 112, 128, 113, 130, 114, 182, 115, 132, /* 224 - 231 */
+ 116, 200, 117, 132, 118, 158, 119, 206, /* 232 - 239 */
+ 120, 240, 121, 162, 122, 147, 123, 152, /* 240 - 247 */
+ 124, 166, 125, 214, 126, 138, 127, 153, /* 248 - 255 */
+ };
! /* If n is larger than lookup table's max index, we use "window method". */
! #define POWI_WINDOW_SIZE 3
! /* Recursive function to expand power operator. The temporary values are put
! in tmpvar. The function return tmpvar[1] ** n. */
! static tree
! gfc_expand_powi (gfc_se * se, int n, tree * tmpvar)
{
! tree op0;
! tree op1;
tree tmp;
! int digit;
! if (n < POWI_TABLE_SIZE)
{
! if (tmpvar[n])
! return tmpvar[n];
! op0 = gfc_expand_powi (se, n - powi_table[n], tmpvar);
! op1 = gfc_expand_powi (se, powi_table[n], tmpvar);
! }
! else if (n & 1)
! {
! digit = n & ((1 << POWI_WINDOW_SIZE) - 1);
! op0 = gfc_expand_powi (se, n - digit, tmpvar);
! op1 = gfc_expand_powi (se, digit, tmpvar);
}
else
{
! op0 = gfc_expand_powi (se, n >> 1, tmpvar);
! op1 = op0;
! }
! tmp = fold (build (MULT_EXPR, TREE_TYPE (op0), op0, op1));
! tmp = gfc_evaluate_now (tmp, &se->pre);
! if (n < POWI_TABLE_SIZE)
! tmpvar[n] = tmp;
! return tmp;
! }
! /* Expand lhs ** rhs. rhs is an constant integer. If expand successfully,
! return 1. Else return 0 and will call runtime library functions. */
! static int
! gfc_expand_cnst_int_power (gfc_se * se, tree lhs, tree rhs)
! {
! tree cond;
! tree tmp;
! tree type;
! tree vartmp[POWI_TABLE_SIZE];
! int n;
! int sgn;
! type = TREE_TYPE (lhs);
! n = abs (TREE_INT_CST_LOW (rhs));
! sgn = tree_int_cst_sgn (rhs);
! if ((!flag_unsafe_math_optimizations || optimize_size) && (n > 2 || n < -1))
! return 0;
! /* rhs == 0 */
! if (sgn == 0)
! {
! se->expr = gfc_build_const (type, integer_one_node);
! return 1;
}
+ /* If rhs < 0 and lhs is an integer, the result is -1, 0 or 1. */
+ if ((sgn == -1) && (TREE_CODE (type) == INTEGER_TYPE))
+ {
+ tmp = build (EQ_EXPR, boolean_type_node, lhs,
+ integer_minus_one_node);
+ cond = build (EQ_EXPR, boolean_type_node, lhs,
+ integer_one_node);
! /* If rhs is an even,
! result = (lhs == 1 || lhs == -1) ? 1 : 0. */
! if ((n & 1) == 0)
! {
! tmp = build (TRUTH_OR_EXPR, boolean_type_node, tmp, cond);
! se->expr = build (COND_EXPR, type, tmp, integer_one_node,
! integer_zero_node);
! return 1;
! }
! /* If rhs is an odd,
! result = (lhs == 1) ? 1 : (lhs == -1) ? -1 : 0. */
! tmp = build (COND_EXPR, type, tmp, integer_minus_one_node,
! integer_zero_node);
! se->expr = build (COND_EXPR, type, cond, integer_one_node,
! tmp);
! return 1;
! }
! memset (vartmp, 0, sizeof (vartmp));
! vartmp[1] = lhs;
! se->expr = gfc_expand_powi (se, n, vartmp);
! if (sgn == -1)
! {
! tmp = gfc_build_const (type, integer_one_node);
! se->expr = build (RDIV_EXPR, type, tmp, se->expr);
! }
! return 1;
}
! /* Power op (**). Constant integer exponent has special handling. */
static void
gfc_conv_power_op (gfc_se * se, gfc_expr * expr)
*************** gfc_conv_power_op (gfc_se * se, gfc_expr
*** 564,570 ****
gfc_se rse;
tree fndecl;
tree tmp;
- tree type;
gfc_init_se (&lse, se);
gfc_conv_expr_val (&lse, expr->op1);
--- 546,551 ----
*************** gfc_conv_power_op (gfc_se * se, gfc_expr
*** 574,597 ****
gfc_conv_expr_val (&rse, expr->op2);
gfc_add_block_to_block (&se->pre, &rse.pre);
! type = TREE_TYPE (lse.expr);
kind = expr->op1->ts.kind;
switch (expr->op2->ts.type)
{
case BT_INTEGER:
! /* Integer powers are expanded inline as multiplications. */
! gfc_conv_integer_power (se, lse.expr, rse.expr);
! return;
!
case BT_REAL:
switch (kind)
{
case 4:
! fndecl = gfor_fndecl_math_powf;
break;
case 8:
! fndecl = gfor_fndecl_math_pow;
break;
default:
abort ();
--- 555,626 ----
gfc_conv_expr_val (&rse, expr->op2);
gfc_add_block_to_block (&se->pre, &rse.pre);
! if (expr->op2->ts.type == BT_INTEGER
! && expr->op2->expr_type == EXPR_CONSTANT)
! if (gfc_expand_cnst_int_power (se, lse.expr, rse.expr))
! return;
kind = expr->op1->ts.kind;
switch (expr->op2->ts.type)
{
case BT_INTEGER:
! rse.expr = convert (integer_type_node, rse.expr);
! switch (expr->op1->ts.type)
! {
! case BT_INTEGER:
! switch (kind)
! {
! case 1:
! case 2:
! case 4:
! lse.expr = convert (integer_type_node, lse.expr);
! fndecl = gfor_fndecl_math_powi4i;
! break;
! case 8:
! fndecl = gfor_fndecl_math_powi8i;
! break;
! default:
! abort ();
! }
! break;
! case BT_REAL:
! switch (kind)
! {
! case 4:
! fndecl = gfor_fndecl_math_powr4i;
! break;
! case 8:
! fndecl = gfor_fndecl_math_powr8i;
! break;
! default:
! abort ();
! }
! break;
! case BT_COMPLEX:
! switch (kind)
! {
! case 4:
! fndecl = gfor_fndecl_math_powc4i;
! break;
! case 8:
! fndecl = gfor_fndecl_math_powc8i;
! break;
! default:
! abort ();
! }
! break;
! default:
! abort ();
! }
! break;
case BT_REAL:
switch (kind)
{
case 4:
! fndecl = built_in_decls[BUILT_IN_POWF];
break;
case 8:
! fndecl = built_in_decls[BUILT_IN_POW];
break;
default:
abort ();
*************** gfc_conv_power_op (gfc_se * se, gfc_expr
*** 619,625 ****
tmp = gfc_chainon_list (NULL_TREE, lse.expr);
tmp = gfc_chainon_list (tmp, rse.expr);
! se->expr = gfc_build_function_call (fndecl, tmp);
}
--- 648,654 ----
tmp = gfc_chainon_list (NULL_TREE, lse.expr);
tmp = gfc_chainon_list (tmp, rse.expr);
! se->expr = fold (gfc_build_function_call (fndecl, tmp));
}
diff -c3p ../fortran/trans.h ../gcc/gcc/fortran/trans.h
*** ../fortran/trans.h 2004-04-30 19:50:59.000000000 -0800
--- ../gcc/gcc/fortran/trans.h 2004-04-30 20:10:08.000000000 -0800
*************** extern GTY(()) tree gfor_fndecl_associat
*** 428,435 ****
/* Math functions. Many other math functions are handled in
trans-intrinsic.c. */
! extern GTY(()) tree gfor_fndecl_math_powf;
! extern GTY(()) tree gfor_fndecl_math_pow;
extern GTY(()) tree gfor_fndecl_math_cpowf;
extern GTY(()) tree gfor_fndecl_math_cpow;
extern GTY(()) tree gfor_fndecl_math_cabsf;
--- 428,440 ----
/* Math functions. Many other math functions are handled in
trans-intrinsic.c. */
!
! extern GTY(()) tree gfor_fndecl_math_powi4i;
! extern GTY(()) tree gfor_fndecl_math_powi8i;
! extern GTY(()) tree gfor_fndecl_math_powr4i;
! extern GTY(()) tree gfor_fndecl_math_powr8i;
! extern GTY(()) tree gfor_fndecl_math_powc4i;
! extern GTY(()) tree gfor_fndecl_math_powc8i;
extern GTY(()) tree gfor_fndecl_math_cpowf;
extern GTY(()) tree gfor_fndecl_math_cpow;
extern GTY(()) tree gfor_fndecl_math_cabsf;
*** ../gcc/libgfortran/Makefile.am 2004-04-30 20:31:02.000000000 -0800
--- ../Makefile.am 2004-04-30 20:30:41.000000000 -0800
*************** runtime/pause.c \
*** 63,68 ****
--- 63,69 ----
runtime/stop.c \
runtime/string.c \
runtime/select.c \
+ runtime/powi.c \
gfortypes.h \
libgfortran.h
/* Implementation of raising a value to an integer power.
Copyright (C) 2004 Free Software Foundation, Inc.
Contributed by Feng Wang <fengwang@nudt.edu.cn>
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 General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
Libgfor 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.
You should have received a copy of the GNU General Public License
along with libgfor; see the file COPYING. If not, write to
the Free Software Foundation, 59 Temple Place - Suite 330,
Boston, MA 02111-1307, USA. */
#include "config.h"
#include <string.h>
#include "libgfortran.h"
#define pow_i4i prefix(pow_i4i)
#define pow_i8i prefix(pow_i8i)
#define pow_r4i prefix(pow_r4i)
#define pow_r8i prefix(pow_r8i)
#define pow_c4i prefix(pow_c4i)
#define pow_c8i prefix(pow_c8i)
/* We 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. */
GFC_INTEGER_4
pow_i4i (GFC_INTEGER_4 a, GFC_INTEGER_4 b)
{
GFC_INTEGER_4 pow, x;
GFC_INTEGER_4 n;
unsigned long u;
x = a;
n = b;
if (n <= 0)
{
if (n == 0 || x == 1)
return 1;
if (x != -1)
return x == 0 ? 1 / x : 0;
n = -n;
}
u = n;
for (pow = 1;;)
{
if (u & 01)
pow *= x;
if (u >>= 1)
x *= x;
else
break;
}
return (pow);
}
GFC_INTEGER_8
pow_i8i (GFC_INTEGER_8 a, GFC_INTEGER_4 b)
{
GFC_INTEGER_8 pow, x;
GFC_INTEGER_4 n;
unsigned long u;
x = a;
n = b;
if (n <= 0)
{
if (n == 0 || x == 1)
return 1;
if (x != -1)
return x == 0 ? 1 / x : 0;
n = -n;
}
u = n;
for (pow = 1;;)
{
if (u & 01)
pow *= x;
if (u >>= 1)
x *= x;
else
break;
}
return (pow);
}
GFC_REAL_4
pow_r4i (GFC_REAL_4 a, GFC_INTEGER_4 b)
{
GFC_REAL_4 pow, x;
GFC_INTEGER_4 n;
unsigned long u;
pow = 1;
x = a;
n = b;
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);
}
GFC_REAL_8
pow_r8i (GFC_REAL_8 a, GFC_INTEGER_4 b)
{
GFC_REAL_8 pow, x;
GFC_INTEGER_4 n;
unsigned long u;
pow = 1;
x = a;
n = b;
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);
}
GFC_COMPLEX_4
pow_c4i (GFC_COMPLEX_4 a, GFC_INTEGER_4 b)
{
GFC_COMPLEX_4 pow, x;
GFC_INTEGER_4 n;
unsigned long u;
n = b;
COMPLEX_ASSIGN(pow, 1.0, 0);
if (n == 0)
goto done;
if (n < 0)
{
n = -n;
x = pow / a;
}
else
x = a;
for (u = n;;)
{
if (u & 01)
pow *= x;
if (u >>= 1)
x *= x;
else
break;
}
done:
return (pow);
}
GFC_COMPLEX_8
pow_c8i (GFC_COMPLEX_8 a, GFC_INTEGER_4 b)
{
GFC_COMPLEX_8 pow, x;
GFC_INTEGER_4 n;
unsigned long u;
n = b;
COMPLEX_ASSIGN(pow, 1.0, 0);
if (n == 0)
goto done;
if (n < 0)
{
n = -n;
x = pow / a;
}
else
x = a;
for (u = n;;)
{
if (u & 01)
pow *= x;
if (u >>= 1)
x *= x;
else
break;
}
done:
return (pow);
}