Stirling numbers and bernoulli numbers for summing up n numbers to the kth power

Question

I am currently working on problem 487 on project euler. I did some research and I only see 2 possibilities to solve this problem: 1. By using faulhabers formula 2. by using the formula featuring the stirling numbers of the second kind (cant remember the formulas name) Because n is 10^12 and k is 10000 I think the only possibility is to use the 2. option and an approximation of the factorial. If there is a better option than 1 and 2 please let me know. If there isnt please help me what kind of approximation I should use for the factorial in order to be accurate enough. Thanks for your help and effort

Answer 1

So I wrote a solution to Euler487 in C#. I went about it in a rather brute force way, so I wanted a language that would run faster than python. I would've used C but I don't have a good IDE installed at the moment.

The range of primes to run through falls just under the maximum size of positive 32 bit signed int. Which means that its square falls just under the size of positive 64 bit signed long. Which makes these numbers a perfect candidate for exponentiation by squaring with modulation. Furthermore the summation of powers doesn't have to range through all trillion numbers since J[p] only has p representatives. This cuts computation down by a factor of about 500. I ran this code for 1 value of p and it took, approximately 17 min. There are 101 primes in the range specified, giving a synchronous run time of a little over a day. I'll let it run on my computer, just out of curiosity. There's got to be a more elegant solution out there. If you find one I would love to see it.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace Euler487
{
    class Program
    {
        static void Main(string[] args)
        {
            for (long p = 2000000000, Euler487 = 0; p < 2000002000; p = findNextPrime(p)) {
                Euler487 += sumPowersMod(1000000000000, 10000, p);
            }
            Console.ReadLine();
        }

        static long findNextPrime(long n)
        {
            n++;
            if (n % 2 == 0) n++;
            for (; !isPrime(n); n++) { }
            return n;
        }

        static bool isPrime(long n)
        {
            if (n % 2 == 0) return false;
            for (long k = n, d = 3; k >= d; d += 2)
            {
                k = n / d;
                if (n % d == 0) return false;
            }
            return true;
        }

        static long modExp(long number, long exponent, long modulus)
        {
            long result = 1;
            long currentPower = number;
            for (; exponent > 0; exponent >>= 1)
            {
                if ((exponent & 1) == 1)
                {
                    result *= currentPower;
                    result = result % modulus;
                }
                currentPower = (currentPower * currentPower) % modulus;
            }
            return result % modulus;
        }

        static long sumPowersMod(long n, long k, long p)
        {
            long quotient = n / p;
            long remainder = n % p;
            long totalUpToP = 0;
            long totalUpToRem = 0;
            for (long i = 1; i < p; i++)
            {
                totalUpToP += modExp(i, k, p);
                if (i == remainder) { totalUpToRem = totalUpToP; }
            }
            return (totalUpToP * quotient + totalUpToRem) % p;
        }
    }
}