Skip to content

Instantly share code, notes, and snippets.

@knowitkodekalender
Created December 22, 2019 20:21
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save knowitkodekalender/65d7a798df1b121148bade6d50b29bde to your computer and use it in GitHub Desktop.
Save knowitkodekalender/65d7a798df1b121148bade6d50b29bde to your computer and use it in GitHub Desktop.

Harshadprimtall

Av: Didrik Pemmer Aalen

Et Harshadtall er definert som et positivt tall som er delelig med summen av alle sifferne i tallet. I noen tilfeller er summen av alle sifferne i tallet et primtall. Disse kaller vi Harshadprimtall.

Eksempel

1729 er et Harshadtall fordi 1 + 7 + 2 + 9 = 19 og 1729 % 19 = 0. Dette er også et Harshadprimtall, fordi 19 er et primtall.

1730 er ikke et Harshadtall fordi 1 + 7 + 3 + 0 = 11 og 1730 % 11 = 3.

Oppgave

Hvor mange tall fra 1 til og med 98765432 er Harshadprimtall?

@terjew
Copy link

terjew commented Jan 2, 2020

@michaelo @bobcat4

Jeg testet litt med pregenerering av tverrsummer (droppet constexpr for å gøre det litt mer "fair", pregenereringen tar jo uansett bare noen mikrosekunder). Med min siste versjon av loopingen ser det faktisk ut som det er mindre effektivt å bruke pregenererte tverrsummer enn å summere kun de sifrene som endrer seg.

Jeg har gjort det konfigurerbart fra kommandolinjen, så dere kan jo teste på deres systemer og se om dere observerer det samme.

Mine resultater er omtrent som følger:

./knowit23 --small --count 1000 
cl: 6600ms
gcc: 6870ms
clang: 6750ms

./knowit23 --small --count 1000 --precalc
cl: 6800ms
gcc: 7100ms
clang: 6800ms
#include <iostream>
#include <thread>
#include <algorithm>
#include <inttypes.h>
#include <climits>
#include <math.h>

#include "cxxopts.hpp"
#include "fastmod.h"

typedef uint64_t digit_t;

bool* buildPrimeSieve(size_t maxPrime, bool inverted = true)
{
  auto maxSquareRoot = (int)sqrt(maxPrime);
  bool* eliminated = new bool[maxPrime + 1]{ false };
  eliminated[0] = true;
  eliminated[1] = true;
  for (int j = 4; j <= maxPrime; j += 2)
  {
    eliminated[j] = true;
  }

  for (int i = 3; i <= maxPrime; i += 2)
  {
    if (!eliminated[i])
    {
      if (i <= maxSquareRoot)
      {
        for (int j = i * i; j <= maxPrime; j += i)
        {
          eliminated[j] = true;
        }
      }
    }
  }
  if (!inverted) {
    for (int i = 0; i < maxPrime; i++) eliminated[i] = !eliminated[i];
  }
  return eliminated;
}

inline void check(uint32_t digitsum, uint32_t& count, bool* prime, uint32_t i, uint64_t* M)
{
  if (prime[digitsum] && fastmod::is_divisible(i, M[digitsum]))
  {
    count++;
  }
}

inline void check2(uint32_t& count, bool* prime, uint64_t* M, uint32_t i, uint32_t digitsum)
{
  if (prime[digitsum] && fastmod::is_divisible(i, M[digitsum]))
  {
    count++;
  }
}

inline void check100rest0(uint32_t digitsum, uint32_t& count, bool* prime, uint32_t i, uint64_t* M)
{
  check2(count, prime, M, i + 1, digitsum + 1);
  check2(count, prime, M, i + 2, digitsum + 2);
  check2(count, prime, M, i + 3, digitsum + 3);
  check2(count, prime, M, i + 5, digitsum + 5);
  check2(count, prime, M, i + 7, digitsum + 7);
  check2(count, prime, M, i + 10, digitsum + 1);
  check2(count, prime, M, i + 11, digitsum + 2);
  check2(count, prime, M, i + 12, digitsum + 3);
  check2(count, prime, M, i + 14, digitsum + 5);
  check2(count, prime, M, i + 16, digitsum + 7);
  check2(count, prime, M, i + 20, digitsum + 2);
  check2(count, prime, M, i + 21, digitsum + 3);
  check2(count, prime, M, i + 23, digitsum + 5);
  check2(count, prime, M, i + 25, digitsum + 7);
  check2(count, prime, M, i + 29, digitsum + 11);
  check2(count, prime, M, i + 30, digitsum + 3);
  check2(count, prime, M, i + 32, digitsum + 5);
  check2(count, prime, M, i + 34, digitsum + 7);
  check2(count, prime, M, i + 38, digitsum + 11);
  check2(count, prime, M, i + 41, digitsum + 5);
  check2(count, prime, M, i + 43, digitsum + 7);
  check2(count, prime, M, i + 47, digitsum + 11);
  check2(count, prime, M, i + 49, digitsum + 13);
  check2(count, prime, M, i + 50, digitsum + 5);
  check2(count, prime, M, i + 52, digitsum + 7);
  check2(count, prime, M, i + 56, digitsum + 11);
  check2(count, prime, M, i + 58, digitsum + 13);
  check2(count, prime, M, i + 61, digitsum + 7);
  check2(count, prime, M, i + 65, digitsum + 11);
  check2(count, prime, M, i + 67, digitsum + 13);
  check2(count, prime, M, i + 70, digitsum + 7);
  check2(count, prime, M, i + 74, digitsum + 11);
  check2(count, prime, M, i + 76, digitsum + 13);
  check2(count, prime, M, i + 83, digitsum + 11);
  check2(count, prime, M, i + 85, digitsum + 13);
  check2(count, prime, M, i + 89, digitsum + 17);
  check2(count, prime, M, i + 92, digitsum + 11);
  check2(count, prime, M, i + 94, digitsum + 13);
  check2(count, prime, M, i + 98, digitsum + 17);
}

inline void check100rest1(uint32_t digitsum, uint32_t& count, bool* prime, uint32_t i, uint64_t* M)
{
  check2(count, prime, M, i + 0, digitsum + 0);
  check2(count, prime, M, i + 1, digitsum + 1);
  check2(count, prime, M, i + 2, digitsum + 2);
  check2(count, prime, M, i + 4, digitsum + 4);
  check2(count, prime, M, i + 6, digitsum + 6);
  check2(count, prime, M, i + 10, digitsum + 1);
  check2(count, prime, M, i + 11, digitsum + 2);
  check2(count, prime, M, i + 13, digitsum + 4);
  check2(count, prime, M, i + 15, digitsum + 6);
  check2(count, prime, M, i + 19, digitsum + 10);
  check2(count, prime, M, i + 20, digitsum + 2);
  check2(count, prime, M, i + 22, digitsum + 4);
  check2(count, prime, M, i + 24, digitsum + 6);
  check2(count, prime, M, i + 28, digitsum + 10);
  check2(count, prime, M, i + 31, digitsum + 4);
  check2(count, prime, M, i + 33, digitsum + 6);
  check2(count, prime, M, i + 37, digitsum + 10);
  check2(count, prime, M, i + 39, digitsum + 12);
  check2(count, prime, M, i + 40, digitsum + 4);
  check2(count, prime, M, i + 42, digitsum + 6);
  check2(count, prime, M, i + 46, digitsum + 10);
  check2(count, prime, M, i + 48, digitsum + 12);
  check2(count, prime, M, i + 51, digitsum + 6);
  check2(count, prime, M, i + 55, digitsum + 10);
  check2(count, prime, M, i + 57, digitsum + 12);
  check2(count, prime, M, i + 60, digitsum + 6);
  check2(count, prime, M, i + 64, digitsum + 10);
  check2(count, prime, M, i + 66, digitsum + 12);
  check2(count, prime, M, i + 73, digitsum + 10);
  check2(count, prime, M, i + 75, digitsum + 12);
  check2(count, prime, M, i + 79, digitsum + 16);
  check2(count, prime, M, i + 82, digitsum + 10);
  check2(count, prime, M, i + 84, digitsum + 12);
  check2(count, prime, M, i + 88, digitsum + 16);
  check2(count, prime, M, i + 91, digitsum + 10);
  check2(count, prime, M, i + 93, digitsum + 12);
  check2(count, prime, M, i + 97, digitsum + 16);
  check2(count, prime, M, i + 99, digitsum + 18);

}


inline void check100rest2(uint32_t digitsum, uint32_t& count, bool* prime, uint32_t i, uint64_t* M)
{
  check2(count, prime, M, i + 0, digitsum + 0);
  check2(count, prime, M, i + 1, digitsum + 1);
  check2(count, prime, M, i + 3, digitsum + 3);
  check2(count, prime, M, i + 5, digitsum + 5);
  check2(count, prime, M, i + 9, digitsum + 9);
  check2(count, prime, M, i + 10, digitsum + 1);
  check2(count, prime, M, i + 12, digitsum + 3);
  check2(count, prime, M, i + 14, digitsum + 5);
  check2(count, prime, M, i + 18, digitsum + 9);
  check2(count, prime, M, i + 21, digitsum + 3);
  check2(count, prime, M, i + 23, digitsum + 5);
  check2(count, prime, M, i + 27, digitsum + 9);
  check2(count, prime, M, i + 29, digitsum + 11);
  check2(count, prime, M, i + 30, digitsum + 3);
  check2(count, prime, M, i + 32, digitsum + 5);
  check2(count, prime, M, i + 36, digitsum + 9);
  check2(count, prime, M, i + 38, digitsum + 11);
  check2(count, prime, M, i + 41, digitsum + 5);
  check2(count, prime, M, i + 45, digitsum + 9);
  check2(count, prime, M, i + 47, digitsum + 11);
  check2(count, prime, M, i + 50, digitsum + 5);
  check2(count, prime, M, i + 54, digitsum + 9);
  check2(count, prime, M, i + 56, digitsum + 11);
  check2(count, prime, M, i + 63, digitsum + 9);
  check2(count, prime, M, i + 65, digitsum + 11);
  check2(count, prime, M, i + 69, digitsum + 15);
  check2(count, prime, M, i + 72, digitsum + 9);
  check2(count, prime, M, i + 74, digitsum + 11);
  check2(count, prime, M, i + 78, digitsum + 15);
  check2(count, prime, M, i + 81, digitsum + 9);
  check2(count, prime, M, i + 83, digitsum + 11);
  check2(count, prime, M, i + 87, digitsum + 15);
  check2(count, prime, M, i + 89, digitsum + 17);
  check2(count, prime, M, i + 90, digitsum + 9);
  check2(count, prime, M, i + 92, digitsum + 11);
  check2(count, prime, M, i + 96, digitsum + 15);
  check2(count, prime, M, i + 98, digitsum + 17);
}
inline void check100rest3(uint32_t digitsum, uint32_t& count, bool* prime, uint32_t i, uint64_t* M)
{
  check2(count, prime, M, i + 0, digitsum + 0);
  check2(count, prime, M, i + 2, digitsum + 2);
  check2(count, prime, M, i + 4, digitsum + 4);
  check2(count, prime, M, i + 8, digitsum + 8);
  check2(count, prime, M, i + 11, digitsum + 2);
  check2(count, prime, M, i + 13, digitsum + 4);
  check2(count, prime, M, i + 17, digitsum + 8);
  check2(count, prime, M, i + 19, digitsum + 10);
  check2(count, prime, M, i + 20, digitsum + 2);
  check2(count, prime, M, i + 22, digitsum + 4);
  check2(count, prime, M, i + 26, digitsum + 8);
  check2(count, prime, M, i + 28, digitsum + 10);
  check2(count, prime, M, i + 31, digitsum + 4);
  check2(count, prime, M, i + 35, digitsum + 8);
  check2(count, prime, M, i + 37, digitsum + 10);
  check2(count, prime, M, i + 40, digitsum + 4);
  check2(count, prime, M, i + 44, digitsum + 8);
  check2(count, prime, M, i + 46, digitsum + 10);
  check2(count, prime, M, i + 53, digitsum + 8);
  check2(count, prime, M, i + 55, digitsum + 10);
  check2(count, prime, M, i + 59, digitsum + 14);
  check2(count, prime, M, i + 62, digitsum + 8);
  check2(count, prime, M, i + 64, digitsum + 10);
  check2(count, prime, M, i + 68, digitsum + 14);
  check2(count, prime, M, i + 71, digitsum + 8);
  check2(count, prime, M, i + 73, digitsum + 10);
  check2(count, prime, M, i + 77, digitsum + 14);
  check2(count, prime, M, i + 79, digitsum + 16);
  check2(count, prime, M, i + 80, digitsum + 8);
  check2(count, prime, M, i + 82, digitsum + 10);
  check2(count, prime, M, i + 86, digitsum + 14);
  check2(count, prime, M, i + 88, digitsum + 16);
  check2(count, prime, M, i + 91, digitsum + 10);
  check2(count, prime, M, i + 95, digitsum + 14);
  check2(count, prime, M, i + 97, digitsum + 16);
}

inline void check100rest4(uint32_t digitsum, uint32_t& count, bool* prime, uint32_t i, uint64_t* M)
{
  check2(count, prime, M, i + 1, digitsum + 1);
  check2(count, prime, M, i + 3, digitsum + 3);
  check2(count, prime, M, i + 7, digitsum + 7);
  check2(count, prime, M, i + 9, digitsum + 9);
  check2(count, prime, M, i + 10, digitsum + 1);
  check2(count, prime, M, i + 12, digitsum + 3);
  check2(count, prime, M, i + 16, digitsum + 7);
  check2(count, prime, M, i + 18, digitsum + 9);
  check2(count, prime, M, i + 21, digitsum + 3);
  check2(count, prime, M, i + 25, digitsum + 7);
  check2(count, prime, M, i + 27, digitsum + 9);
  check2(count, prime, M, i + 30, digitsum + 3);
  check2(count, prime, M, i + 34, digitsum + 7);
  check2(count, prime, M, i + 36, digitsum + 9);
  check2(count, prime, M, i + 43, digitsum + 7);
  check2(count, prime, M, i + 45, digitsum + 9);
  check2(count, prime, M, i + 49, digitsum + 13);
  check2(count, prime, M, i + 52, digitsum + 7);
  check2(count, prime, M, i + 54, digitsum + 9);
  check2(count, prime, M, i + 58, digitsum + 13);
  check2(count, prime, M, i + 61, digitsum + 7);
  check2(count, prime, M, i + 63, digitsum + 9);
  check2(count, prime, M, i + 67, digitsum + 13);
  check2(count, prime, M, i + 69, digitsum + 15);
  check2(count, prime, M, i + 70, digitsum + 7);
  check2(count, prime, M, i + 72, digitsum + 9);
  check2(count, prime, M, i + 76, digitsum + 13);
  check2(count, prime, M, i + 78, digitsum + 15);
  check2(count, prime, M, i + 81, digitsum + 9);
  check2(count, prime, M, i + 85, digitsum + 13);
  check2(count, prime, M, i + 87, digitsum + 15);
  check2(count, prime, M, i + 90, digitsum + 9);
  check2(count, prime, M, i + 94, digitsum + 13);
  check2(count, prime, M, i + 96, digitsum + 15);
}

inline void check100rest5(uint32_t digitsum, uint32_t& count, bool* prime, uint32_t i, uint64_t* M)
{
  check2(count, prime, M, i + 0, digitsum + 0);
  check2(count, prime, M, i + 2, digitsum + 2);
  check2(count, prime, M, i + 6, digitsum + 6);
  check2(count, prime, M, i + 8, digitsum + 8);
  check2(count, prime, M, i + 11, digitsum + 2);
  check2(count, prime, M, i + 15, digitsum + 6);
  check2(count, prime, M, i + 17, digitsum + 8);
  check2(count, prime, M, i + 20, digitsum + 2);
  check2(count, prime, M, i + 24, digitsum + 6);
  check2(count, prime, M, i + 26, digitsum + 8);
  check2(count, prime, M, i + 33, digitsum + 6);
  check2(count, prime, M, i + 35, digitsum + 8);
  check2(count, prime, M, i + 39, digitsum + 12);
  check2(count, prime, M, i + 42, digitsum + 6);
  check2(count, prime, M, i + 44, digitsum + 8);
  check2(count, prime, M, i + 48, digitsum + 12);
  check2(count, prime, M, i + 51, digitsum + 6);
  check2(count, prime, M, i + 53, digitsum + 8);
  check2(count, prime, M, i + 57, digitsum + 12);
  check2(count, prime, M, i + 59, digitsum + 14);
  check2(count, prime, M, i + 60, digitsum + 6);
  check2(count, prime, M, i + 62, digitsum + 8);
  check2(count, prime, M, i + 66, digitsum + 12);
  check2(count, prime, M, i + 68, digitsum + 14);
  check2(count, prime, M, i + 71, digitsum + 8);
  check2(count, prime, M, i + 75, digitsum + 12);
  check2(count, prime, M, i + 77, digitsum + 14);
  check2(count, prime, M, i + 80, digitsum + 8);
  check2(count, prime, M, i + 84, digitsum + 12);
  check2(count, prime, M, i + 86, digitsum + 14);
  check2(count, prime, M, i + 93, digitsum + 12);
  check2(count, prime, M, i + 95, digitsum + 14);
  check2(count, prime, M, i + 99, digitsum + 18);
}

inline void addToDigit(digit_t* digits, uint32_t& digitsum, uint32_t digitno, digit_t amount)
{
  digits[digitno] += amount;
  digit_t diff = 0;
  diff += amount;
  if (digits[digitno] >= 10) {
    digits[digitno] -= 10;
    diff -= 9;
    uint32_t i = digitno + 1;
    while (true) {
      digits[i]++;
      if (digits[i] == 10) {
        diff -= 9;
        digits[i] = 0;
        i++;
      }
      else {
        break;
      }
    }
  }
  digitsum += (uint32_t)diff;
}

inline void loadDigits(digit_t* digits, uint32_t number)
{
  int rest = number;
  for (uint32_t i = 0; i < 10; i++) {
    digits[i] = rest % 10;
    rest /= 10;
  }
}

inline void sumDigits(digit_t* digits, uint32_t& sum)
{
  digit_t sumt = digits[0];
  sumt += digits[1];
  sumt += digits[2];
  sumt += digits[3];
  sumt += digits[4];
  sumt += digits[5];
  sumt += digits[6];
  sumt += digits[7];
  sumt += digits[8];
  sumt += digits[9];
  sum = (uint32_t)sumt;
}

bool* prime;
uint64_t * M;
uint32_t * digitsums;
uint64_t digitsums_M;

uint32_t countHarshadPrimes10(uint32_t start, uint32_t end)
{
  uint32_t count = 0;
  uint32_t digitsum = 0;
  digit_t digits[10];

  //head:
  uint32_t i = start;
  loadDigits(digits, start);
  sumDigits(digits, digitsum);

  uint32_t headEnd = start + (10 - start % 10);
  for (; i < headEnd; i++) {
    check(digitsum++, count, prime, i, M);
  }

  //body:
  loadDigits(digits, i);
  sumDigits(digits, digitsum);
  uint32_t bodyEnd = (end / 10) * 10;
  for (; i < bodyEnd; i += 10)
  {
    check(digitsum, count, prime, i, M);
    if (digitsum % 2 != 0) {
      check(digitsum + 2, count, prime, i + 2, M);
      check(digitsum + 4, count, prime, i + 4, M);
      check(digitsum + 6, count, prime, i + 6, M);
      check(digitsum + 8, count, prime, i + 8, M);
    }
    else {
      check(digitsum + 1, count, prime, i + 1, M);
      check(digitsum + 3, count, prime, i + 3, M);
      check(digitsum + 5, count, prime, i + 5, M);
      check(digitsum + 7, count, prime, i + 7, M);
      check(digitsum + 9, count, prime, i + 9, M);
    }
    addToDigit(digits, digitsum, 1, 1);
  }

  //tail:
  while (i < end) {
    check(digitsum++, count, prime, i++, M);
  }

  return count;
}

inline uint32_t getDigitSum(uint32_t num, const uint32_t * precalc, const uint64_t precalc_M, const uint32_t count)
{
  uint32_t upper = fastmod::fastdiv_u32(num, precalc_M);
  uint32_t lower = num - (upper * count);
  return precalc[lower] + precalc[upper];
}

uint32_t countHarshadPrimes100(uint32_t start, uint32_t end)
{
  uint32_t count = 0;
  uint32_t digitsum = 0;
  digit_t digits[10];

  //head:
  uint32_t i = std::min(start + (100 - start % 100), end);
  count = countHarshadPrimes10(start, i);

  //body:
  uint32_t bodyEnd = (end / 100) * 100;
  loadDigits(digits, i - 100);
  sumDigits(digits, digitsum);
  for (; i < bodyEnd; i += 100)
  {
    addToDigit(digits, digitsum, 2, 1);
    uint32_t mod = fastmod::fastmod_u32(digitsum, M[6], 6);
    switch (mod)
    {
    case 0:
      check100rest0(digitsum, count, prime, i, M);
      break;
    case 1:
      check100rest1(digitsum, count, prime, i, M);
      break;
    case 2:
      check100rest2(digitsum, count, prime, i, M);
      break;
    case 3:
      check100rest3(digitsum, count, prime, i, M);
      break;
    case 4:
      check100rest4(digitsum, count, prime, i, M);
      break;
    case 5:
      check100rest5(digitsum, count, prime, i, M);
      break;
    }
  }

  //tail:
  if (i < end) count += countHarshadPrimes10(i, end);

  return count;
}

uint32_t countHarshadPrimes100Precalc(uint32_t start, uint32_t end)
{
  uint32_t count = 0;
  uint32_t digitsum = 0;

  const uint32_t* precalc = digitsums;
  const uint64_t precalc_M = digitsums_M;
  const uint32_t precalc_count = 100000;

  //head:
  uint32_t i = std::min(start + (100 - start % 100), end);
  count = countHarshadPrimes10(start, i);

  //body:
  uint32_t bodyEnd = (end / 100) * 100;
  for (; i < bodyEnd; i += 100)
  {
    digitsum = getDigitSum(i, precalc, precalc_M, precalc_count);
    uint32_t mod = fastmod::fastmod_u32(digitsum, M[6], 6);
    switch (mod)
    {
    case 0:
      check100rest0(digitsum, count, prime, i, M);
      break;
    case 1:
      check100rest1(digitsum, count, prime, i, M);
      break;
    case 2:
      check100rest2(digitsum, count, prime, i, M);
      break;
    case 3:
      check100rest3(digitsum, count, prime, i, M);
      break;
    case 4:
      check100rest4(digitsum, count, prime, i, M);
      break;
    case 5:
      check100rest5(digitsum, count, prime, i, M);
      break;
    }
  }

  //tail:
  if (i < end) count += countHarshadPrimes10(i, end);

  return count;
}

int countHarshadPrimesMultithreaded(uint32_t start, uint32_t end, uint32_t func(uint32_t,uint32_t), uint32_t numthreads)
{
  if (numthreads == 0) {
    unsigned concurentThreadsSupported = std::thread::hardware_concurrency();
    numthreads = std::max(concurentThreadsSupported, 1U);
  }

  if (numthreads == 1) return func(start, end);

  uint32_t partitionsize = (end + numthreads - start) / numthreads;
  std::thread* threads = new std::thread[numthreads];
  uint32_t* sums = new uint32_t[numthreads];

  for (uint32_t threadNo = 0; threadNo < numthreads; threadNo++)
  {
    uint32_t threadBegin = start;
    uint32_t threadEnd = std::min(start + partitionsize, end);
    uint32_t currentThread = threadNo;
    threads[threadNo] = std::thread(
      [threadBegin, threadEnd, currentThread, &sums, func] {
        sums[currentThread] = func(threadBegin, threadEnd);
      }
    );
    start += partitionsize;
  }

  uint32_t sum = 0;
  for (uint32_t threadNo = 0; threadNo < numthreads; threadNo++)
  {
    threads[threadNo].join();
    sum += sums[threadNo];
  }

  delete[] threads;
  delete[] sums;
  return sum;
}

void precalcDigitSums(uint32_t count = 100000)
{
  digitsums = new uint32_t[count];

  uint32_t digitsum = 0;
  digit_t digits[10];
  loadDigits(digits, 0);

  for (uint32_t i = 0; i < count; i++)
  {
    digitsums[i] = digitsum;
    addToDigit(digits, digitsum, 0, 1);
  }

  digitsums_M = fastmod::computeM_u32(count);
  std::cout << "Precalculated " << count << " digitsums" << std::endl;
}


int main(int argc, char** argv)
{
  cxxopts::Options options("knowit23", "Count number of Harshad numbers whose digit sums are primes");
  options.add_options()
    ("h,help", "Print this help message")
    ("s,small", "Run small set (2 to 98765432)")
    ("p,precalc", "Precalc array of digitsums and use in calculation")
    ("t,threads", "Thread count (0 to auto detect)", cxxopts::value<std::uint32_t>()->default_value("0"))
    ("c,count", "Run count", cxxopts::value<std::uint32_t>()->default_value("10"))
    ;

  try {
    auto opt = options.parse(argc, argv);

    if (opt.count("help") > 0) {
      std::cout << options.help();
      return 0;
    }

    int count = opt["count"].as<uint32_t>();
    bool small = opt["small"].as<bool>();
    int numthreads = opt["threads"].as<uint32_t>();

    uint32_t end = small ? 98765432 : INT_MAX - 256;

    prime = buildPrimeSieve(90, false);
    M = new uint64_t[91];
    uint32_t (*func)(uint32_t, uint32_t) = countHarshadPrimes100;
    if (opt["precalc"].as<bool>()) {
      precalcDigitSums();
      func = countHarshadPrimes100Precalc;
    }
    for (uint32_t d = 1; d < 91; d++) M[d] = fastmod::computeM_u32(d);

    std::cout << "From 2 to " << end << std::endl;
    int result = 0;
    for (int i = 0; i < count; i++)
      result = countHarshadPrimesMultithreaded(2, end, func, numthreads);
    std::cout << result << std::endl;
  }
  catch (const cxxopts::OptionException & e)
  {
    std::cout << "error parsing options: " << e.what() << std::endl;
    exit(1);
  }
}

Et lite screenshot fra codegen-verktøyet mitt:

image

@bobcat4
Copy link

bobcat4 commented Jan 2, 2020

@terjew

-s -t 1 -c 100

          3,466.97 msec task-clock                #    1.000 CPUs utilized
                 3      context-switches          #    0.001 K/sec
                 0      cpu-migrations            #    0.000 K/sec
               140      page-faults               #    0.040 K/sec
    13,495,975,038      cycles                    #    3.893 GHz
    27,387,408,960      instructions              #    2.03  insn per cycle
     6,710,139,681      branches                  # 1935.450 M/sec
       135,834,587      branch-misses             #    2.02% of all branches

       3.467150855 seconds time elapsed

       3.467190000 seconds user
       0.000000000 seconds sys

Samme jobb med @michaelo sin siste
(countHarshadPrimesOptimized(2, 98765432) i en loop fra 0-99)

          4,965.99 msec task-clock                #    1.000 CPUs utilized
                63      context-switches          #    0.013 K/sec
                 0      cpu-migrations            #    0.000 K/sec
               136      page-faults               #    0.027 K/sec
    19,330,579,736      cycles                    #    3.893 GHz
    49,240,330,840      instructions              #    2.55  insn per cycle
     7,902,949,202      branches                  # 1591.416 M/sec
        13,909,003      branch-misses             #    0.18% of all branches

       4.966442363 seconds time elapsed

       4.964241000 seconds user
       0.002000000 seconds sys

Litt tregere, men mange færre branch-misses.

@terjew
Copy link

terjew commented Jan 2, 2020

@bobcat4

Interessante tall. Du kan gjerne også sammenlikne med min med --precalc I argumentlisten.

@bobcat4
Copy link

bobcat4 commented Jan 2, 2020

terjew

Interessante tall. Du kan gjerne også sammenlikne med min med --precalc I argumentlisten.

Det glemte jeg å si ja. 😄

--precalc ga noe dårligere resultat. Ikke mye men målbart (mindre enn 2%)

@terjew
Copy link

terjew commented Jan 2, 2020

@bobcat4
Skjønner, men den er egentlig mer relevant å sammenlikne med @michaelo sin variant mtp branching. Jeg vil nemlig gjette at det meste av misprediction i den raskeste løsningen skyldes loopingen i addToDigit.

@terjew
Copy link

terjew commented Jan 3, 2020

Jeg klarer fortsatt ikke helt gi slipp på denne, og trengte litt motivasjon, så jeg laget følgende:

image

EDIT: Ser ut som oppløsningen ble strupet så jeg lastet det opp her i tillegg

@michaelo
Copy link

michaelo commented Jan 3, 2020

@terjew

Jeg klarer fortsatt ikke helt gi slipp på denne

Hehe. Kjenner følelsen!

så jeg laget følgende:

Den var fin! 👍

Jeg nevnte forøvrig i en tidligere post at jeg hadde noen spor på å utlede rekka av summer direkte, men uansett hvor lur og/eller enkel jeg har forsøkt å være så ender det opp i for mange løkker til å bli ordentlig effektivt (det ble fort noe ala kjøretidskompleksiteten til posten til @DarioSucic lenger opp her) selv om jeg jobbet direkte på summen og rullet ut det dypeste nivået. Fikk den vel ned mot 300ms st, så ikke noe å hoppe i taket av.

@terjew
Copy link

terjew commented Jan 3, 2020

@michaelo

Jeg har vært inne på tanken selv. Rekka i seg selv er jo egentlig ganske rett frem, men 10 nøstede for-løkker er jo litt voldsomt. Kunne kanskje rullet ut 2 og 2 siffer, og nøste 5 slike rekker?

Min strategi er jo egentlig å droppe løkker og loop-variabler mest mulig slik at man kan slippe branchingen det fører til. Unntaket er addToDigit som jeg ikke fikk mer effektiv ved å unrolle, men denne kalles uansett bare 1 gang for hvert 100 tall.

@terjew
Copy link

terjew commented Jan 3, 2020

@michaelo
Hmm, fikk akkurat en ny idé... Hvis jeg kan splitte arbeidet i litt penere aligned chunker før trådene setter i gang blir jo egentlig alt mye lettere. Da kan jeg kanskje bake inn noen antakelser om modulo 6 også... Fristende! Men det må bli en annen dag.

@terjew
Copy link

terjew commented Jan 5, 2020

@michaelo @bobcat4 @DarioSucic

Jeg har tenkt litt... i rettferdighetens navn må man kunne argumentere for at når man sammenlikner benchmarks kjørt på et skjermkort til 9000 kroner så må det være rimelig å sammenlikne dem med tilsvarende benchmarks kjørt på en prosessor til samme pris.

Så, jeg fant en kollega med en 3950X og tigget meg til litt kjøretid på hans private boks. Her er resultatet:

Lite, ST: 34,5 ms
Lite, MT (precalc): 3,78ms
Lite, MT: 3,74ms

Stort, ST: 712ms
Stort. MT (precalc): 57ms
Stort. MT: 57ms

@michaelo
Copy link

michaelo commented Jan 6, 2020

@terjew cc: @DarioSucic @bobcat4

Hehe. Bra bra. Gyldig argument det!

Så er da de neste spørsmålene:

  • Hvor mye kan CUDA/shader-koden optimaliseres?
  • Varer også kodekalenderjula helt til påske?

?

@DarioSucic
Copy link

Fikk ned CUDA fra 5.4ms til 4.1ms ved å bytte fra Atomic Add til Parallell blokk-reduksjon slik det er beskrevet her.

Tror det er mye mer å hente, men det er i alle fall en start. 😅

#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>
#include <math.h>

#define u32 uint32_t
#define u8  uint8_t

#define f32 float
#define f64 double

#define BLOCK_SIZE 128

__device__ void warp_reduce(volatile u32 *block_data, u32 tid)
{
    block_data[tid] += block_data[tid + 32];
    block_data[tid] += block_data[tid + 16];
    block_data[tid] += block_data[tid +  8];
    block_data[tid] += block_data[tid +  4];
    block_data[tid] += block_data[tid +  2];
    block_data[tid] += block_data[tid +  1];
}

__device__ u32 compute_single(u32 k, u8 *prime)
{
    const u32 pow_table[9] = {1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000};
    u32 digit_sum = 0;
    for (u32 i = 0; i < 9; i++)
        digit_sum += (k / pow_table[i]) % 10;
    
    return (k % digit_sum == 0) & prime[digit_sum];
}

__global__ void kernel(u32 *block_counts, u8 *prime, u32 n)
{
    extern __shared__ u32 block_data[];

    u32 tid = threadIdx.x;
    u32 idx = blockIdx.x * blockDim.x + tid;
    if (idx > n) block_data[tid] = 0;
    else         block_data[tid] = compute_single(idx + 1, prime);
    __syncthreads();

    for(u32 s = blockDim.x / 2; s > 32; s >>= 1)
    {
        if(tid < s) block_data[tid] += block_data[tid + s];
        __syncthreads();
    }
    
    if(tid < 32) warp_reduce(block_data, tid);
    if(tid == 0) block_counts[blockIdx.x] = block_data[0];
}

u8 *sieve(u32 n)
{
    u32 lim = n + 1;
    u8 *prime = (u8 *)malloc(lim);
    memset(prime, 1, lim);
    prime[0] = 0;
    prime[1] = 0;

    for (u32 k = 2; k < lim; k++)
    {
        if (!prime[k])
            continue;
        for (u32 i = k * k; i < lim; i += k)
            prime[i] = 0;
    }

    return prime;
}

int main()
{
    const u32 n = 98765432;
    const u32 prime_limit = 9 * ceil(log10(n));

    cudaEvent_t st, et;
    cudaEventCreate(&st);
    cudaEventCreate(&et);

    u8 *prime = sieve(prime_limit);
    u8 *prime_device;
    cudaMalloc(&prime_device, prime_limit);
    cudaMemcpy(prime_device, prime, prime_limit, cudaMemcpyHostToDevice);


    u32 num_blocks = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;
    u32 shared_mem_size = BLOCK_SIZE * sizeof(u32);
    

    u32 *counts = (u32 *) malloc(num_blocks * sizeof(u32));
    u32 *counts_device;
    cudaMalloc(&counts_device, num_blocks * sizeof(u32));
    
    cudaEventRecord(st);
    kernel<<<num_blocks, BLOCK_SIZE, shared_mem_size>>>(counts_device, prime_device, n);
    cudaEventRecord(et);
    cudaEventSynchronize(et);
    float ms = 0;
    cudaEventElapsedTime(&ms, st, et);
    printf("Kernel execution: %.3f ms\n", ms);


    cudaMemcpy((void *)counts, counts_device, num_blocks * sizeof(u32), cudaMemcpyDeviceToHost);

    u32 total_count = 0;
    for(u32 i = 0; i < num_blocks; i++)
        total_count += counts[i];

    printf("Count: %d\n", total_count);


    cudaFree(counts_device);
    cudaFree(prime_device);
    free(prime);

    return 0;
}

@terjew
Copy link

terjew commented Jan 6, 2020

Fikk ned CUDA fra 5.4ms til 4.1ms ved å bytte fra Atomic Add til Parallell blokk-reduksjon slik det er beskrevet her.

Tror det er mye mer å hente, men det er i alle fall en start. 😅

Nais! 👏

Tar du en benchmark på stort sett også?

@terjew
Copy link

terjew commented Jan 6, 2020

@michaelo

  • Varer også kodekalenderjula helt til påske?

Minst! Og da må den vel avløses av noen kodepåskenøtter?

@DarioSucic
Copy link

DarioSucic commented Jan 8, 2020

Tar du en benchmark på stort sett også?

Endte opp med 85ms på det store datasettet, men egentlig jukser jeg litt siden mesteparten av summeringen nå er flyttet over til CPU utenfor kernelen. Om jeg får tid kan jeg se om jeg ikke får slengt inn litt bedre timing slik at resultatet er sammenlignbart med resten av tidene deres.

@terjew
Copy link

terjew commented Jan 16, 2020

Tar du en benchmark på stort sett også?

Endte opp med 85ms på det store datasettet, men egentlig jukser jeg litt siden mesteparten av summeringen nå er flyttet over til CPU utenfor kernelen. Om jeg får tid kan jeg se om jeg ikke får slengt inn litt bedre timing slik at resultatet er sammenlignbart med resten av tidene deres.

Tenkte på det samme når jeg så over koden. Du gjemmer bort et par ting i benchmarken, deriblant den siste summeringen men også memcpy til og fra CUDA-devicen. Tenker kanskje det hadde vært riktigere å la programmet loope f.eks 1000 ganger på lite datasett, og måle total tid brukt på dette inkludert de punktene, men uten å ta med opprettelse og teardown av CUDA-konteksten. Det utgjør vel ikke all verden, men når vi er nede på noen få millisekunder for hele kjøringen begynner alle sånne småting å utgjøre en relevant del av tiden.

@terjew
Copy link

terjew commented Jan 16, 2020

Evt time hele programmet med 1000 interne kjøringer fra shellet slik vi gjør for de andre variantene. Da får du jo med de 200 ms på opprettelsen av CUDA-kontekst etc, men dette fordeles jo utover de 1000 kjøringene så det skal vel ikke utgjøre mer enn 0,2ms per kjøring.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment