swampyankee has asked for the
wisdom of the Perl Monks concerning the following question:
In another forum (one that's science, not software, oriented), there was a poster asking for programming advice. He's trying to port some code to a PLC, the development environment of which can't do exponentiation.
Alas, he needs to do noninteger exponents. So, my question for the Mathematical Monks is does anybody know of an exponentiation function? Because of the way exponens work, the only concern is fractional powers, i.e. those with absolute values less than one. I'm digging through my copy of Knuth and Numerical Recipes, but since real languages can do $x = $y ** $z; with some restrictions, I've never looked into that deeply into the numerical libraries.
Information about American English usage here and here. Floating point issues? Please read this before posting. — emc
Re: [OT] Function to do x**y for y<1 w/o ** operator. by BrowserUk (Pope) on Feb 19, 2013 at 19:06 UTC 
print exp( log( 100 ) *.5 );;
10
print exp( log( 1000 ) *.3333333333333333333333 );;
10
Ie. sub pow{ exp( log( $_[0] ) * $_[1] ) };;
printf "%2d : %20.17f\n", $_, pow( 10**$_, 1/$_ ) for 1 .. 20;;
1 : 10.00000000000000200
2 : 10.00000000000000200
3 : 9.99999999999999820
4 : 10.00000000000000200
5 : 10.00000000000000200
6 : 9.99999999999999820
7 : 9.99999999999999820
8 : 10.00000000000000200
9 : 9.99999999999999820
10 : 10.00000000000000200
11 : 10.00000000000000200
12 : 9.99999999999999820
13 : 10.00000000000000200
14 : 9.99999999999999820
15 : 9.99999999999999820
16 : 10.00000000000000200
17 : 10.00000000000000200
18 : 9.99999999999999820
19 : 10.00000000000000200
20 : 10.00000000000000200
printf "%2d : %20.17f\n", $_, pow( 12345**$_, 1/$_ ) for 1 .. 20;;
1 : 12345.00000000000500000
2 : 12345.00000000000500000
3 : 12345.00000000000500000
4 : 12345.00000000000500000
5 : 12345.00000000000500000
6 : 12345.00000000000500000
7 : 12345.00000000000500000
8 : 12345.00000000000500000
9 : 12345.00000000000500000
10 : 12345.00000000000500000
11 : 12345.00000000000500000
12 : 12345.00000000000500000
13 : 12345.00000000000500000
14 : 12345.00000000000500000
15 : 12344.99999999998400000
16 : 12345.00000000000500000
17 : 12344.99999999998400000
18 : 12345.00000000000500000
19 : 12344.99999999998400000
20 : 12345.00000000000500000
With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
Examine what is said, not who speaks  Silence betokens consent  Love the truth but pardon error.
"Science is about questioning the status quo. Questioning authority".
In the absence of evidence, opinion is indistinguishable from prejudice.
 [reply] [d/l] [select] 

He has log, but not exp.
Information about American English usage here and here. Floating point issues? Please read this before posting. — emc
 [reply] 

#! perl slw
use strict;
sub myExp{
my $e = 1 + ( shift() / 10000000 );
my $r = $e;
$r *= $e for 2 .. 10000000;
$r;
}
printf "%30.20g : %30.20g\n", exp( $_ ), myExp( $_ )
for map log( 10**$_ ), 10 .. 10;
__END__
C:\test>exp
9.9999999999999965e011 : 9.9997349036884958e011
1.0000000000000007e009 : 9.9997852777366009e010
9.999999999999982e009 : 9.9998303391048947e009
9.9999999999999943e008 : 9.9998700988200641e008
1.0000000000000004e006 : 9.9999045679358269e007
9.9999999999999974e006 : 9.9999337241854084e006
0.00010000000000000009 : 9.9999575897369183e005
0.0010000000000000002 : 0.00099999761423453593
0.010000000000000005 : 0.0099999893930861859
0.10000000000000002 : 0.099999973530420908
1 : 1
10.000000000000002 : 9.9999973561653093
100.00000000000004 : 99.999893882330838
999.99999999999977 : 999.99761406302275
10000.000000000009 : 9999.9575911775646
100000.00000000001 : 99999.337176962814
999999.99999999953 : 999990.45646420098
10000000.000000006 : 9999870.1098560225
100000000.00000018 : 99998303.311974689
999999999.99999928 : 999978527.31280696
10000000000.000004 : 9999734913.5273266
With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
Examine what is said, not who speaks  Silence betokens consent  Love the truth but pardon error.
"Science is about questioning the status quo. Questioning authority".
In the absence of evidence, opinion is indistinguishable from prejudice.
 [reply] [d/l] 


 
Update:Changed constants for slightly greater accuracy for the sake of another couple of iterations.
Here's a somewhat more accurate, and much faster version:
#! perl slw
use strict;
sub myExp{
my $e = 1 + ( shift() / 4294967296 );
$e *= $e for 1 .. 32;
$e;
}
printf "%30.20g : %30.20g\n", exp( $_ ), myExp( $_ )
for map log( 10**$_ ), 10 .. 10;
__END__
C:\test>exp
9.9999999999999965e011 : 1.0000000285332257e010
1.0000000000000007e009 : 1.0000000829163549e009
9.999999999999982e009 : 1.0000002130576865e008
9.9999999999999943e008 : 9.9999976329917627e008
1.0000000000000004e006 : 9.9999969112813618e007
9.9999999999999974e006 : 9.9999977582604716e006
0.00010000000000000009 : 9.999998681099902e005
0.0010000000000000002 : 0.00099999984556405617
0.010000000000000005 : 0.0099999993405499293
0.10000000000000002 : 0.099999996702749588
1 : 1
10.000000000000002 : 9.9999952742988558
100.00000000000004 : 100.00000084944983
999.99999999999977 : 999.99959325720954
10000.000000000009 : 10000.000169889967
100000.00000000001 : 99999.966690555739
999999.99999999953 : 1000000.1401405634
10000000.000000006 : 9999996.1589527112
100000000.00000018 : 100000003.39779937
999999999.99999928 : 999999662.81900096
10000000000.000004 : 10000002874.419859
With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
Examine what is said, not who speaks  Silence betokens consent  Love the truth but pardon error.
"Science is about questioning the status quo. Questioning authority".
In the absence of evidence, opinion is indistinguishable from prejudice.
 [reply] [d/l] 

He has log, but not exp.
Perhaps there's also something helpful in Math::BigFloat's _pow() sub:
# Taylor:  u u^2 u^3 
# x ** y = 1 +   +  +  + ... 
# _ 1 1*2 1*2*3 _
where (I gather)
u = y * ln x
Cheers, Rob  [reply] [d/l] [select] 


 [reply] 
Re: [OT] Function to do x**y for y<1 w/o ** operator. by moritz (Cardinal) on Feb 19, 2013 at 21:43 UTC 
 [reply] 
Re: [OT] Function to do x**y for y<1 w/o ** operator. by salva (Abbot) on Feb 20, 2013 at 09:18 UTC 
If you have log available, you can use Newton method to easyly find the solution to the equation f(z) = log(z)  y * log x = 0 :
#!/usr/bin/perl
use strict;
use warnings;
my ($X, $Y, $n) = @ARGV;
$n = 20;
my $Z = $X ** $Y;
# f(z) = log z  y * log x = 0;
# f'(z) = f1(z) = 1/z
my $ylogx = $Y * log $X;
my $z = 1;
for my $i (1..$n) {
my $fz = log($z)  $ylogx;
my $f1z = 1 / $z;
$z = $z  $fz / $f1z;
my $e = abs($z  $Z);
printf "%3d: z=%8f, e=%8f\n", $i, $z, $e;
}
This method does not handle correctly values of y < 0 but this is not a problem because z(y) = 1 / z(y).
Also, note that even if the algorithm converges in very few iterations, it uses log that is usually an expensive operation and so, other methods requiring more iterations may be actually faster.  [reply] [d/l] [select] 
Re: [OT] Function to do x**y for y<1 w/o ** operator. by salva (Abbot) on Feb 20, 2013 at 11:37 UTC 
Another solution based on the fact that a floating point number is actually a sum of powers of 2.
For 0 <= y <= 0, it holds that y = b_{0}*2^{0} + b_{1}*2^{1} + b_{2}*2^{2} + b_{3}*2^{3} + ...
And so, x^{y} = x^{b0*20 + b1*21 + b2*22 + b3*23 + ...} = x^{b0*20} * x^{b1*21} * x^{b2*22} * x^{b3*23} * ... = Π_{i, bi≠0} x^{2i}
Which in Perl becomes...
#!/usr/bin/perl
use strict;
use warnings;
my ($X, $Y, $n) = @ARGV;
$n = 20;
my $Z = $X ** $Y;
my ($neg, $y) = ($Y < 0 ? (1, $Y) : (0, $Y));
$y > 1 and die "Y is out of range";
my $bit = 1;
my $x = $X;
my $z = 1;
while ($y) {
if ($y >= $bit) {
$y = $bit;
$z *= $x;
}
$bit /= 2;
$x = sqrt($x);
}
$z = 1/$z if $neg;
my $e = abs($z  $X ** $Y);
print "z=$z, e=$e\n";
sqrt is also an expensive operation. If you have to calculate x^{y} for different values of y while x stays constant, you may be able to speed up the process creating a table with the values of x^{2i}.
Alternatively, you can write a function to calculate e^{x} (using a table with the values of e^{2i}), and then calculate x^{y} as e^{log(x)*y}.
update: there are also several implementations of exp(x) freely available, for instance, the one in OpenBSD is here.  [reply] [d/l] [select] 

