Commits

Shlomi Fish  committed 8b4d59b

Add the euler-101 - is not working.

  • Participants
  • Parent commits 126743b

Comments (0)

Files changed (1)

File project-euler/101/euler-101.pl

+#!/usr/bin/perl
+
+use strict;
+use warnings;
+
+use Math::MatrixReal;
+
+use List::Util;
+
+=head1 PROBLEM
+
+If we are presented with the first k terms of a sequence it is impossible to say with certainty the value of the next term, as there are infinitely many polynomial functions that can model the sequence.
+
+As an example, let us consider the sequence of cube numbers. This is defined by the generating function,
+un = n3: 1, 8, 27, 64, 125, 216, ...
+
+Suppose we were only given the first two terms of this sequence. Working on the principle that "simple is best" we should assume a linear relationship and predict the next term to be 15 (common difference 7). Even if we were presented with the first three terms, by the same principle of simplicity, a quadratic relationship should be assumed.
+
+We shall define OP(k, n) to be the nth term of the optimum polynomial generating function for the first k terms of a sequence. It should be clear that OP(k, n) will accurately generate the terms of the sequence for n ≤ k, and potentially the first incorrect term (FIT) will be OP(k, k+1); in which case we shall call it a bad OP (BOP).
+
+As a basis, if we were only given the first term of sequence, it would be most sensible to assume constancy; that is, for n ≥ 2, OP(1, n) = u1.
+
+Hence we obtain the following OPs for the cubic sequence:
+OP(1, n) = 1 	1, 1, 1, 1, ...
+OP(2, n) = 7n−6 	1, 8, 15, ...
+OP(3, n) = 6n2−11n+6      	1, 8, 27, 58, ...
+OP(4, n) = n3 	1, 8, 27, 64, 125, ...
+
+Clearly no BOPs exist for k ≥ 4.
+
+By considering the sum of FITs generated by the BOPs (indicated in red above), we obtain 1 + 15 + 58 = 74.
+
+Consider the following tenth degree polynomial generating function:
+
+un = 1 − n + n2 − n3 + n4 − n5 + n6 − n7 + n8 − n9 + n10
+
+Find the sum of FITs for the BOPs.
+
+=head1 Planning.
+
+Suppose we have n terms for which the polynomial yields (a1, a2, a3 ... an).
+Then we can have:
+
+c_0 + c_1 * 1 + c_2 * 1 ** 2 ... + c_n-1 * 1 ** n-1 = a_1
+c_0 + c_1 * 2 + c_2 * 2 ** 2 ... + c_n-1 * 2 ** n-1 = a_2
+c_0 + c_1 * 3 + c_2 * 3 ** 2 ... + c_n-1 * 3 ** n-1 = a_3
+
+So we can find c_0 ... c_n-1 by solving the co-efficients problem.
+
+=cut
+
+sub find_coeff
+{
+    my ($a_s) = @_;
+
+    my $coeffs = Math::MatrixReal->new_from_rows(
+        [ map { my $x = $_; [map { $x ** $_ } (0 .. $#$a_s) ] } 
+            (1 .. scalar(@$a_s))
+        ]
+    );
+    
+    $coeffs->transpose($coeffs);
+    
+    my $a_mat = Math::MatrixReal->new_from_rows([map { [$_] } @$a_s]);
+
+    return $coeffs->inverse() * $a_mat;
+}
+
+sub get_bop
+{
+    my ($a_s) = @_;
+
+    return 
+    (
+        Math::MatrixReal->new_from_rows([[map { @$a_s ** $_ } 1 .. scalar(@$a_s) ]]) 
+        * find_coeff($a_s)
+    )->element(1,1);
+}
+
+my @u_coeffs = (map { (-1) ** $_ } (0 .. 10));
+
+sub calc_u_result_vec
+{
+    my $x = shift;
+
+    return [map { $x ** $_ * $u_coeffs[$_] } (0 .. $#u_coeffs)];
+}
+
+sub calc_u_result
+{
+    my $x = shift;
+    return List::Util::sum (@{calc_u_result_vec($x)});
+}
+
+my @u_results = (map { calc_u_result($_) } (1 .. scalar(@u_coeffs)));
+
+sub get_u_bop
+{
+    my $i = shift;
+
+    return get_bop( [ @u_results[0 .. $i-1] ]);
+}
+
+
+my $s = 0;
+foreach my $i (1 .. 10)
+{
+    my $val = get_u_bop($i);
+    print "$i : $val\n";
+    $s += $val;
+}
+print "$s\n";