Problem 776: Digit Sum Division

Problem 776 source code

C++

#include <cassert>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>

#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/multiprecision/cpp_int.hpp>

namespace {

using boost::multiprecision::cpp_dec_float_100;
using boost::multiprecision::cpp_int;

struct Stats {
    cpp_int count = 0;
    cpp_int sum = 0;
};

cpp_dec_float_100 F_of_u64(std::uint64_t N) {
    const std::string s = std::to_string(N);
    const int L = static_cast<int>(s.size());
    const int MAX_SUM = 9 * L;

    std::vector<Stats> tight(MAX_SUM + 1), loose(MAX_SUM + 1);
    tight[0].count = 1;

    for (int pos = 0; pos < L; ++pos) {
        const int limit = s[static_cast<std::size_t>(pos)] - '0';
        std::vector<Stats> ntight(MAX_SUM + 1), nloose(MAX_SUM + 1);

        for (int sum = 0; sum <= MAX_SUM; ++sum) {
            if (tight[sum].count != 0) {
                for (int d = 0; d <= limit; ++d) {
                    const int ns = sum + d;
                    Stats& dst = (d == limit) ? ntight[ns] : nloose[ns];
                    dst.count += tight[sum].count;
                    dst.sum += tight[sum].sum * 10 + tight[sum].count * d;
                }
            }

            if (loose[sum].count != 0) {
                for (int d = 0; d <= 9; ++d) {
                    const int ns = sum + d;
                    Stats& dst = nloose[ns];
                    dst.count += loose[sum].count;
                    dst.sum += loose[sum].sum * 10 + loose[sum].count * d;
                }
            }
        }

        tight.swap(ntight);
        loose.swap(nloose);
    }

    cpp_dec_float_100 ans = 0;
    for (int sum = 1; sum <= MAX_SUM; ++sum) {
        const cpp_int total_sum = tight[sum].sum + loose[sum].sum;
        if (total_sum == 0) {
            continue;
        }
        ans += cpp_dec_float_100(total_sum) / cpp_dec_float_100(sum);
    }
    return ans;
}

void assert_close(const cpp_dec_float_100& got, const cpp_dec_float_100& expected,
                  const cpp_dec_float_100& tol) {
    cpp_dec_float_100 diff = got - expected;
    if (diff < 0) {
        diff = -diff;
    }
    assert(diff <= tol);
}

}  // namespace

int main() {
    assert(F_of_u64(10) == cpp_dec_float_100(19));
    assert_close(F_of_u64(123), cpp_dec_float_100("1.187764610390e3"),
                 cpp_dec_float_100("1e-9"));
    assert_close(F_of_u64(12345), cpp_dec_float_100("4.855801996238e6"),
                 cpp_dec_float_100("1e-6"));

    const cpp_dec_float_100 ans = F_of_u64(1'234'567'890'123'456'789ULL);
    std::cout << std::scientific << std::setprecision(12) << ans << '\n';
    return 0;
}

Python

import decimal

def f_of_u64(N):
    decimal.getcontext().prec = 100
    s = str(N)
    L = len(s)
    max_sum = 9 * L
    
    tight = [{'count': 0, 'sum': 0} for _ in range(max_sum + 1)]
    loose = [{'count': 0, 'sum': 0} for _ in range(max_sum + 1)]
    
    tight[0]['count'] = 1
    
    for pos in range(L):
        limit = int(s[pos])
        ntight = [{'count': 0, 'sum': 0} for _ in range(max_sum + 1)]
        nloose = [{'count': 0, 'sum': 0} for _ in range(max_sum + 1)]
        
        for sum_val in range(max_sum + 1):
            t_cnt = tight[sum_val]['count']
            t_sum = tight[sum_val]['sum']
            if t_cnt != 0:
                for d in range(limit + 1):
                    ns = sum_val + d
                    if d == limit:
                        ntight[ns]['count'] += t_cnt
                        ntight[ns]['sum'] += t_sum * 10 + t_cnt * d
                    else:
                        nloose[ns]['count'] += t_cnt
                        nloose[ns]['sum'] += t_sum * 10 + t_cnt * d
                        
            l_cnt = loose[sum_val]['count']
            l_sum = loose[sum_val]['sum']
            if l_cnt != 0:
                for d in range(10):
                    ns = sum_val + d
                    nloose[ns]['count'] += l_cnt
                    nloose[ns]['sum'] += l_sum * 10 + l_cnt * d
                    
        tight = ntight
        loose = nloose
        
    ans = decimal.Decimal(0)
    for sum_val in range(1, max_sum + 1):
        total_sum = tight[sum_val]['sum'] + loose[sum_val]['sum']
        if total_sum == 0:
            continue
        ans += decimal.Decimal(total_sum) / decimal.Decimal(sum_val)
        
    # Return formatted as e.g. 4.855801996238e+06
    # Note that scientific format with 12 decimals exactly matches C++'s setprecision(12)
    s_val = f"{ans:.12e}"
    # Replace e+06 with e+06 so it looks natural, python does e+06, C++ does e+06 or e+006
    return s_val

def solve():
    return f_of_u64(1234567890123456789)

if __name__ == "__main__":
    print(solve())

Java

import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;

public class Euler776 {

    static class Stats {
        BigInteger count = BigInteger.ZERO;
        BigInteger sum = BigInteger.ZERO;
    }

    public static String fOfU64(long n) {
        String s = Long.toString(n);
        int L = s.length();
        int maxSum = 9 * L;

        Stats[] tight = new Stats[maxSum + 1];
        Stats[] loose = new Stats[maxSum + 1];
        for (int i = 0; i <= maxSum; i++) {
            tight[i] = new Stats();
            loose[i] = new Stats();
        }

        tight[0].count = BigInteger.ONE;
        BigInteger ten = BigInteger.TEN;

        for (int pos = 0; pos < L; ++pos) {
            int limit = s.charAt(pos) - '0';
            Stats[] ntight = new Stats[maxSum + 1];
            Stats[] nloose = new Stats[maxSum + 1];
            for (int i = 0; i <= maxSum; i++) {
                ntight[i] = new Stats();
                nloose[i] = new Stats();
            }

            for (int sum = 0; sum <= maxSum; ++sum) {
                if (tight[sum].count.compareTo(BigInteger.ZERO) != 0) {
                    for (int d = 0; d <= limit; ++d) {
                        int ns = sum + d;
                        Stats dst = (d == limit) ? ntight[ns] : nloose[ns];
                        dst.count = dst.count.add(tight[sum].count);
                        BigInteger addSum = tight[sum].sum.multiply(ten)
                                .add(tight[sum].count.multiply(BigInteger.valueOf(d)));
                        dst.sum = dst.sum.add(addSum);
                    }
                }

                if (loose[sum].count.compareTo(BigInteger.ZERO) != 0) {
                    for (int d = 0; d <= 9; ++d) {
                        int ns = sum + d;
                        Stats dst = nloose[ns];
                        dst.count = dst.count.add(loose[sum].count);
                        BigInteger addSum = loose[sum].sum.multiply(ten)
                                .add(loose[sum].count.multiply(BigInteger.valueOf(d)));
                        dst.sum = dst.sum.add(addSum);
                    }
                }
            }

            tight = ntight;
            loose = nloose;
        }

        MathContext mc = new MathContext(100);
        BigDecimal ans = BigDecimal.ZERO;

        for (int sum = 1; sum <= maxSum; ++sum) {
            BigInteger totalSum = tight[sum].sum.add(loose[sum].sum);
            if (totalSum.equals(BigInteger.ZERO))
                continue;

            BigDecimal term = new BigDecimal(totalSum).divide(new BigDecimal(sum), mc);
            ans = ans.add(term, mc);
        }

        // Format to scientific notation with 12 decimal places
        // The format %e inherently uses 6 decimal places by default, so %.12e is used
        return String.format(java.util.Locale.US, "%.12e", ans);
    }

    public static String solve() {
        return fOfU64(1234567890123456789L);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}