A weird recurrence relation
Author: Shlomi Fish
http://projecteuler.net/problem=463
The function f
is defined for all positive integers as follows:
f(1)=1 f(3)=3 f(2n)=f(n) f(4n + 1)=2f(2n + 1) - f(n) f(4n + 3)=3f(2n + 1) - 2f(n)
The function S(n)
is defined as ∑_i=1..q f(i)
.
S(8)=22 and S(100)=3604.
Find S(3**37)
. Give the last 9 digits of your answer.
Dedicated to the memory of Christina Grimmie (https://en.wikipedia.org/wiki/Christina_Grimmie), see: https://www.youtube.com/watch?v=kYX8sjIzjGw. -- shlomif
Source code: prob463-shlomif.p6
use v6; my $debug = False; class Poly { has $.mult; # Coefficients. has @.c; method inc() { return Poly.new(mult => $.mult, c => [@.c[0]+@.c[1], @.c[1]]); } method double() { return Poly.new(mult => $.mult, c => [@.c[0], @.c[1]*2]); } } my %cache; { my $mult = 4; my @polys = ( Poly.new(mult => 6, c => [1,2]), Poly.new(mult => -2, c => [0,1]), ); %cache = ("$mult" => @polys); } sub lookup_proto($mult) { return %cache{$mult} //= sub { my $polys = lookup_proto($mult / 2); my @extend = @$polys.map( { .inc }); my @new_polys = flat(@$polys, @extend).map({ .double }); my $x2_p_1_coeff = 0; my $x_coeff = 0; for @new_polys -> $p { my $m = $p.mult; my @c = @($p.c); while (@c[1] > 1 and @c[0] % 2 == 0) { for @c -> $x is rw { $x /= 2; } } if @c[1] == 1 { if @c[0] != 0 { die "Dont know how to handle."; } $x_coeff += $m; } elsif @c[1] == 2 { if @c[0] != 1 { die "Dont know how to handle."; } $x2_p_1_coeff += $m; } elsif @c[1] == 4 { if @c[0] == 1 { $x2_p_1_coeff += ($m +< 1); $x_coeff -= $m; } elsif @c[0] == 3 { $x2_p_1_coeff += ($m * 3); $x_coeff -= ($m +< 1); } else { die "Dont know how to handle."; } } else { die "Dont know how to handle."; } } my @ret = ( Poly.new(mult => $x2_p_1_coeff, c => [1,2]), Poly.new(mult => $x_coeff, c => [0,1]), ); return @ret; }.(); } sub lookup($mult) { my $ret = lookup_proto($mult); return @($ret).map({ .mult }); } if $debug { my $mult = 4; while (1) { say lookup($mult); $mult +<= 1; } } my $MOD = 1_000_000_000; sub _cache(%h, $key, $promise) { my $ret = %h{$key}; if !defined($ret) { $ret = ($promise.() % $MOD); %h{$key} = $ret; } return $ret; } my %f_cache; sub f_mod($n) { if $n < 1 { die "Foo"; } return _cache(%f_cache, $n, sub { if $n == 1 { return 1; } elsif $n == 3 { return 3; } elsif ($n +& 1) == 0 { return f_mod($n +> 1); } elsif ($n +& 3) == 1 { return (2 * f_mod(($n +> 1) + 1) - f_mod($n +> 2)); } else { return (3 * f_mod(($n +> 1)) - 2 * f_mod($n +> 2)); } } ); } sub s_bruteforce { my ($n) = @_; my $s = 0; for 1 .. $n -> $i { ($s += f_mod($i)) %= $MOD; } return $s; } sub s_smart($start, $end) { state %s_cache; # say "s->e : $start->$end"; return _cache(%s_cache, "$start|$end", sub { if $start > $end { return 0; } if $start == $end { return f_mod($start); } if $end <= 8 { return s_bruteforce($end) - s_bruteforce($start - 1); } if ($start +& 0b11) != 0 { return (f_mod($start) + s_smart($start+1, $end)); } if ($end +& 0b11) != 0b11 { return (f_mod($end) + s_smart($start, $end-1)); } my $power2 = ((($start +& ($start-1)) != 0) ?? 1+(($start-1)+^$start) !! $start); my $new_end = $start + $power2 - 1; while ($new_end > $end) { $new_end = $start + ($power2 +>= 1) - 1; } my @c = lookup($power2); my $m = $start / $power2; return (@c[0] * f_mod($m*2+1) + @c[1] * f_mod($m) + s_smart($new_end+1, $end)); }, ); } if $debug { my $want = 0; for 1 .. 100_000 -> $n { if $n % 1_000 == 0 { say "Reached n=$n"; } ($want += f_mod($n)) %= $MOD; my $have = s_smart(1, $n); if $want != $have { die "want=$want have=$have n=$n!"; } } } { say "S(3 ** 37) = ", s_smart(1, 3 ** 37); }