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*2-1 + b2*2-2 + b3*2-3 + ...} = x^{b0*20} * x^{b1*2-1} * x^{b2*2-2} * x^{b3*2-3} * ... = Π_{i, bi≠0} x^{2-i}

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^{2-i}.

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.

Comment onRe: [OT] Function to do x**y for |y|<1 w/o ** operator.SelectorDownloadCode