Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Add ipow() method and a test in t/ops.t
ipow() calculates the integer exponentiation of a value
by successive multiplications.  This allows one to avoid
loss of significance in integer exponents.  pow() converts
the value to double and will always have <52 precision.
  • Loading branch information
devel-chm committed Oct 6, 2015
1 parent fca5668 commit 7e2a426
Show file tree
Hide file tree
Showing 2 changed files with 60 additions and 1 deletion.
50 changes: 50 additions & 0 deletions Basic/Ops/ops.pd
Expand Up @@ -408,6 +408,56 @@ pp_def(
'Plain numerical assignment. This is used to implement the ".=" operator',
); # pp_def assgn

pp_def('ipow',
Doc => qq{
=for ref

raise piddle C<\$a> to integer power C<\$b>

=for example

\$c = \$a->ipow(\$b,0); # explicit function call
\$c = ipow \$a, \$b;
\$a->inplace->ipow(\$b,0); # modify \$a inplace

It can be made to work inplace with the C<\$a-E<gt>inplace> syntax.
Note that when calling this function explicitly you need to supply
a third argument that should generally be zero (see first example).
This restriction is expected to go away in future releases.

Algorithm from L<Wikipedia|http://en.wikipedia.org/wiki/Exponentiation_by_squaring>

=cut

},
Pars => 'a(); b(); [o] ans()',
Code => pp_line_numbers(__LINE__, q{
PDL_Indx n = $b();
$GENERIC() y = 1;
$GENERIC() x = $a();
if (n < 0) {
x = 1 / x;
n = -n;
}
if (n == 0) {
$ans() = 1;
} else {
while (n > 1) {
if (n % 2 == 0) {
x = x * x;
n = n / 2;
} else {
y = x * y;
x = x * x;
n = (n - 1) / 2;
}
}
$ans() = x * y;
}
})
);


#pp_export_nothing();

pp_done();
Expand Down
11 changes: 10 additions & 1 deletion t/ops.t
@@ -1,4 +1,4 @@
use Test::More tests => 59;
use Test::More tests => 60;
use PDL::LiteF;
use Config;
kill INT,$$ if $ENV{UNDER_DEBUGGER}; # Useful for debugging.
Expand Down Expand Up @@ -135,6 +135,15 @@ ok(all($data == 1), 'or assign');

ok(all($data eq $data), 'eq'); # check eq operator

SKIP:
{
skip("ipow skipped without 64bit support", 1) if $Config{ivsize} < 8;
# check ipow routine
my $xdata = indx(0xeb * ones(8));
my $n = sequence(indx,8);
my $exact = indx(1,235,55225,12977875,3049800625,716703146875,168425239515625,39579931286171875);
ok(all($exact - ipow($xdata,$n) == indx(0)), 'ipow');
}

#### Modulus checks ####

Expand Down

0 comments on commit 7e2a426

Please sign in to comment.