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);
}