Add modular exponentiation for UNSIGNED.

gcc/fortran/ChangeLog:

	* arith.cc (arith_power): Handle modular arithmetic for
	BT_UNSIGNED.
	(eval_intrinsic):  Error for unsigned exponentiation with
	-pedantic.
	* expr.cc (gfc_type_convert_binary): Use type of first
	argument for unsigned exponentiation.
	* gfortran.texi: Mention arithmetic exponentiation.
	* resolve.cc (resolve_operator): Allow unsigned exponentiation.
	* trans-decl.cc (gfc_build_intrinsic_function_decls): Build
	declarations for unsigned exponentiation.
	* trans-expr.cc (gfc_conv_cst_uint_power): New function.
	(gfc_conv_power_op): Call it.  Handle unsigned exponentiation.
	* trans.h (gfor_fndecl_unsigned_pow_list):  Add declaration.

libgfortran/ChangeLog:

	* Makefile.am: Add files for unsigned exponentiation.
	* Makefile.in: Regenerate.
	* gfortran.map: Add functions for unsigned exponentiation.
	* generated/pow_m16_m1.c: New file.
	* generated/pow_m16_m16.c: New file.
	* generated/pow_m16_m2.c: New file.
	* generated/pow_m16_m4.c: New file.
	* generated/pow_m16_m8.c: New file.
	* generated/pow_m1_m1.c: New file.
	* generated/pow_m1_m16.c: New file.
	* generated/pow_m1_m2.c: New file.
	* generated/pow_m1_m4.c: New file.
	* generated/pow_m1_m8.c: New file.
	* generated/pow_m2_m1.c: New file.
	* generated/pow_m2_m16.c: New file.
	* generated/pow_m2_m2.c: New file.
	* generated/pow_m2_m4.c: New file.
	* generated/pow_m2_m8.c: New file.
	* generated/pow_m4_m1.c: New file.
	* generated/pow_m4_m16.c: New file.
	* generated/pow_m4_m2.c: New file.
	* generated/pow_m4_m4.c: New file.
	* generated/pow_m4_m8.c: New file.
	* generated/pow_m8_m1.c: New file.
	* generated/pow_m8_m16.c: New file.
	* generated/pow_m8_m2.c: New file.
	* generated/pow_m8_m4.c: New file.
	* generated/pow_m8_m8.c: New file.
	* m4/powu.m4: New file.

gcc/testsuite/ChangeLog:

	* gfortran.dg/unsigned_15.f90: Adjust error messages.
	* gfortran.dg/unsigned_43.f90: New test.
	* gfortran.dg/unsigned_44.f90: New test.
This commit is contained in:
Thomas Koenig 2025-01-27 18:43:44 +01:00
parent 5b46c01c50
commit c2a0ee5886
39 changed files with 22762 additions and 46 deletions

View file

@ -1143,6 +1143,20 @@ arith_power (gfc_expr *op1, gfc_expr *op2, gfc_expr **resultp)
op2->value.complex, GFC_MPC_RND_MODE);
}
break;
case BT_UNSIGNED:
{
int k;
mpz_t x;
gcc_assert (op1->ts.type == BT_UNSIGNED);
k = gfc_validate_kind (BT_UNSIGNED, op1->ts.kind, false);
/* Exponentiation is performed modulo x = 2**n. */
mpz_init (x);
mpz_add_ui (x, gfc_unsigned_kinds[k].huge, 1);
mpz_powm (result->value.integer, op1->value.integer,
op2->value.integer, x);
mpz_clear (x);
}
break;
default:
gfc_internal_error ("arith_power(): unknown type");
}
@ -1827,10 +1841,11 @@ eval_intrinsic (gfc_intrinsic_op op,
gcc_fallthrough ();
/* Numeric binary */
case INTRINSIC_POWER:
if (flag_unsigned && op == INTRINSIC_POWER)
if (pedantic && (op1->ts.type == BT_UNSIGNED || op2->ts.type == BT_UNSIGNED))
{
if (op1->ts.type == BT_UNSIGNED || op2->ts.type == BT_UNSIGNED)
goto runtime;
gfc_error ("Unsigned exponentiation not permitted with -pedantic "
"at %L", &op1->where);
goto runtime;
}
gcc_fallthrough ();
@ -1940,7 +1955,6 @@ runtime:
/* Create a run-time expression. */
result = gfc_get_operator_expr (&op1->where, op, op1, op2);
result->ts = temp.ts;
return result;
}

View file

@ -981,6 +981,14 @@ gfc_type_convert_binary (gfc_expr *e, int wconversion)
goto done;
}
/* Unsigned exponentiation is special, we need the type of the first
argument here because of modulo arithmetic. */
if (op1->ts.type == BT_UNSIGNED && e->value.op.op == INTRINSIC_POWER)
{
e->ts = op1->ts;
goto done;
}
if (op1->ts.kind > op2->ts.kind)
gfc_convert_type_warn (op2, &op1->ts, 2, wconversion);
else

View file

@ -2772,9 +2772,10 @@ end program main
@noindent
which outputs the number 18446744073709551583.
Arithmetic operations work on unsigned integers, except for exponentiation,
which is prohibited. Unary minus is not permitted when @code{-pedantic}
is in force; this prohibition is part of J3/24-116.txt.
Arithmetic operations work on unsigned integers, also for
exponentiation. As an extension to J3/24-116.txt, unary minus
and exponentiation of unsigned integers are permitted unless
@code{-pedantic} is in force.
In intrinsic procedures, unsigned arguments are typically permitted
for arguments for the data to be processed, analogous to the

View file

@ -4310,19 +4310,6 @@ resolve_operator (gfc_expr *e)
return false;
case INTRINSIC_POWER:
if (flag_unsigned)
{
if (op1->ts.type == BT_UNSIGNED || op2->ts.type == BT_UNSIGNED)
{
CHECK_INTERFACES
gfc_error ("Exponentiation not valid at %L for %s and %s",
&e->where, gfc_typename (op1), gfc_typename (op2));
return false;
}
}
gcc_fallthrough ();
case INTRINSIC_PLUS:
case INTRINSIC_MINUS:
case INTRINSIC_TIMES:

View file

@ -189,6 +189,8 @@ tree gfor_fndecl_caf_random_init;
trans-intrinsic.cc. */
gfc_powdecl_list gfor_fndecl_math_powi[4][3];
tree gfor_fndecl_unsigned_pow_list[5][5];
tree gfor_fndecl_math_ishftc4;
tree gfor_fndecl_math_ishftc8;
tree gfor_fndecl_math_ishftc16;
@ -3746,8 +3748,10 @@ gfc_build_intrinsic_function_decls (void)
int rkind, ikind, jkind;
#define NIKINDS 3
#define NRKINDS 4
static int ikinds[NIKINDS] = {4, 8, 16};
static int rkinds[NRKINDS] = {4, 8, 10, 16};
#define NUKINDS 5
static const int ikinds[NIKINDS] = {4, 8, 16};
static const int rkinds[NRKINDS] = {4, 8, 10, 16};
static const int ukinds[NUKINDS] = {1, 2, 4, 8, 16};
char name[PREFIX_LEN + 12]; /* _gfortran_pow_?n_?n */
for (ikind=0; ikind < NIKINDS; ikind++)
@ -3797,9 +3801,29 @@ gfc_build_intrinsic_function_decls (void)
TREE_NOTHROW (gfor_fndecl_math_powi[rkind][ikind].cmplx) = 1;
}
}
/* For unsigned types, we have every power for every type. */
for (int base = 0; base < NUKINDS; base++)
{
tree base_type = gfc_get_unsigned_type (ukinds[base]);
for (int expon = 0; expon < NUKINDS; expon++)
{
tree expon_type = gfc_get_unsigned_type (ukinds[base]);
if (base_type && expon_type)
{
sprintf (name, PREFIX("pow_m%d_m%d"), ukinds[base],
ukinds[expon]);
gfor_fndecl_unsigned_pow_list [base][expon] =
gfc_build_library_function_decl (get_identifier (name),
base_type, 2, base_type, expon_type);
TREE_READONLY (gfor_fndecl_unsigned_pow_list[base][expon]) = 1;
TREE_NOTHROW (gfor_fndecl_unsigned_pow_list[base][expon]) = 1;
}
}
}
}
#undef NIKINDS
#undef NRKINDS
#undef NUKINDS
}
gfor_fndecl_math_ishftc4 = gfc_build_library_function_decl (

View file

@ -3595,8 +3595,93 @@ gfc_conv_cst_int_power (gfc_se * se, tree lhs, tree rhs)
return 1;
}
/* Convert lhs**rhs, for constant rhs, when both are unsigned.
Method:
if (rhs == 0) ! Checked here.
return 1;
if (lhs & 1 == 1) ! odd_cnd
{
if (bit_size(rhs) < bit_size(lhs)) ! Checked here.
return lhs ** rhs;
/* Power op (**). Constant integer exponent has special handling. */
mask = (1 < bit_size(a) - 1) / 2;
return lhs ** (n & rhs);
}
if (rhs > bit_size(lhs)) ! Checked here.
return 0;
return lhs ** rhs;
*/
static int
gfc_conv_cst_uint_power (gfc_se * se, tree lhs, tree rhs)
{
tree type = TREE_TYPE (lhs);
tree tmp, is_odd, odd_branch, even_branch;
unsigned HOST_WIDE_INT lhs_prec, rhs_prec;
wi::tree_to_wide_ref wrhs = wi::to_wide (rhs);
unsigned HOST_WIDE_INT n, n_odd;
tree vartmp_odd[POWI_TABLE_SIZE], vartmp_even[POWI_TABLE_SIZE];
/* Anything ** 0 is one. */
if (tree_int_cst_sgn (rhs) == 0)
{
se->expr = build_int_cst (type, 1);
return 1;
}
if (!wi::fits_shwi_p (wrhs))
return 0;
n = wrhs.to_uhwi ();
/* tmp = a & 1; . */
tmp = fold_build2_loc (input_location, BIT_AND_EXPR, type,
lhs, build_int_cst (type, 1));
is_odd = fold_build2_loc (input_location, EQ_EXPR, logical_type_node,
tmp, build_int_cst (type, 1));
lhs_prec = TYPE_PRECISION (type);
rhs_prec = TYPE_PRECISION (TREE_TYPE(rhs));
if (rhs_prec >= lhs_prec)
{
unsigned HOST_WIDE_INT mask;
mask = (((unsigned HOST_WIDE_INT) 1) << (lhs_prec - 1)) - 1;
n_odd = n & mask;
}
else
n_odd = n;
memset (vartmp_odd, 0, sizeof (vartmp_odd));
vartmp_odd[0] = build_int_cst(type, 1);
vartmp_odd[1] = lhs;
odd_branch = gfc_conv_powi (se, n_odd, vartmp_odd);
even_branch = NULL_TREE;
if (n > lhs_prec)
even_branch = build_int_cst (type, 0);
else
{
if (n_odd != n)
{
memset (vartmp_even, 0, sizeof (vartmp_even));
vartmp_even[0] = build_int_cst(type, 1);
vartmp_even[1] = lhs;
even_branch = gfc_conv_powi (se, n, vartmp_even);
}
}
if (even_branch != NULL_TREE)
se->expr = fold_build3_loc (input_location, COND_EXPR, type, is_odd,
odd_branch, even_branch);
else
se->expr = odd_branch;
return 1;
}
/* Power op (**). Constant integer exponent and powers of 2 have special
handling. */
static void
gfc_conv_power_op (gfc_se * se, gfc_expr * expr)
@ -3618,13 +3703,29 @@ gfc_conv_power_op (gfc_se * se, gfc_expr * expr)
gfc_conv_expr_val (&rse, expr->value.op.op2);
gfc_add_block_to_block (&se->pre, &rse.pre);
if (expr->value.op.op2->ts.type == BT_INTEGER
if (expr->value.op.op2->expr_type == EXPR_CONSTANT)
{
if (expr->value.op.op2->ts.type == BT_INTEGER)
{
if (gfc_conv_cst_int_power (se, lse.expr, rse.expr))
return;
}
else if (expr->value.op.op2->ts.type == BT_UNSIGNED)
{
if (gfc_conv_cst_uint_power (se, lse.expr, rse.expr))
return;
}
}
if ((expr->value.op.op2->ts.type == BT_INTEGER
|| expr->value.op.op2->ts.type == BT_UNSIGNED)
&& expr->value.op.op2->expr_type == EXPR_CONSTANT)
if (gfc_conv_cst_int_power (se, lse.expr, rse.expr))
return;
if (INTEGER_CST_P (lse.expr)
&& TREE_CODE (TREE_TYPE (rse.expr)) == INTEGER_TYPE)
&& TREE_CODE (TREE_TYPE (rse.expr)) == INTEGER_TYPE
&& expr->value.op.op2->ts.type == BT_INTEGER)
{
wi::tree_to_wide_ref wlhs = wi::to_wide (lse.expr);
HOST_WIDE_INT v;
@ -3724,6 +3825,49 @@ gfc_conv_power_op (gfc_se * se, gfc_expr * expr)
return;
}
}
/* Handle unsigned separate from signed above, things would be too
complicated otherwise. */
if (INTEGER_CST_P (lse.expr) && expr->value.op.op1->ts.type == BT_UNSIGNED)
{
gfc_expr * op1 = expr->value.op.op1;
tree type;
type = TREE_TYPE (lse.expr);
if (mpz_cmp_ui (op1->value.integer, 1) == 0)
{
/* 1**something is always 1. */
se->expr = build_int_cst (type, 1);
return;
}
/* Simplify 2u**x to a shift, with the value set to zero if it falls
outside the range. */
if (mpz_popcount (op1->value.integer) == 1)
{
tree prec_m1, lim, shift, lshift, cond, tmp;
tree rtype = TREE_TYPE (rse.expr);
int e = mpz_scan1 (op1->value.integer, 0);
shift = fold_build2_loc (input_location, MULT_EXPR,
rtype, build_int_cst (rtype, e),
rse.expr);
lshift = fold_build2_loc (input_location, LSHIFT_EXPR, type,
build_int_cst (type, 1), shift);
prec_m1 = fold_build2_loc (input_location, MINUS_EXPR, rtype,
build_int_cst (rtype, TYPE_PRECISION (type)),
build_int_cst (rtype, 1));
lim = fold_build2_loc (input_location, TRUNC_DIV_EXPR, rtype,
prec_m1, build_int_cst (rtype, e));
cond = fold_build2_loc (input_location, GT_EXPR, logical_type_node,
rse.expr, lim);
tmp = fold_build3_loc (input_location, COND_EXPR, type, cond,
build_int_cst (type, 0), lshift);
se->expr = tmp;
return;
}
}
gfc_int4_type_node = gfc_get_int_type (4);
@ -3856,6 +4000,16 @@ gfc_conv_power_op (gfc_se * se, gfc_expr * expr)
fndecl = gfc_builtin_decl_for_float_kind (BUILT_IN_CPOW, kind);
break;
case BT_UNSIGNED:
{
/* Valid kinds for unsigned are 1, 2, 4, 8, 16. Instead of using a
large switch statement, let's just use __builtin_ctz. */
int base = __builtin_ctz (expr->value.op.op1->ts.kind);
int expon = __builtin_ctz (expr->value.op.op2->ts.kind);
fndecl = gfor_fndecl_unsigned_pow_list[base][expon];
}
break;
default:
gcc_unreachable ();
break;

View file

@ -942,6 +942,8 @@ typedef struct GTY(()) gfc_powdecl_list {
gfc_powdecl_list;
extern GTY(()) gfc_powdecl_list gfor_fndecl_math_powi[4][3];
extern GTY(()) tree gfor_fndecl_unsigned_pow_list[5][5];
extern GTY(()) tree gfor_fndecl_math_ishftc4;
extern GTY(()) tree gfor_fndecl_math_ishftc8;
extern GTY(()) tree gfor_fndecl_math_ishftc16;

View file

@ -6,8 +6,8 @@ program main
unsigned :: u
print *,1 + 2u ! { dg-error "Operands of binary numeric operator" }
print *,2u + 1 ! { dg-error "Operands of binary numeric operator" }
print *,2u ** 1 ! { dg-error "Exponentiation not valid" }
print *,2u ** 1u ! { dg-error "Exponentiation not valid" }
print *,2u ** 1 ! { dg-error "Operands of binary numeric operator" }
print *,2u ** 1u
print *,1u < 2 ! { dg-error "Inconsistent types" }
print *,int(1u) < 2
end program main

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

View file

@ -874,6 +874,33 @@ generated/pow_c10_i16.c \
generated/pow_c16_i16.c \
generated/pow_c17_i16.c
i_powu_c = \
generated/pow_m1_m1.c \
generated/pow_m1_m2.c \
generated/pow_m1_m4.c \
generated/pow_m1_m8.c \
generated/pow_m1_m16.c \
generated/pow_m2_m1.c \
generated/pow_m2_m2.c \
generated/pow_m2_m4.c \
generated/pow_m2_m8.c \
generated/pow_m2_m16.c \
generated/pow_m4_m1.c \
generated/pow_m4_m2.c \
generated/pow_m4_m4.c \
generated/pow_m4_m8.c \
generated/pow_m4_m16.c \
generated/pow_m8_m1.c \
generated/pow_m8_m2.c \
generated/pow_m8_m4.c \
generated/pow_m8_m8.c \
generated/pow_m8_m16.c \
generated/pow_m16_m1.c \
generated/pow_m16_m2.c \
generated/pow_m16_m4.c \
generated/pow_m16_m8.c \
generated/pow_m16_m16.c
i_pack_c = \
generated/pack_i1.c \
generated/pack_i2.c \
@ -956,7 +983,7 @@ gfor_built_src= $(i_all_c) $(i_any_c) $(i_count_c) $(i_maxloc0_c) \
$(i_minloc1s_c) $(i_maxloc2s_c) $(i_minloc2s_c) $(i_maxvals_c) \
$(i_maxval0s_c) $(i_minval0s_c) $(i_maxval1s_c) $(i_minval1s_c) \
$(i_findloc0_c) $(i_findloc0s_c) $(i_findloc1_c) $(i_findloc1s_c) \
$(i_findloc2s_c) $(i_isobinding_c)
$(i_findloc2s_c) $(i_isobinding_c) $(i_powu_c)
# Machine generated specifics
gfor_built_specific_src= \
@ -1376,6 +1403,9 @@ $(in_unpack_c): m4/in_unpack.m4 $(I_M4_DEPS)
$(i_pow_c): m4/pow.m4 $(I_M4_DEPS)
$(M4) -Dfile=$@ -I$(srcdir)/m4 pow.m4 > $@
$(i_powu_c): m4/powu.m4 $(I_M4_DEPS)
$(M4) -Dfile=$@ -I$(srcdir)/m4 powu.m4 > $@
$(i_pack_c): m4/pack.m4 $(I_M4_DEPS)
$(M4) -Dfile=$@ -I$(srcdir)/m4 pack.m4 > $@

View file

@ -540,7 +540,20 @@ am__objects_49 = generated/findloc1_i1.lo generated/findloc1_i2.lo \
am__objects_50 = generated/findloc1_s1.lo generated/findloc1_s4.lo
am__objects_51 = generated/findloc2_s1.lo generated/findloc2_s4.lo
am__objects_52 = runtime/ISO_Fortran_binding.lo
am__objects_53 = $(am__objects_4) $(am__objects_5) $(am__objects_6) \
am__objects_53 = generated/pow_m1_m1.lo generated/pow_m1_m2.lo \
generated/pow_m1_m4.lo generated/pow_m1_m8.lo \
generated/pow_m1_m16.lo generated/pow_m2_m1.lo \
generated/pow_m2_m2.lo generated/pow_m2_m4.lo \
generated/pow_m2_m8.lo generated/pow_m2_m16.lo \
generated/pow_m4_m1.lo generated/pow_m4_m2.lo \
generated/pow_m4_m4.lo generated/pow_m4_m8.lo \
generated/pow_m4_m16.lo generated/pow_m8_m1.lo \
generated/pow_m8_m2.lo generated/pow_m8_m4.lo \
generated/pow_m8_m8.lo generated/pow_m8_m16.lo \
generated/pow_m16_m1.lo generated/pow_m16_m2.lo \
generated/pow_m16_m4.lo generated/pow_m16_m8.lo \
generated/pow_m16_m16.lo
am__objects_54 = $(am__objects_4) $(am__objects_5) $(am__objects_6) \
$(am__objects_7) $(am__objects_8) $(am__objects_9) \
$(am__objects_10) $(am__objects_11) $(am__objects_12) \
$(am__objects_13) $(am__objects_14) $(am__objects_15) \
@ -556,16 +569,16 @@ am__objects_53 = $(am__objects_4) $(am__objects_5) $(am__objects_6) \
$(am__objects_43) $(am__objects_44) $(am__objects_45) \
$(am__objects_46) $(am__objects_47) $(am__objects_48) \
$(am__objects_49) $(am__objects_50) $(am__objects_51) \
$(am__objects_52)
@LIBGFOR_MINIMAL_FALSE@am__objects_54 = io/close.lo io/file_pos.lo \
$(am__objects_52) $(am__objects_53)
@LIBGFOR_MINIMAL_FALSE@am__objects_55 = io/close.lo io/file_pos.lo \
@LIBGFOR_MINIMAL_FALSE@ io/format.lo io/inquire.lo \
@LIBGFOR_MINIMAL_FALSE@ io/intrinsics.lo io/list_read.lo \
@LIBGFOR_MINIMAL_FALSE@ io/lock.lo io/open.lo io/read.lo \
@LIBGFOR_MINIMAL_FALSE@ io/transfer.lo io/transfer128.lo \
@LIBGFOR_MINIMAL_FALSE@ io/unit.lo io/unix.lo io/write.lo \
@LIBGFOR_MINIMAL_FALSE@ io/fbuf.lo io/async.lo
am__objects_55 = io/size_from_kind.lo $(am__objects_54)
@LIBGFOR_MINIMAL_FALSE@am__objects_56 = intrinsics/access.lo \
am__objects_56 = io/size_from_kind.lo $(am__objects_55)
@LIBGFOR_MINIMAL_FALSE@am__objects_57 = intrinsics/access.lo \
@LIBGFOR_MINIMAL_FALSE@ intrinsics/c99_functions.lo \
@LIBGFOR_MINIMAL_FALSE@ intrinsics/chdir.lo intrinsics/chmod.lo \
@LIBGFOR_MINIMAL_FALSE@ intrinsics/clock.lo \
@ -589,8 +602,8 @@ am__objects_55 = io/size_from_kind.lo $(am__objects_54)
@LIBGFOR_MINIMAL_FALSE@ intrinsics/system_clock.lo \
@LIBGFOR_MINIMAL_FALSE@ intrinsics/time.lo intrinsics/umask.lo \
@LIBGFOR_MINIMAL_FALSE@ intrinsics/unlink.lo
@IEEE_SUPPORT_TRUE@am__objects_57 = ieee/ieee_helper.lo
am__objects_58 = intrinsics/associated.lo intrinsics/abort.lo \
@IEEE_SUPPORT_TRUE@am__objects_58 = ieee/ieee_helper.lo
am__objects_59 = intrinsics/associated.lo intrinsics/abort.lo \
intrinsics/args.lo intrinsics/cshift0.lo \
intrinsics/eoshift0.lo intrinsics/eoshift2.lo \
intrinsics/erfc_scaled.lo intrinsics/extends_type_of.lo \
@ -605,12 +618,12 @@ am__objects_58 = intrinsics/associated.lo intrinsics/abort.lo \
intrinsics/selected_real_kind.lo intrinsics/trigd.lo \
intrinsics/unpack_generic.lo runtime/in_pack_generic.lo \
runtime/in_unpack_generic.lo runtime/in_pack_class.lo \
runtime/in_unpack_class.lo $(am__objects_56) $(am__objects_57)
@IEEE_SUPPORT_TRUE@am__objects_59 = ieee/ieee_arithmetic.lo \
runtime/in_unpack_class.lo $(am__objects_57) $(am__objects_58)
@IEEE_SUPPORT_TRUE@am__objects_60 = ieee/ieee_arithmetic.lo \
@IEEE_SUPPORT_TRUE@ ieee/ieee_exceptions.lo \
@IEEE_SUPPORT_TRUE@ ieee/ieee_features.lo
am__objects_60 =
am__objects_61 = generated/_abs_c4.lo generated/_abs_c8.lo \
am__objects_61 =
am__objects_62 = generated/_abs_c4.lo generated/_abs_c8.lo \
generated/_abs_c10.lo generated/_abs_c16.lo \
generated/_abs_c17.lo generated/_abs_i4.lo \
generated/_abs_i8.lo generated/_abs_i16.lo \
@ -679,7 +692,7 @@ am__objects_61 = generated/_abs_c4.lo generated/_abs_c8.lo \
generated/_aint_r17.lo generated/_anint_r4.lo \
generated/_anint_r8.lo generated/_anint_r10.lo \
generated/_anint_r16.lo generated/_anint_r17.lo
am__objects_62 = generated/_sign_i4.lo generated/_sign_i8.lo \
am__objects_63 = generated/_sign_i4.lo generated/_sign_i8.lo \
generated/_sign_i16.lo generated/_sign_r4.lo \
generated/_sign_r8.lo generated/_sign_r10.lo \
generated/_sign_r16.lo generated/_sign_r17.lo \
@ -694,13 +707,13 @@ am__objects_62 = generated/_sign_i4.lo generated/_sign_i8.lo \
generated/_mod_r4.lo generated/_mod_r8.lo \
generated/_mod_r10.lo generated/_mod_r16.lo \
generated/_mod_r17.lo
am__objects_63 = generated/misc_specifics.lo
am__objects_64 = $(am__objects_61) $(am__objects_62) $(am__objects_63) \
am__objects_64 = generated/misc_specifics.lo
am__objects_65 = $(am__objects_62) $(am__objects_63) $(am__objects_64) \
intrinsics/dprod_r8.lo intrinsics/f2c_specifics.lo \
intrinsics/random_init.lo
am_libgfortran_la_OBJECTS = $(am__objects_3) $(am__objects_53) \
$(am__objects_55) $(am__objects_58) $(am__objects_59) \
$(am__objects_60) $(am__objects_64)
am_libgfortran_la_OBJECTS = $(am__objects_3) $(am__objects_54) \
$(am__objects_56) $(am__objects_59) $(am__objects_60) \
$(am__objects_61) $(am__objects_65)
libgfortran_la_OBJECTS = $(am_libgfortran_la_OBJECTS)
AM_V_P = $(am__v_P_@AM_V@)
am__v_P_ = $(am__v_P_@AM_DEFAULT_V@)
@ -1667,6 +1680,33 @@ generated/pow_c10_i16.c \
generated/pow_c16_i16.c \
generated/pow_c17_i16.c
i_powu_c = \
generated/pow_m1_m1.c \
generated/pow_m1_m2.c \
generated/pow_m1_m4.c \
generated/pow_m1_m8.c \
generated/pow_m1_m16.c \
generated/pow_m2_m1.c \
generated/pow_m2_m2.c \
generated/pow_m2_m4.c \
generated/pow_m2_m8.c \
generated/pow_m2_m16.c \
generated/pow_m4_m1.c \
generated/pow_m4_m2.c \
generated/pow_m4_m4.c \
generated/pow_m4_m8.c \
generated/pow_m4_m16.c \
generated/pow_m8_m1.c \
generated/pow_m8_m2.c \
generated/pow_m8_m4.c \
generated/pow_m8_m8.c \
generated/pow_m8_m16.c \
generated/pow_m16_m1.c \
generated/pow_m16_m2.c \
generated/pow_m16_m4.c \
generated/pow_m16_m8.c \
generated/pow_m16_m16.c
i_pack_c = \
generated/pack_i1.c \
generated/pack_i2.c \
@ -1749,7 +1789,7 @@ gfor_built_src = $(i_all_c) $(i_any_c) $(i_count_c) $(i_maxloc0_c) \
$(i_minloc1s_c) $(i_maxloc2s_c) $(i_minloc2s_c) $(i_maxvals_c) \
$(i_maxval0s_c) $(i_minval0s_c) $(i_maxval1s_c) $(i_minval1s_c) \
$(i_findloc0_c) $(i_findloc0s_c) $(i_findloc1_c) $(i_findloc1s_c) \
$(i_findloc2s_c) $(i_isobinding_c)
$(i_findloc2s_c) $(i_isobinding_c) $(i_powu_c)
# Machine generated specifics
@ -3315,6 +3355,56 @@ generated/findloc2_s4.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
runtime/ISO_Fortran_binding.lo: runtime/$(am__dirstamp) \
runtime/$(DEPDIR)/$(am__dirstamp)
generated/pow_m1_m1.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m1_m2.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m1_m4.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m1_m8.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m1_m16.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m2_m1.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m2_m2.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m2_m4.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m2_m8.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m2_m16.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m4_m1.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m4_m2.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m4_m4.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m4_m8.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m4_m16.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m8_m1.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m8_m2.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m8_m4.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m8_m8.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m8_m16.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m16_m1.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m16_m2.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m16_m4.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m16_m8.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
generated/pow_m16_m16.lo: generated/$(am__dirstamp) \
generated/$(DEPDIR)/$(am__dirstamp)
io/$(am__dirstamp):
@$(MKDIR_P) io
@: > io/$(am__dirstamp)
@ -4347,6 +4437,31 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_i8_i16.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_i8_i4.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_i8_i8.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m16_m1.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m16_m16.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m16_m2.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m16_m4.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m16_m8.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m1_m1.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m1_m16.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m1_m2.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m1_m4.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m1_m8.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m2_m1.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m2_m16.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m2_m2.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m2_m4.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m2_m8.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m4_m1.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m4_m16.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m4_m2.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m4_m4.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m4_m8.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m8_m1.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m8_m16.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m8_m2.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m8_m4.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_m8_m8.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_r10_i16.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_r10_i8.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@generated/$(DEPDIR)/pow_r16_i16.Plo@am__quote@
@ -5084,6 +5199,9 @@ include/ISO_Fortran_binding.h: $(srcdir)/ISO_Fortran_binding.h
@MAINTAINER_MODE_TRUE@$(i_pow_c): m4/pow.m4 $(I_M4_DEPS)
@MAINTAINER_MODE_TRUE@ $(M4) -Dfile=$@ -I$(srcdir)/m4 pow.m4 > $@
@MAINTAINER_MODE_TRUE@$(i_powu_c): m4/powu.m4 $(I_M4_DEPS)
@MAINTAINER_MODE_TRUE@ $(M4) -Dfile=$@ -I$(srcdir)/m4 powu.m4 > $@
@MAINTAINER_MODE_TRUE@$(i_pack_c): m4/pack.m4 $(I_M4_DEPS)
@MAINTAINER_MODE_TRUE@ $(M4) -Dfile=$@ -I$(srcdir)/m4 pack.m4 > $@

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_16) && defined (HAVE_GFC_UINTEGER_1)
GFC_UINTEGER_16 pow_m16_m1 (GFC_UINTEGER_16 x, GFC_UINTEGER_1 n);
export_proto(pow_m16_m1);
inline static GFC_UINTEGER_16
power_simple_m16_m1 (GFC_UINTEGER_16 x, GFC_UINTEGER_1 n)
{
GFC_UINTEGER_16 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_16
pow_m16_m1 (GFC_UINTEGER_16 x, GFC_UINTEGER_1 n)
{
const GFC_UINTEGER_16 mask = (GFC_UINTEGER_16) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m16_m1 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m16_m1 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_16) && defined (HAVE_GFC_UINTEGER_16)
GFC_UINTEGER_16 pow_m16_m16 (GFC_UINTEGER_16 x, GFC_UINTEGER_16 n);
export_proto(pow_m16_m16);
inline static GFC_UINTEGER_16
power_simple_m16_m16 (GFC_UINTEGER_16 x, GFC_UINTEGER_16 n)
{
GFC_UINTEGER_16 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_16
pow_m16_m16 (GFC_UINTEGER_16 x, GFC_UINTEGER_16 n)
{
const GFC_UINTEGER_16 mask = (GFC_UINTEGER_16) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m16_m16 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m16_m16 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_16) && defined (HAVE_GFC_UINTEGER_2)
GFC_UINTEGER_16 pow_m16_m2 (GFC_UINTEGER_16 x, GFC_UINTEGER_2 n);
export_proto(pow_m16_m2);
inline static GFC_UINTEGER_16
power_simple_m16_m2 (GFC_UINTEGER_16 x, GFC_UINTEGER_2 n)
{
GFC_UINTEGER_16 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_16
pow_m16_m2 (GFC_UINTEGER_16 x, GFC_UINTEGER_2 n)
{
const GFC_UINTEGER_16 mask = (GFC_UINTEGER_16) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m16_m2 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m16_m2 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_16) && defined (HAVE_GFC_UINTEGER_4)
GFC_UINTEGER_16 pow_m16_m4 (GFC_UINTEGER_16 x, GFC_UINTEGER_4 n);
export_proto(pow_m16_m4);
inline static GFC_UINTEGER_16
power_simple_m16_m4 (GFC_UINTEGER_16 x, GFC_UINTEGER_4 n)
{
GFC_UINTEGER_16 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_16
pow_m16_m4 (GFC_UINTEGER_16 x, GFC_UINTEGER_4 n)
{
const GFC_UINTEGER_16 mask = (GFC_UINTEGER_16) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m16_m4 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m16_m4 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_16) && defined (HAVE_GFC_UINTEGER_8)
GFC_UINTEGER_16 pow_m16_m8 (GFC_UINTEGER_16 x, GFC_UINTEGER_8 n);
export_proto(pow_m16_m8);
inline static GFC_UINTEGER_16
power_simple_m16_m8 (GFC_UINTEGER_16 x, GFC_UINTEGER_8 n)
{
GFC_UINTEGER_16 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_16
pow_m16_m8 (GFC_UINTEGER_16 x, GFC_UINTEGER_8 n)
{
const GFC_UINTEGER_16 mask = (GFC_UINTEGER_16) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m16_m8 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m16_m8 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_1) && defined (HAVE_GFC_UINTEGER_1)
GFC_UINTEGER_1 pow_m1_m1 (GFC_UINTEGER_1 x, GFC_UINTEGER_1 n);
export_proto(pow_m1_m1);
inline static GFC_UINTEGER_1
power_simple_m1_m1 (GFC_UINTEGER_1 x, GFC_UINTEGER_1 n)
{
GFC_UINTEGER_1 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_1
pow_m1_m1 (GFC_UINTEGER_1 x, GFC_UINTEGER_1 n)
{
const GFC_UINTEGER_1 mask = (GFC_UINTEGER_1) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m1_m1 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m1_m1 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_1) && defined (HAVE_GFC_UINTEGER_16)
GFC_UINTEGER_1 pow_m1_m16 (GFC_UINTEGER_1 x, GFC_UINTEGER_16 n);
export_proto(pow_m1_m16);
inline static GFC_UINTEGER_1
power_simple_m1_m16 (GFC_UINTEGER_1 x, GFC_UINTEGER_16 n)
{
GFC_UINTEGER_1 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_1
pow_m1_m16 (GFC_UINTEGER_1 x, GFC_UINTEGER_16 n)
{
const GFC_UINTEGER_1 mask = (GFC_UINTEGER_1) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m1_m16 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m1_m16 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_1) && defined (HAVE_GFC_UINTEGER_2)
GFC_UINTEGER_1 pow_m1_m2 (GFC_UINTEGER_1 x, GFC_UINTEGER_2 n);
export_proto(pow_m1_m2);
inline static GFC_UINTEGER_1
power_simple_m1_m2 (GFC_UINTEGER_1 x, GFC_UINTEGER_2 n)
{
GFC_UINTEGER_1 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_1
pow_m1_m2 (GFC_UINTEGER_1 x, GFC_UINTEGER_2 n)
{
const GFC_UINTEGER_1 mask = (GFC_UINTEGER_1) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m1_m2 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m1_m2 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_1) && defined (HAVE_GFC_UINTEGER_4)
GFC_UINTEGER_1 pow_m1_m4 (GFC_UINTEGER_1 x, GFC_UINTEGER_4 n);
export_proto(pow_m1_m4);
inline static GFC_UINTEGER_1
power_simple_m1_m4 (GFC_UINTEGER_1 x, GFC_UINTEGER_4 n)
{
GFC_UINTEGER_1 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_1
pow_m1_m4 (GFC_UINTEGER_1 x, GFC_UINTEGER_4 n)
{
const GFC_UINTEGER_1 mask = (GFC_UINTEGER_1) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m1_m4 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m1_m4 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_1) && defined (HAVE_GFC_UINTEGER_8)
GFC_UINTEGER_1 pow_m1_m8 (GFC_UINTEGER_1 x, GFC_UINTEGER_8 n);
export_proto(pow_m1_m8);
inline static GFC_UINTEGER_1
power_simple_m1_m8 (GFC_UINTEGER_1 x, GFC_UINTEGER_8 n)
{
GFC_UINTEGER_1 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_1
pow_m1_m8 (GFC_UINTEGER_1 x, GFC_UINTEGER_8 n)
{
const GFC_UINTEGER_1 mask = (GFC_UINTEGER_1) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m1_m8 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m1_m8 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_2) && defined (HAVE_GFC_UINTEGER_1)
GFC_UINTEGER_2 pow_m2_m1 (GFC_UINTEGER_2 x, GFC_UINTEGER_1 n);
export_proto(pow_m2_m1);
inline static GFC_UINTEGER_2
power_simple_m2_m1 (GFC_UINTEGER_2 x, GFC_UINTEGER_1 n)
{
GFC_UINTEGER_2 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_2
pow_m2_m1 (GFC_UINTEGER_2 x, GFC_UINTEGER_1 n)
{
const GFC_UINTEGER_2 mask = (GFC_UINTEGER_2) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m2_m1 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m2_m1 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_2) && defined (HAVE_GFC_UINTEGER_16)
GFC_UINTEGER_2 pow_m2_m16 (GFC_UINTEGER_2 x, GFC_UINTEGER_16 n);
export_proto(pow_m2_m16);
inline static GFC_UINTEGER_2
power_simple_m2_m16 (GFC_UINTEGER_2 x, GFC_UINTEGER_16 n)
{
GFC_UINTEGER_2 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_2
pow_m2_m16 (GFC_UINTEGER_2 x, GFC_UINTEGER_16 n)
{
const GFC_UINTEGER_2 mask = (GFC_UINTEGER_2) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m2_m16 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m2_m16 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_2) && defined (HAVE_GFC_UINTEGER_2)
GFC_UINTEGER_2 pow_m2_m2 (GFC_UINTEGER_2 x, GFC_UINTEGER_2 n);
export_proto(pow_m2_m2);
inline static GFC_UINTEGER_2
power_simple_m2_m2 (GFC_UINTEGER_2 x, GFC_UINTEGER_2 n)
{
GFC_UINTEGER_2 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_2
pow_m2_m2 (GFC_UINTEGER_2 x, GFC_UINTEGER_2 n)
{
const GFC_UINTEGER_2 mask = (GFC_UINTEGER_2) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m2_m2 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m2_m2 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_2) && defined (HAVE_GFC_UINTEGER_4)
GFC_UINTEGER_2 pow_m2_m4 (GFC_UINTEGER_2 x, GFC_UINTEGER_4 n);
export_proto(pow_m2_m4);
inline static GFC_UINTEGER_2
power_simple_m2_m4 (GFC_UINTEGER_2 x, GFC_UINTEGER_4 n)
{
GFC_UINTEGER_2 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_2
pow_m2_m4 (GFC_UINTEGER_2 x, GFC_UINTEGER_4 n)
{
const GFC_UINTEGER_2 mask = (GFC_UINTEGER_2) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m2_m4 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m2_m4 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_2) && defined (HAVE_GFC_UINTEGER_8)
GFC_UINTEGER_2 pow_m2_m8 (GFC_UINTEGER_2 x, GFC_UINTEGER_8 n);
export_proto(pow_m2_m8);
inline static GFC_UINTEGER_2
power_simple_m2_m8 (GFC_UINTEGER_2 x, GFC_UINTEGER_8 n)
{
GFC_UINTEGER_2 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_2
pow_m2_m8 (GFC_UINTEGER_2 x, GFC_UINTEGER_8 n)
{
const GFC_UINTEGER_2 mask = (GFC_UINTEGER_2) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m2_m8 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m2_m8 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_4) && defined (HAVE_GFC_UINTEGER_1)
GFC_UINTEGER_4 pow_m4_m1 (GFC_UINTEGER_4 x, GFC_UINTEGER_1 n);
export_proto(pow_m4_m1);
inline static GFC_UINTEGER_4
power_simple_m4_m1 (GFC_UINTEGER_4 x, GFC_UINTEGER_1 n)
{
GFC_UINTEGER_4 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_4
pow_m4_m1 (GFC_UINTEGER_4 x, GFC_UINTEGER_1 n)
{
const GFC_UINTEGER_4 mask = (GFC_UINTEGER_4) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m4_m1 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m4_m1 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_4) && defined (HAVE_GFC_UINTEGER_16)
GFC_UINTEGER_4 pow_m4_m16 (GFC_UINTEGER_4 x, GFC_UINTEGER_16 n);
export_proto(pow_m4_m16);
inline static GFC_UINTEGER_4
power_simple_m4_m16 (GFC_UINTEGER_4 x, GFC_UINTEGER_16 n)
{
GFC_UINTEGER_4 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_4
pow_m4_m16 (GFC_UINTEGER_4 x, GFC_UINTEGER_16 n)
{
const GFC_UINTEGER_4 mask = (GFC_UINTEGER_4) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m4_m16 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m4_m16 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_4) && defined (HAVE_GFC_UINTEGER_2)
GFC_UINTEGER_4 pow_m4_m2 (GFC_UINTEGER_4 x, GFC_UINTEGER_2 n);
export_proto(pow_m4_m2);
inline static GFC_UINTEGER_4
power_simple_m4_m2 (GFC_UINTEGER_4 x, GFC_UINTEGER_2 n)
{
GFC_UINTEGER_4 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_4
pow_m4_m2 (GFC_UINTEGER_4 x, GFC_UINTEGER_2 n)
{
const GFC_UINTEGER_4 mask = (GFC_UINTEGER_4) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m4_m2 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m4_m2 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_4) && defined (HAVE_GFC_UINTEGER_4)
GFC_UINTEGER_4 pow_m4_m4 (GFC_UINTEGER_4 x, GFC_UINTEGER_4 n);
export_proto(pow_m4_m4);
inline static GFC_UINTEGER_4
power_simple_m4_m4 (GFC_UINTEGER_4 x, GFC_UINTEGER_4 n)
{
GFC_UINTEGER_4 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_4
pow_m4_m4 (GFC_UINTEGER_4 x, GFC_UINTEGER_4 n)
{
const GFC_UINTEGER_4 mask = (GFC_UINTEGER_4) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m4_m4 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m4_m4 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_4) && defined (HAVE_GFC_UINTEGER_8)
GFC_UINTEGER_4 pow_m4_m8 (GFC_UINTEGER_4 x, GFC_UINTEGER_8 n);
export_proto(pow_m4_m8);
inline static GFC_UINTEGER_4
power_simple_m4_m8 (GFC_UINTEGER_4 x, GFC_UINTEGER_8 n)
{
GFC_UINTEGER_4 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_4
pow_m4_m8 (GFC_UINTEGER_4 x, GFC_UINTEGER_8 n)
{
const GFC_UINTEGER_4 mask = (GFC_UINTEGER_4) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m4_m8 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m4_m8 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_8) && defined (HAVE_GFC_UINTEGER_1)
GFC_UINTEGER_8 pow_m8_m1 (GFC_UINTEGER_8 x, GFC_UINTEGER_1 n);
export_proto(pow_m8_m1);
inline static GFC_UINTEGER_8
power_simple_m8_m1 (GFC_UINTEGER_8 x, GFC_UINTEGER_1 n)
{
GFC_UINTEGER_8 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_8
pow_m8_m1 (GFC_UINTEGER_8 x, GFC_UINTEGER_1 n)
{
const GFC_UINTEGER_8 mask = (GFC_UINTEGER_8) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m8_m1 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m8_m1 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_8) && defined (HAVE_GFC_UINTEGER_16)
GFC_UINTEGER_8 pow_m8_m16 (GFC_UINTEGER_8 x, GFC_UINTEGER_16 n);
export_proto(pow_m8_m16);
inline static GFC_UINTEGER_8
power_simple_m8_m16 (GFC_UINTEGER_8 x, GFC_UINTEGER_16 n)
{
GFC_UINTEGER_8 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_8
pow_m8_m16 (GFC_UINTEGER_8 x, GFC_UINTEGER_16 n)
{
const GFC_UINTEGER_8 mask = (GFC_UINTEGER_8) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m8_m16 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m8_m16 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_8) && defined (HAVE_GFC_UINTEGER_2)
GFC_UINTEGER_8 pow_m8_m2 (GFC_UINTEGER_8 x, GFC_UINTEGER_2 n);
export_proto(pow_m8_m2);
inline static GFC_UINTEGER_8
power_simple_m8_m2 (GFC_UINTEGER_8 x, GFC_UINTEGER_2 n)
{
GFC_UINTEGER_8 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_8
pow_m8_m2 (GFC_UINTEGER_8 x, GFC_UINTEGER_2 n)
{
const GFC_UINTEGER_8 mask = (GFC_UINTEGER_8) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m8_m2 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m8_m2 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_8) && defined (HAVE_GFC_UINTEGER_4)
GFC_UINTEGER_8 pow_m8_m4 (GFC_UINTEGER_8 x, GFC_UINTEGER_4 n);
export_proto(pow_m8_m4);
inline static GFC_UINTEGER_8
power_simple_m8_m4 (GFC_UINTEGER_8 x, GFC_UINTEGER_4 n)
{
GFC_UINTEGER_8 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_8
pow_m8_m4 (GFC_UINTEGER_8 x, GFC_UINTEGER_4 n)
{
const GFC_UINTEGER_8 mask = (GFC_UINTEGER_8) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m8_m4 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m8_m4 (x, n);
}
#endif

View file

@ -0,0 +1,79 @@
/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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_UINTEGER_8) && defined (HAVE_GFC_UINTEGER_8)
GFC_UINTEGER_8 pow_m8_m8 (GFC_UINTEGER_8 x, GFC_UINTEGER_8 n);
export_proto(pow_m8_m8);
inline static GFC_UINTEGER_8
power_simple_m8_m8 (GFC_UINTEGER_8 x, GFC_UINTEGER_8 n)
{
GFC_UINTEGER_8 pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler's theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
GFC_UINTEGER_8
pow_m8_m8 (GFC_UINTEGER_8 x, GFC_UINTEGER_8 n)
{
const GFC_UINTEGER_8 mask = (GFC_UINTEGER_8) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_m8_m8 (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_m8_m8 (x, n);
}
#endif

View file

@ -1998,4 +1998,29 @@ GFORTRAN_15 {
_gfortran_sminloc1_8_m4;
_gfortran_sminloc1_8_m8;
_gfortran_report_exception;
_gfortran_pow_m1_m1;
_gfortran_pow_m1_m2;
_gfortran_pow_m1_m4;
_gfortran_pow_m1_m8;
_gfortran_pow_m1_m16;
_gfortran_pow_m2_m1;
_gfortran_pow_m2_m2;
_gfortran_pow_m2_m4;
_gfortran_pow_m2_m8;
_gfortran_pow_m2_m16;
_gfortran_pow_m4_m1;
_gfortran_pow_m4_m2;
_gfortran_pow_m4_m4;
_gfortran_pow_m4_m8;
_gfortran_pow_m4_m16;
_gfortran_pow_m8_m1;
_gfortran_pow_m8_m2;
_gfortran_pow_m8_m4;
_gfortran_pow_m8_m8;
_gfortran_pow_m8_m16;
_gfortran_pow_m16_m1;
_gfortran_pow_m16_m2;
_gfortran_pow_m16_m4;
_gfortran_pow_m16_m8;
_gfortran_pow_m16_m16;
} GFORTRAN_14;

80
libgfortran/m4/powu.m4 Normal file
View file

@ -0,0 +1,80 @@
`/* Support routines for the intrinsic power (**) operator
for UNSIGNED, using modulo arithmetic.
Copyright (C) 2025 Free Software Foundation, Inc.
Contributed by Thomas Koenig.
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"'
include(iparm.m4)dnl
/* 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_'rtype_name`) && defined (HAVE_'atype_name`)'
rtype_name `pow_'rtype_code`_'atype_code` ('rtype_name` x, 'atype_name` n);
export_proto(pow_'rtype_code`_'atype_code`);
inline static 'rtype_name`
power_simple_'rtype_code`_'atype_code` ('rtype_name` x, 'atype_name` n)
{
'rtype_name` pow = 1;
for (;;)
{
if (n & 1)
pow *= x;
n >>= 1;
if (n)
x *= x;
else
break;
}
return pow;
}
/* For odd x, Euler''`s theorem tells us that x**(2^(m-1)) = 1 mod 2^m.
For even x, we use the fact that (2*x)^m = 0 mod 2^m. */
'rtype_name`
pow_'rtype_code`_'atype_code` ('rtype_name` x, 'atype_name` n)
{
const 'rtype_name` mask = ('rtype_name`) (-1) / 2;
if (n == 0)
return 1;
if (x == 0)
return 0;
if (x & 1)
return power_simple_'rtype_code`_'atype_code` (x, n & mask);
if (n > sizeof (x) * 8)
return 0;
return power_simple_'rtype_code`_'atype_code` (x, n);
}
#endif'