Acme-Tools
view release on metacpan or search on metacpan
mins
maxs
sum
avg
geomavg
harmonicavg
stddev
rstddev
median
percentile
$Resolve_iterations
$Resolve_last_estimate
$Resolve_time
resolve
resolve_equation
conv
rank
rankstr
egrep
eqarr
sorted
B<Input:> 1-6 arguments. At least one argument.
First argument: must be a coderef to a subroutine (a function)
Second argument: if present, the target, f(x)=target. Default 0.
Third argument: a start position for x. Default 0.
Fourth argument: a small delta value. Default 1e-4 (0.0001).
Fifth argument: a maximum number of iterations before resolve gives up
and carps. Default 100 (if fifth argument is not given or is
undef). The number 0 means infinite here. If the derivative of the
start position is zero or close to zero more iterations are typically
needed.
Sixth argument: A number of seconds to run before giving up. If both
fifth and sixth argument is given and > 0, C<resolve> stops at
whichever comes first.
B<Output:> returns the number C<x> for C<f(x)> = 0
...or equal to the second input argument if present.
B<Example:>
The equation C<< x^2 - 4x - 21 = 0 >> has two solutions: -3 and 7.
The result of C<resolve> will depend on the start position:
print resolve(sub{ $_**2 - 4*$_ - 21 }); # -3 with $_ as your x
print resolve(sub{ my $x=shift; $x**2 - 4*$x - 21 }); # -3 more elaborate call
print resolve(sub{ my $x=shift; $x**2 - 4*$x - 21 },0,3); # 7 with start position 3
print "Iterations: $Acme::Tools::Resolve_iterations\n"; # 3 or larger, about 10-15 is normal
The variable C< $Acme::Tools::Resolve_iterations > (which is exported) will be set
to the last number of iterations C<resolve> used. Also if C<resolve> dies (carps).
The variable C< $Acme::Tools::Resolve_last_estimate > (which is exported) will be
set to the last estimate. This number will often be close to the solution and can
be used even if C<resolve> dies (carps).
B<BigFloat-example:>
If either second, third or fourth argument is an instance of L<Math::BigFloat>, so will the result be:
use Acme::Tools;
See:
L<http://en.wikipedia.org/wiki/Newtons_method>
L<Math::BigFloat>
L<http://en.wikipedia.org/wiki/Golden_ratio>
=cut
our $Resolve_iterations;
our $Resolve_last_estimate;
our $Resolve_time;
#sub resolve(\[&$]@) {
#sub resolve(&@) { <=0.17
#todo: perl -MAcme::Tools -le'print resolve(sub{$_[0]**2-9431**2});print$Acme::Tools::Resolve_iterations'
#todo: perl -MAcme::Tools -le'sub d{5.3*1.0094**$_[0]-10.2*1.0072**$_[0]} print resolve(\&d)' #err, pop norway vs sweden
#todo: perl -MAcme::Tools -le' print resolve(sub{5.3*1.0094**$_[0]-10.2*1.0072**$_[0]})' #err, pop norway vs sweden
# =>Div by zero: df(x) = 0 at n'th iteration, n=0, delta=0.0001, fx=CODE(0xc81d470) at -e line 1
#todo: ren solve?
sub resolve {
my($f,$goal,$start,$delta,$iters,$sec)=@_;
$goal=0 if!defined$goal;
$start=0 if!defined$start;
$delta=1e-4 if!defined$delta;
$iters=100 if!defined$iters;
$sec=0 if!defined$sec;
$iters=13e13 if $iters==0;
croak "Iterations ($iters) or seconds ($sec) can not be a negative number" if $iters<0 or $sec<0;
$Resolve_iterations=undef;
$Resolve_last_estimate=undef;
croak "Should have at least 1 argument, a coderef" if !@_;
croak "First argument should be a coderef" if ref($f) ne 'CODE';
my @x=($start);
my $time_start=$sec>0?time_fp():undef;
my $ds=ref($start) eq 'Math::BigFloat' ? Math::BigFloat->div_scale() : undef;
my $fx=sub{
local$_=$_[0];
my $fx=&$f($_);
#warn "delta=$delta\n";
my $n=0;
while($n<=$iters-1){
my $fd= &$fx($x[$n]+$delta*0.5) - &$fx($x[$n]-$delta*0.5);
$fd = &$fx($x[$n]+$delta*0.7) - &$fx($x[$n]-$delta*0.3) if $fd==0;# and warn"wigle 1\n";
$fd = &$fx($x[$n]+$delta*0.2) - &$fx($x[$n]-$delta*0.8) if $fd==0;# and warn"wigle 2\n";
croak "Div by zero: df(x) = $x[$n] at n'th iteration, n=$n, delta=$delta, fx=$fx" if $fd==0;
$x[$n+1]=$x[$n]-(&$fx($x[$n])-$goal)/($fd/$delta);
$Resolve_last_estimate=$x[$n+1];
#warn "n=$n fd=$fd x=$x[$n+1]\n";
$Resolve_iterations=$n;
last if $n>3 and $x[$n+1]==$x[$n] and $x[$n]==$x[$n-1];
last if $n>4 and $x[$n]!=0 and abs(1-$x[$n+1]/$x[$n])<1e-13; #sub{3*$_+$_**4-12}
last if $n>3 and ref($x[$n+1]) eq 'Math::BigFloat' and substr($x[$n+1],0,$ds) eq substr($x[$n],0,$ds); #hm
croak "Could not resolve, perhaps too little time given ($sec), iteratons=$n"
if $sec>0 and ($Resolve_time=time_fp()-$time_start)>$sec;
#warn "$n: ".$x[$n+1]."\n";
$n++;
}
croak "Could not resolve, perhaps too few iterations ($iters)" if @x>=$iters;
return $x[-1];
}
=head2 resolve_equation
This prints 2:
print resolve_equation "x + 13*(3-x) = 17 - x"
A string containing at least one x is converted into a perl function.
t/04_resolve.t view on Meta::CPAN
# perl Makefile.PL && make && perl -Iblib/lib t/04_resolve.t
use lib '.'; BEGIN{require 't/common.pl'}
use Test::More tests => 17;
if($ENV{ATDEBUG}){
deb "Resolve: ".resolve(sub{my($x)=(@_); $x**2 - 4*$x -1},20,2)."\n";
deb "Resolve: ".resolve(sub{my($x)=@_; $x**log($x)-$x},0,3)."\n";
deb "Resolve: ".resolve(sub{$_[0]})." iters=$Acme::Tools::Resolve_iterations\n";
}
my $e;
ok(resolve(sub{my($x)=@_; $x**2 - 4*$x -21}) == -3 ,'first solution');
ok(($e=resolve(sub{ $_**2 - 4*$_ - 21 })) == -3 ,"first solution with \$_ (=$e)");
ok(resolve(sub{$_**2 - 4*$_ -21},0,3) == 7 ,'second solution, start 3');
ok(resolve(sub{my($x)=@_; $x**2 - 4*$x -21},0,2) == 7 ,'second solution, start 2');
my $f=sub{ $_**2 - 4*$_ - 21 };
ok(do{my$r=resolve($f,0,2); $r== 7} ,'second solution, start 2');
ok(resolve($f,0,2) == 7 ,'second solution, start 2');
ok(resolve($f,0,2) == 7 ,'second solution, start 2');
ok($Resolve_iterations > 1 ,"iterations=$Resolve_iterations");
ok($Resolve_last_estimate == 7 ,"last_estimate=$Resolve_last_estimate (should be 7)");
eval{ resolve(sub{1}) }; # 1=0
ok($@=~/Div by zero/);
ok(!defined $Resolve_iterations);
ok(!defined $Resolve_last_estimate);
my $c;
eval{$e=resolve(sub{$c++; sleep_fp(0.02); $_**2 - 4*$_ -21},0,.02,undef,undef,0.05)};
deb "x=$e, est=$Resolve_last_estimate, iters=$Resolve_iterations, time=$Resolve_time, c=$c -- $@\n";
ok($@=~/Could not resolve, perhaps too little time given/,'ok $@');
my$no=0;sub isr{is( ($e=$_[0]), $_[1], "r".(++$no).": e=$e, iters=$Resolve_iterations")}
isr( sprintf("%.12f",resolve(sub{3*$_ + $_**4 - 12})), '1.632498783713' ); #*)
isr( log(resolve(sub{ $_**log($_)-$_},0,2)), 1);
isr( resolve(sub{$_**2+7*$_-60},0,1), 5);
isr( resolve_equation("x^2+7x-60"), 5);
#*) http://www.quickmath.com/webMathematica3/quickmath/equations/solve/basic.jsp#c=solve_stepssolveequation&v1=3x%2Bx%5E4-12%3D0&v2=x
( run in 0.576 second using v1.01-cache-2.11-cpan-96521ef73a4 )