1package Math::Utils; 2 3use 5.010001; 4use strict; 5use warnings; 6use Carp; 7 8use Exporter; 9our @ISA = qw(Exporter); 10 11our %EXPORT_TAGS = ( 12 compare => [ qw(generate_fltcmp generate_relational) ], 13 fortran => [ qw(log10 copysign) ], 14 utility => [ qw(log10 log2 copysign flipsign 15 sign floor ceil fsum 16 gcd hcf lcm moduli softmax 17 uniform_scaling uniform_01scaling) ], 18 polynomial => [ qw(pl_evaluate pl_dxevaluate pl_translate 19 pl_add pl_sub pl_div pl_mult 20 pl_derivative pl_antiderivative) ], 21); 22 23our @EXPORT_OK = ( 24 @{ $EXPORT_TAGS{compare} }, 25 @{ $EXPORT_TAGS{utility} }, 26 @{ $EXPORT_TAGS{polynomial} }, 27); 28 29# 30# Add an :all tag automatically. 31# 32$EXPORT_TAGS{all} = [@EXPORT_OK]; 33 34our $VERSION = '1.14'; 35 36=head1 NAME 37 38Math::Utils - Useful mathematical functions not in Perl. 39 40=head1 SYNOPSIS 41 42 use Math::Utils qw(:utility); # Useful functions 43 44 # 45 # Base 10 and base 2 logarithms. 46 # 47 $scale = log10($pagewidth); 48 $bits = log2(1/$probability); 49 50 # 51 # Two uses of sign(). 52 # 53 $d = sign($z - $w); 54 55 @ternaries = sign(@coefficients); 56 57 # 58 # Using copysign(), $dist will be doubled negative or 59 # positive $offest, depending upon whether ($from - $to) 60 # is positive or negative. 61 # 62 my $dist = copysign(2 * $offset, $from - $to); 63 64 # 65 # Change increment direction if goal is negative. 66 # 67 $incr = flipsign($incr, $goal); 68 69 # 70 # floor() and ceil() functions. 71 # 72 $point = floor($goal); 73 $limit = ceil($goal); 74 75 # 76 # gcd() and lcm() functions. 77 # 78 $divisor = gcd(@multipliers); 79 $numerator = lcm(@multipliers); 80 81 # 82 # Safer summation. 83 # 84 $tot = fsum(@inputs); 85 86 # 87 # The remainders of n after successive divisions of b, or 88 # remainders after a set of divisions. 89 # 90 @rems = moduli($n, $b); 91 92or 93 94 use Math::Utils qw(:compare); # Make comparison functions with tolerance. 95 96 # 97 # Floating point comparison function. 98 # 99 my $fltcmp = generate_fltmcp(1.0e-7); 100 101 if (&$fltcmp($x0, $x1) < 0) 102 { 103 add_left($data); 104 } 105 else 106 { 107 add_right($data); 108 } 109 110 # 111 # Or we can create single-operation comparison functions. 112 # 113 # Here we are only interested in the greater than and less than 114 # comparison functions. 115 # 116 my(undef, undef, 117 $approx_gt, undef, $approx_lt) = generate_relational(1.5e-5); 118 119or 120 121 use Math::Utils qw(:polynomial); # Basic polynomial ops 122 123 # 124 # Coefficient lists run from 0th degree upward, left to right. 125 # 126 my @c1 = (1, 3, 5, 7, 11, 13, 17, 19); 127 my @c2 = (1, 3, 1, 7); 128 my @c3 = (1, -1, 1) 129 130 my $c_ref = pl_mult(\@c1, \@c2); 131 $c_ref = pl_add($c_ref, \@c3); 132 133=head1 EXPORT 134 135All functions can be exported by name, or by using the tag that they're 136grouped under. 137 138=cut 139 140=head2 utility tag 141 142Useful, general-purpose functions, including those that originated in 143FORTRAN and were implemented in Perl in the module L<Math::Fortran>, 144by J. A. R. Williams. 145 146There is a name change -- copysign() was known as sign() 147in Math::Fortran. 148 149=head3 log10() 150 151 $xlog10 = log10($x); 152 @xlog10 = log10(@x); 153 154Return the log base ten of the argument. A list form of the function 155is also provided. 156 157=cut 158 159sub log10 160{ 161 my $log10 = log(10); 162 return wantarray? map(log($_)/$log10, @_): log($_[0])/$log10; 163} 164 165=head3 log2() 166 167 $xlog2 = log2($x); 168 @xlog2 = log2(@x); 169 170Return the log base two of the argument. A list form of the function 171is also provided. 172 173=cut 174 175sub log2 176{ 177 my $log2 = log(2); 178 return wantarray? map(log($_)/$log2, @_): log($_[0])/$log2; 179} 180 181=head3 sign() 182 183 $s = sign($x); 184 @valsigns = sign(@values); 185 186Returns -1 if the argument is negative, 0 if the argument is zero, and 1 187if the argument is positive. 188 189In list form it applies the same operation to each member of the list. 190 191=cut 192 193sub sign 194{ 195 return wantarray? map{($_ < 0)? -1: (($_ > 0)? 1: 0)} @_: 196 ($_[0] < 0)? -1: (($_[0] > 0)? 1: 0); 197} 198 199=head3 copysign() 200 201 $ms = copysign($m, $n); 202 $s = copysign($x); 203 204Take the sign of the second argument and apply it to the first. Zero 205is considered part of the positive signs. 206 207 copysign(-5, 0); # Returns 5. 208 copysign(-5, 7); # Returns 5. 209 copysign(-5, -7); # Returns -5. 210 copysign(5, -7); # Returns -5. 211 212If there is only one argument, return -1 if the argument is negative, 213otherwise return 1. For example, copysign(1, -4) and copysign(-4) both 214return -1. 215 216=cut 217 218sub copysign 219{ 220 return ($_[1] < 0)? -abs($_[0]): abs($_[0]) if (@_ == 2); 221 return ($_[0] < 0)? -1: 1; 222} 223 224=head3 flipsign() 225 226 $ms = flipsign($m, $n); 227 228Multiply the signs of the arguments and apply it to the first. As 229with copysign(), zero is considered part of the positive signs. 230 231Effectively this means change the sign of the first argument if 232the second argument is negative. 233 234 flipsign(-5, 0); # Returns -5. 235 flipsign(-5, 7); # Returns -5. 236 flipsign(-5, -7); # Returns 5. 237 flipsign(5, -7); # Returns -5. 238 239If for some reason flipsign() is called with a single argument, 240that argument is returned unchanged. 241 242=cut 243 244sub flipsign 245{ 246 return -$_[0] if (@_ == 2 and $_[1] < 0); 247 return $_[0]; 248} 249 250=head3 floor() 251 252 $b = floor($a/2); 253 254 @ilist = floor(@numbers); 255 256Returns the greatest integer less than or equal to its argument. 257A list form of the function also exists. 258 259 floor(1.5, 1.87, 1); # Returns (1, 1, 1) 260 floor(-1.5, -1.87, -1); # Returns (-2, -2, -1) 261 262=cut 263 264sub floor 265{ 266 return wantarray? map(($_ < 0 and int($_) != $_)? int($_ - 1): int($_), @_): 267 ($_[0] < 0 and int($_[0]) != $_[0])? int($_[0] - 1): int($_[0]); 268} 269 270=head3 ceil() 271 272 $b = ceil($a/2); 273 274 @ilist = ceil(@numbers); 275 276Returns the lowest integer greater than or equal to its argument. 277A list form of the function also exists. 278 279 ceil(1.5, 1.87, 1); # Returns (2, 2, 1) 280 ceil(-1.5, -1.87, -1); # Returns (-1, -1, -1) 281 282=cut 283 284sub ceil 285{ 286 return wantarray? map(($_ > 0 and int($_) != $_)? int($_ + 1): int($_), @_): 287 ($_[0] > 0 and int($_[0]) != $_[0])? int($_[0] + 1): int($_[0]); 288} 289 290=head3 fsum() 291 292Return a sum of the values in the list, done in a manner to avoid rounding 293and cancellation errors. Currently this is done via 294L<Kahan's summation algorithm|https://en.wikipedia.org/wiki/Kahan_summation_algorithm>. 295 296=cut 297 298sub fsum 299{ 300 my($sum, $c) = (0, 0); 301 302 for my $v (@_) 303 { 304 my $y = $v - $c; 305 my $t = $sum + $y; 306 307 # 308 # If we lost low-order bits of $y (usually because 309 # $sum is much larger than $y), save them in $c 310 # for the next loop iteration. 311 # 312 $c = ($t - $sum) - $y; 313 $sum = $t; 314 } 315 316 return $sum; 317} 318 319=head3 softmax() 320 321Return a list of values as probabilities. 322 323The function takes the list, and creates a new list by raising I<e> to 324each value. The function then returns each value divided by the sum of 325the list. Each value in the new list is now a set of probabilities that 326sum to 1.0. 327 328The summation is performed using I<fsum()> above. 329 330See L<Softmax function|https://en.wikipedia.org/wiki/Softmax_function> at 331Wikipedia. 332 333=cut 334 335sub softmax 336{ 337 my @nlist = @_; 338 339 # 340 # There's a nice trick where you find the maximum value in 341 # the list, and subtract it from every number in the list. 342 # This renders everything zero or negative, which makes 343 # exponentation safe from overflow, but doesn't affect 344 # the end result. 345 # 346 # If we weren't using this trick, then we'd start with 347 # the 'my @explist' line, feeding it '@_' instead. 348 # 349 my $listmax = $nlist[0]; 350 for (@nlist[1 .. $#nlist]) 351 { 352 $listmax = $_ if ($_ > $listmax); 353 } 354 @nlist = map{$_ - $listmax} @nlist if ($listmax > 0); 355 356 my @explist = map{exp($_)} @nlist; 357 my $sum = fsum(@explist); 358 return map{$_/$sum} @explist; 359} 360 361=head3 uniform_scaling 362 363=head3 uniform_01scaling 364 365Uniformly, or linearly, scale a number either from one range to another range 366(C<uniform_scaling()>), or to a default range of [0 .. 1] 367(C<uniform_01scaling()>). 368 369 @v = uniform_scaling(\@original_range, \@new_range, @oldvalues); 370 371For example, these two lines are equivalent, and both return 0: 372 373 $y = uniform_scaling([50, 100], [0, 1], 50); 374 375 $y = uniform_01scaling([50, 100], 50); 376 377They may also be called with a list or array of numbers: 378 379 @cm_measures = uniform_scaling([0, 10000], [0, 25400], @in_measures); 380 381 @melt_centigrade = uniform_scaling([0, 2000], [-273.15, 1726.85], \@melting_points); 382 383A number that is outside the original bounds will be proportionally changed 384to be outside of the new bounds, but then again having a number outside the 385original bounds is probably an error that should be checked before calling 386this function. 387 388L<https://stats.stackexchange.com/q/281164> 389 390=cut 391 392sub uniform_scaling 393{ 394 my @fromrange = @{$_[0]}; 395 my @torange = @{$_[1]}; 396 397 # 398 # The remaining parameters are the numbers to rescale. 399 # 400 # It could happen. Someone might type \$x instead of $x. 401 # 402 my @xvalues = map{(ref $_ eq "ARRAY")? @$_: 403 ((ref $_ eq "SCALAR")? $$_: $_)} @_[2 .. $#_]; 404 405 return map{($_ - $fromrange[0])/($fromrange[1] - $fromrange[0]) * ($torange[1] - $torange[0]) + $torange[0]} @xvalues; 406} 407 408sub uniform_01scaling 409{ 410 my @fromrange = @{$_[0]}; 411 412 # 413 # The remaining parameters are the numbers to rescale. 414 # 415 # It could happen. Someone might type \$x instead of $x. 416 # 417 my @xvalues = map{(ref $_ eq "ARRAY")? @$_: 418 ((ref $_ eq "SCALAR")? $$_: $_)} @_[1 .. $#_]; 419 420 return map{($_ - $fromrange[0]) / ($fromrange[1] - $fromrange[0])} @xvalues; 421} 422 423=head3 gcd 424 425=head3 hcf 426 427Return the greatest common divisor (also known as the highest 428common factor) of a list of integers. These are simply synomyms: 429 430 $factor = gcd(@numbers); 431 $factor = hcf(@numbers); 432 433=cut 434 435sub gcd 436{ 437 use integer; 438 my($x, $y, $r); 439 440 # 441 # It could happen. Someone might type \$x instead of $x. 442 # 443 my @values = map{(ref $_ eq "ARRAY")? @$_: 444 ((ref $_ eq "SCALAR")? $$_: $_)} grep {$_} @_; 445 446 return 0 if (scalar @values == 0); 447 448 $y = abs pop @values; 449 $x = abs pop @values; 450 451 while (1) 452 { 453 ($x, $y) = ($y, $x) if ($y < $x); 454 455 $r = $y % $x; 456 $y = $x; 457 458 if ($r == 0) 459 { 460 return $x if (scalar @values == 0); 461 $r = abs pop @values; 462 } 463 464 $x = $r; 465 } 466 467 return $y; 468} 469 470# 471#sub bgcd 472#{ 473# my($x, $y) = map(abs($_), @_); 474# 475# return $y if ($x == 0); 476# return $x if ($y == 0); 477# 478# my $lsbx = low_set_bit($x); 479# my $lsby = low_set_bit($y); 480# $x >>= $lsbx; 481# $y >>= $lsby; 482# 483# while ($x != $y) 484# { 485# ($x, $y) = ($y, $x) if ($x > $y); 486# 487# $y -= $x; 488# $y >>= low_set_bit($y); 489# } 490# return ($x << (($lsbx > $lsby)? $lsby: $lsbx)); 491#} 492 493*hcf = \&gcd; 494 495=head3 lcm 496 497Return the least common multiple of a list of integers. 498 499 $factor = lcm(@values); 500 501=cut 502 503sub lcm 504{ 505 # 506 # It could happen. Someone might type \$x instead of $x. 507 # 508 my @values = map{(ref $_ eq "ARRAY")? @$_: 509 ((ref $_ eq "SCALAR")? $$_: $_)} @_; 510 511 my $x = pop @values; 512 513 for my $m (@values) 514 { 515 $x *= $m/gcd($m, $x); 516 } 517 518 return abs $x; 519} 520 521=head3 moduli() 522 523Return the moduli of an integer after repeated divisions. The remainders are 524returned in a list from left to right. 525 526 @digits = moduli(1899, 10); # Returns (9, 9, 8, 1) 527 @rems = moduli(29, 3); # Returns (2, 0, 0, 1) 528 529=cut 530 531sub moduli 532{ 533 my($n, $b) = (abs($_[0]), abs($_[1])); 534 my @mlist; 535 use integer; 536 537 for (;;) 538 { 539 push @mlist, $n % $b; 540 $n /= $b; 541 return @mlist if ($n == 0); 542 } 543 return (); 544} 545 546=head2 compare tag 547 548Create comparison functions for floating point (non-integer) numbers. 549 550Since exact comparisons of floating point numbers tend to be iffy, 551the comparison functions use a tolerance chosen by you. You may 552then use those functions from then on confident that comparisons 553will be consistent. 554 555If you do not provide a tolerance, a default tolerance of 1.49012e-8 556(approximately the square root of an Intel Pentium's 557L<machine epsilon|https://en.wikipedia.org/wiki/Machine_epsilon>) 558will be used. 559 560=head3 generate_fltcmp() 561 562Returns a comparison function that will compare values using a tolerance 563that you supply. The generated function will return -1 if the first 564argument compares as less than the second, 0 if the two arguments 565compare as equal, and 1 if the first argument compares as greater than 566the second. 567 568 my $fltcmp = generate_fltcmp(1.5e-7); 569 570 my(@xpos) = grep {&$fltcmp($_, 0) == 1} @xvals; 571 572=cut 573 574my $default_tolerance = 1.49012e-8; 575 576sub generate_fltcmp 577{ 578 my $tol = $_[0] // $default_tolerance; 579 580 return sub { 581 my($x, $y) = @_; 582 return 0 if (abs($x - $y) <= $tol); 583 return -1 if ($x < $y); 584 return 1; 585 } 586} 587 588=head3 generate_relational() 589 590Returns a list of comparison functions that will compare values using a 591tolerance that you supply. The generated functions will be the equivalent 592of the equal, not equal, greater than, greater than or equal, less than, 593and less than or equal operators. 594 595 my($eq, $ne, $gt, $ge, $lt, $le) = generate_relational(1.5e-7); 596 597 my(@approx_5) = grep {&$eq($_, 5)} @xvals; 598 599Of course, if you were only interested in not equal, you could use: 600 601 my(undef, $ne) = generate_relational(1.5e-7); 602 603 my(@not_around5) = grep {&$ne($_, 5)} @xvals; 604 605=cut 606 607sub generate_relational 608{ 609 my $tol = $_[0] // $default_tolerance; 610 611 # 612 # In order: eq, ne, gt, ge, lt, le. 613 # 614 return ( 615 sub {return (abs($_[0] - $_[1]) <= $tol)? 1: 0;}, # eq 616 sub {return (abs($_[0] - $_[1]) > $tol)? 1: 0;}, # ne 617 618 sub {return ((abs($_[0] - $_[1]) > $tol) and ($_[0] > $_[1]))? 1: 0;}, # gt 619 sub {return ((abs($_[0] - $_[1]) <= $tol) or ($_[0] > $_[1]))? 1: 0;}, # ge 620 621 sub {return ((abs($_[0] - $_[1]) > $tol) and ($_[0] < $_[1]))? 1: 0;}, # lt 622 sub {return ((abs($_[0] - $_[1]) <= $tol) or ($_[0] < $_[1]))? 1: 0;} # le 623 ); 624} 625 626=head2 polynomial tag 627 628Perform some polynomial operations on plain lists of coefficients. 629 630 # 631 # The coefficient lists are presumed to go from low order to high: 632 # 633 @coefficients = (1, 2, 4, 8); # 1 + 2x + 4x**2 + 8x**3 634 635In all functions the coeffcient list is passed by reference to the function, 636and the functions that return coefficients all return references to a 637coefficient list. 638 639B<It is assumed that any leading zeros in the coefficient lists have 640already been removed before calling these functions, and that any leading 641zeros found in the returned lists will be handled by the caller.> This caveat 642is particularly important to note in the case of C<pl_div()>. 643 644Although these functions are convenient for simple polynomial operations, 645for more advanced polynonial operations L<Math::Polynomial> is recommended. 646 647=head3 pl_evaluate() 648 649Returns either a y-value for a corresponding x-value, or a list of 650y-values on the polynomial for a corresponding list of x-values, 651using Horner's method. 652 653 $y = pl_evaluate(\@coefficients, $x); 654 @yvalues = pl_evaluate(\@coefficients, @xvalues); 655 656 @ctemperatures = pl_evaluate([-160/9, 5/9], @ftemperatures); 657 658The list of X values may also include X array references: 659 660 @yvalues = pl_evaluate(\@coefficients, @xvalues, \@primes, $x, [-1, -10, -100]); 661 662=cut 663 664sub pl_evaluate 665{ 666 my @coefficients = @{$_[0]}; 667 668 # 669 # It could happen. Someone might type \$x instead of $x. 670 # 671 my @xvalues = map{(ref $_ eq "ARRAY")? @$_: 672 ((ref $_ eq "SCALAR")? $$_: $_)} @_[1 .. $#_]; 673 674 # 675 # Move the leading coefficient off the polynomial list 676 # and use it as our starting value(s). 677 # 678 my @results = (pop @coefficients) x scalar @xvalues; 679 680 for my $c (reverse @coefficients) 681 { 682 for my $j (0..$#xvalues) 683 { 684 $results[$j] = $results[$j] * $xvalues[$j] + $c; 685 } 686 } 687 688 return wantarray? @results: $results[0]; 689} 690 691=head3 pl_dxevaluate() 692 693 ($y, $dy, $ddy) = pl_dxevaluate(\@coefficients, $x); 694 695Returns p(x), p'(x), and p"(x) of the polynomial for an 696x-value, using Horner's method. Note that unlike C<pl_evaluate()> 697above, the function can only use one x-value. 698 699If the polynomial is a linear equation, the second derivative value 700will be zero. Similarly, if the polynomial is a simple constant, 701the first derivative value will be zero. 702 703=cut 704 705sub pl_dxevaluate 706{ 707 my($coef_ref, $x) = @_; 708 my(@coefficients) = @$coef_ref; 709 my $n = $#coefficients; 710 my $val = pop @coefficients; 711 my $d1val = $val * $n; 712 my $d2val = 0; 713 714 # 715 # Special case for the linear eq'n (the y = constant eq'n 716 # takes care of itself). 717 # 718 if ($n == 1) 719 { 720 $val = $val * $x + $coefficients[0]; 721 } 722 elsif ($n >= 2) 723 { 724 my $lastn = --$n; 725 $d2val = $d1val * $n; 726 727 # 728 # Loop through the coefficients, except for 729 # the linear and constant terms. 730 # 731 for my $c (reverse @coefficients[2..$lastn]) 732 { 733 $val = $val * $x + $c; 734 $d1val = $d1val * $x + ($c *= $n--); 735 $d2val = $d2val * $x + ($c * $n); 736 } 737 738 # 739 # Handle the last two coefficients. 740 # 741 $d1val = $d1val * $x + $coefficients[1]; 742 $val = ($val * $x + $coefficients[1]) * $x + $coefficients[0]; 743 } 744 745 return ($val, $d1val, $d2val); 746} 747 748=head3 pl_translate() 749 750 $x = [8, 3, 1]; 751 $y = [3, 1]; 752 753 # 754 # Translating C<x**2 + 3*x + 8> by C<x + 3> returns [26, 9, 1] 755 # 756 $z = pl_translate($x, $y); 757 758Returns a polynomial transformed by substituting a polynomial variable with another polynomial. 759For example, a simple linear translation by 1 to the polynomial C<x**3 + x**2 + 4*x + 4> 760would be accomplished by setting x = (y - 1); resulting in C<x**3 - 2*x**2 + 5*x>. 761 762 $x = [4, 4, 1, 1]; 763 $y = [-1, 1]; 764 $z = pl_translate($x, $y); # Becomes [0, 5, -2, 1] 765 766=cut 767 768sub pl_translate 769{ 770 my($x, $y) = @_; 771 772 my @x_arr = @$x; 773 my @z = pop @x_arr; 774 775 for my $c (reverse @x_arr) 776 { 777 @z = @{ pl_mult(\@z, $y) }; 778 $z[0] += $c; 779 } 780 781 return [@z]; 782} 783 784=head3 pl_add() 785 786 $polyn_ref = pl_add(\@m, \@n); 787 788Add two lists of numbers as though they were polynomial coefficients. 789 790=cut 791 792sub pl_add 793{ 794 my(@av) = @{$_[0]}; 795 my(@bv) = @{$_[1]}; 796 my $ldiff = scalar @av - scalar @bv; 797 798 my @result = ($ldiff < 0)? 799 splice(@bv, scalar @bv + $ldiff, -$ldiff): 800 splice(@av, scalar @av - $ldiff, $ldiff); 801 802 unshift @result, map($av[$_] + $bv[$_], 0.. $#av); 803 804 return \@result; 805} 806 807=head3 pl_sub() 808 809 $polyn_ref = pl_sub(\@m, \@n); 810 811Subtract the second list of numbers from the first as though they 812were polynomial coefficients. 813 814=cut 815 816sub pl_sub 817{ 818 my(@av) = @{$_[0]}; 819 my(@bv) = @{$_[1]}; 820 my $ldiff = scalar @av - scalar @bv; 821 822 my @result = ($ldiff < 0)? 823 map {-$_} splice(@bv, scalar @bv + $ldiff, -$ldiff): 824 splice(@av, scalar @av - $ldiff, $ldiff); 825 826 unshift @result, map($av[$_] - $bv[$_], 0.. $#av); 827 828 return \@result; 829} 830 831=head3 pl_div() 832 833 ($q_ref, $r_ref) = pl_div(\@numerator, \@divisor); 834 835Synthetic division for polynomials. Divides the first list of coefficients 836by the second list. 837 838Returns references to the quotient and the remainder. 839 840Remember to check for leading zeros (which are rightmost in the list) in 841the returned values. For example, 842 843 my @n = (4, 12, 9, 3); 844 my @d = (1, 3, 3, 1); 845 846 my($q_ref, $r_ref) = pl_div(\@n, \@d); 847 848After division you will have returned C<(3)> as the quotient, 849and C<(1, 3, 0)> as the remainder. In general, you will want to remove 850the leading zero, or for that matter values within epsilon of zero, in 851the remainder. 852 853 my($q_ref, $r_ref) = pl_div($f1, $f2); 854 855 # 856 # Remove any leading zeros (i.e., numbers smaller in 857 # magnitude than machine epsilon) in the remainder. 858 # 859 my @remd = @{$r_ref}; 860 pop @remd while (@remd and abs($remd[$#remd]) < $epsilon); 861 862 $f1 = $f2; 863 $f2 = [@remd]; 864 865If C<$f1> and C<$f2> were to go through that bit of code again, not 866removing the leading zeros would lead to a divide-by-zero error. 867 868If either list of coefficients is empty, pl_div() returns undefs for 869both quotient and remainder. 870 871=cut 872 873sub pl_div 874{ 875 my @numerator = @{$_[0]}; 876 my @divisor = @{$_[1]}; 877 878 my @quotient; 879 880 my $n_degree = $#numerator; 881 my $d_degree = $#divisor; 882 883 # 884 # Sanity checks: a numerator less than the divisor 885 # is automatically the remainder; and return a pair 886 # of undefs if either set of coefficients are 887 # empty lists. 888 # 889 return ([0], \@numerator) if ($n_degree < $d_degree); 890 return (undef, undef) if ($d_degree < 0 or $n_degree < 0); 891 892 my $lead_coefficient = $divisor[$#divisor]; 893 894 # 895 # Perform the synthetic division. The remainder will 896 # be what's left in the numerator. 897 # (4, 13, 4, -9, 6) / (1, 2) = (4, 5, -6, 3) 898 # 899 @quotient = reverse map { 900 # 901 # Get the next term for the quotient. We pop 902 # off the lead numerator term, which would become 903 # zero due to subtraction anyway. 904 # 905 my $q = (pop @numerator)/$lead_coefficient; 906 907 for my $k (0..$d_degree - 1) 908 { 909 $numerator[$#numerator - $k] -= $q * $divisor[$d_degree - $k - 1]; 910 } 911 912 $q; 913 } reverse (0 .. $n_degree - $d_degree); 914 915 return (\@quotient, \@numerator); 916} 917 918=head3 pl_mult() 919 920 $m_ref = pl_mult(\@coefficients1, \@coefficients2); 921 922Returns the reference to the product of the two multiplicands. 923 924=cut 925 926sub pl_mult 927{ 928 my($av, $bv) = @_; 929 my $a_degree = $#{$av}; 930 my $b_degree = $#{$bv}; 931 932 # 933 # Rather than multiplying left to right for each element, 934 # sum to each degree of the resulting polynomial (the list 935 # after the map block). Still an O(n**2) operation, but 936 # we don't need separate storage variables. 937 # 938 return [ map { 939 my $a_idx = ($a_degree > $_)? $_: $a_degree; 940 my $b_to = ($b_degree > $_)? $_: $b_degree; 941 my $b_from = $_ - $a_idx; 942 943 my $c = $av->[$a_idx] * $bv->[$b_from]; 944 945 for my $b_idx ($b_from+1 .. $b_to) 946 { 947 $c += $av->[--$a_idx] * $bv->[$b_idx]; 948 } 949 $c; 950 } (0 .. $a_degree + $b_degree) ]; 951} 952 953=head3 pl_derivative() 954 955 $poly_ref = pl_derivative(\@coefficients); 956 957Returns the derivative of a polynomial. 958 959=cut 960 961sub pl_derivative 962{ 963 my @coefficients = @{$_[0]}; 964 my $degree = $#coefficients; 965 966 return [] if ($degree < 1); 967 968 $coefficients[$_] *= $_ for (2..$degree); 969 970 shift @coefficients; 971 return \@coefficients; 972} 973 974=head3 pl_antiderivative() 975 976 $poly_ref = pl_antiderivative(\@coefficients); 977 978Returns the antiderivative of a polynomial. The constant value is 979always set to zero and will need to be changed by the caller if a 980different constant is needed. 981 982 my @coefficients = (1, 2, -3, 2); 983 my $integral = pl_antiderivative(\@coefficients); 984 985 # 986 # Integral needs to be 0 at x = 1. 987 # 988 my @coeff1 = @{$integral}; 989 $coeff1[0] = - pl_evaluate($integral, 1); 990 991=cut 992 993sub pl_antiderivative 994{ 995 my @coefficients = @{$_[0]}; 996 my $degree = scalar @coefficients; 997 998 # 999 # Sanity check if its an empty list. 1000 # 1001 return [0] if ($degree < 1); 1002 1003 $coefficients[$_ - 1] /= $_ for (2..$degree); 1004 1005 unshift @coefficients, 0; 1006 return \@coefficients; 1007} 1008 1009=head1 AUTHOR 1010 1011John M. Gamble, C<< <jgamble at cpan.org> >> 1012 1013=head1 SEE ALSO 1014 1015L<Math::Polynomial> for a complete set of polynomial operations, with the 1016added convenience that objects bring. 1017 1018Among its other functions, L<List::Util> has the mathematically useful 1019functions max(), min(), product(), sum(), and sum0(). 1020 1021L<List::MoreUtils> has the function minmax(). 1022 1023L<Math::Prime::Util> has gcd() and lcm() functions, as well as vecsum(), 1024vecprod(), vecmin(), and vecmax(), which are like the L<List::Util> 1025functions but which can force integer use, and when appropriate use 1026L<Math::BigInt>. 1027 1028L<Math::VecStat> Likewise has min(), max(), sum() (which can take 1029as arguments array references as well as arrays), plus maxabs(), 1030minabs(), sumbyelement(), convolute(), and other functions. 1031 1032=head1 BUGS 1033 1034Please report any bugs or feature requests to C<bug-math-util at rt.cpan.org>, or through 1035the web interface at L<http://rt.cpan.org/NoAuth/ReportBug.html?Queue=Math-Utils>. I will be notified, and then you'll 1036automatically be notified of progress on your bug as I make changes. 1037 1038=head1 SUPPORT 1039 1040This module is on Github at L<https://github.com/jgamble/Math-Utils>. 1041 1042You can also look for information at: 1043 1044=over 4 1045 1046=item * RT: CPAN's request tracker (report bugs here) 1047 1048L<http://rt.cpan.org/NoAuth/Bugs.html?Dist=Math-Utils> 1049 1050=item * MetaCPAN 1051 1052L<https://metacpan.org/release/Math-Utils> 1053 1054=back 1055 1056=head1 ACKNOWLEDGEMENTS 1057 1058To J. A. R. Williams who got the ball rolling with L<Math::Fortran>. 1059 1060=head1 LICENSE AND COPYRIGHT 1061 1062Copyright (c) 2017 John M. Gamble. All rights reserved. This program is 1063free software; you can redistribute it and/or modify it under the same 1064terms as Perl itself. 1065 1066=cut 1067 10681; # End of Math::Utils 1069