view release on metacpan or  search on metacpan

devel/ema-omitted-old.pl  view on Meta::CPAN

    }
    if ($i < 0) {
      last;
    }

    push @ret, ($coeff > 0 ? $CONFIG{PLUS}.$coeff : $CONFIG{MINUS}.-$coeff);
    $need_parens = 1;
  }
  $ret[0] = '(' x $open_parens;
  return join('', @ret);
}

package main;


#------------------------------------------------------------------------------

sub print_ma_terms {
  my ($name) = @_;
  print "$name\n";
  my $term_proc = do { no strict; \&$name };
  foreach my $k (0 .. 20) {
    say "  p[$k] * (@{[$term_proc->($k)]})";
  }
}

my $poly_0 = Math::Polynomial->new;

# {
#   say Math::Polynomial->new()->to_string_for_eval;
#   say Math::Polynomial->new(1)->to_string_for_eval;
#   say Math::Polynomial->new(-1)->to_string_for_eval;
#   say Math::Polynomial->new(1,0)->to_string_for_eval;
#   say Math::Polynomial->new(-1,0)->to_string_for_eval;
#   say Math::Polynomial->new(-1,0,1)->to_string_for_eval;
#   say Math::Polynomial->new(123,0,-456)->to_string_for_eval;
#   say Math::Polynomial->new(-1,-1,-456)->to_string_for_eval;
#   say Math::Polynomial->new(1,1,1,1,1)->to_string_for_eval;
#   say Math::Polynomial->new(1,1,1,1,0)->to_string_for_eval;
#   say Math::Polynomial->new(1,1,1,0,0)->to_string_for_eval;
#   say Math::Polynomial->new(0,1,1,1,0,0)->to_string_for_eval;
#   exit 0;
# }

# per Knuth vol 2 sect 4.3.3
sub interpolate_successive {
  my $start = 0;
  if (! ref $_[0] && $_[0] eq 'start') {
    shift;
    $start = shift;
  }

  my @data = @_;
  my @divisor;
  my $divisor = 1;
  foreach my $i (1 .. $#data) {
    push @divisor, ($divisor *= $i);
    for (my $j = $#data; $j >= $i; $j--) {
      $data[$j] -= $data[$j-1];
    }
    if (DEBUG) { require Data::Dumper;
                 print Data::Dumper->Dump([\@data],['data']); }
  }
  if (DEBUG) { require Data::Dumper;
               print Data::Dumper->Dump([\@divisor],['divisor']); }

  my $ret = Math::Polynomial->new((pop @data) / (pop @divisor));
  while (@data > 2) {
    $ret->mul1c ($#data + $start);
    $ret->[-1] += (pop @data) / (pop @divisor);
  }
  # $data[1]
  $ret->mul1c (1 + $start);
  $ret->[-1] += pop @data;
  # $data[0]
  if ($start) { $ret->mul1c ($start) }
  push @$ret, pop @data;

  return $ret;
}

# sub interpolate_faster {
#   my (@x, @y);
#   while (@_) {
#     push @x, shift;
#     push @y, shift;
#   }
#   my @nums;
#   my $numerator = Math::Polynomial->new(1);
#   for (my $i = $#x; $i > 0; $i--) {
#     $numerator->mul1c($x[$i]);
#     push @nums, $numerator->clone;
#   }
# 
#   my $result = Math::Polynomial->new(0);
#   foreach my $i (0 .. $#x-1) {
#     $numerator = shift @nums;
#         print "num ",$numerator,"\n";
#   my $constant = pop(@y) / $numerator->eval($_);
#     $result += $numerator * $constant;
#   }
#   return $result;
# }


{
  say interpolate_successive(10, 304, 1980, 7084, 18526);
  say interpolate_successive(0,1,4);
  say interpolate_successive(1,2,3);
  say interpolate_successive(1, 4, 9, 16);

  my $triangle = sub {
    my ($n) = @_;
    return List::Util::sum (1 .. $n);
  };
  my $sum_squares = sub {
    my ($n) = @_;
    return List::Util::sum (map {$_**2} (1 .. $n));
  };
  my $yden = sub {
    my ($n) = @_;
    return $n * $sum_squares->($n) - ($triangle->($n))**2;
  };
  require Math::Polynomial;
  require Math::BigRat;