Skip to content

Instantly share code, notes, and snippets.

@knowitkodekalender
Created Dec 22, 2019
Embed
What would you like to do?

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?

@teilin

This comment has been minimized.

Copy link

@teilin teilin commented Dec 23, 2019

Go

import (
	"fmt"
	"math/big"
)

func sumIntDigits(num int64) int64 {
	if num == 0 {
		return 0
	}
	var slice []int64
	for {
		slice = append(slice, num%10)
		num /= 10
		if num <= 0 {
			break
		}
	}
	var sum int64
	for _, v := range slice {
		sum += v
	}
	return sum
}

func main() {
	count := 0
	for i := 1; i <= 98765432; i++ {
		sumDigits := sumIntDigits(int64(i))
		isHarshadNumber := int64(i)%sumDigits == 0
		if isHarshadNumber {
			if big.NewInt(sumDigits).ProbablyPrime(0) {
				count += 1
			}
		}
	}
	fmt.Println(count)
}
@DarioSucic

This comment has been minimized.

Copy link

@DarioSucic DarioSucic commented Dec 23, 2019

Forferdelig 32s løsning i Python woohoo

import numpy as np

def sieve(n):
    lim = n + 1
    not_prime = np.zeros(lim, dtype=np.uint8)
    prime = set()
    for k in range(2, lim):
        if not_prime[k]: continue
        not_prime[k*k:lim:k] = 1
        prime.add(k)
    return prime

ks = np.arange(1, 98765432+1, dtype=np.uint32)

tmp = np.copy(ks)
digit_sums = np.zeros_like(ks)
for _ in range(len(str(98765432)) + 1):
    digit_sums += tmp % 10
    tmp //= 10

prime = sieve(98765432)

print(sum(1 for ds in digit_sums[ks % digit_sums == 0] if ds in prime))
@Nilzone-

This comment has been minimized.

Copy link

@Nilzone- Nilzone- commented Dec 23, 2019

~2s..

const isPrime = num => {
	for (let j = 2, k = Math.sqrt(num); j <= k; j++)
		if(num % j === 0) return false; 
	return num > 1;
}

const N = 98765432;
let [i, _i, s, c] = [0, 0, 0, 0];
for (i = 1; i <= N; i++) {
	[s, _i] = [0, i];
	while (_i) {
		s += _i%10;
		_i = (_i/10) | 0;
	}
	c += i%s === 0 && isPrime(s) ? 1 : 0
}
console.log(c)
@finnmich

This comment has been minimized.

Copy link

@finnmich finnmich commented Dec 23, 2019

C# løsning, kjører på litt over et sekund med AsParallel for litt latmanns multithreading :)
Unngår også strenger i delsum-utregningen

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

namespace Kodekalender
{
    public class Luke23
    {
        public static long Run()
        {
            return Enumerable.Range(1, 98765431).AsParallel().Where(n => IsHarshadPrime(n)).Count();
        }

        public static bool IsHarshadPrime(int i)
        {
            var partsSum = PartSum(i);
            return (i % partsSum == 0 && IsPrime(partsSum));
        }

        private static int PartSum(int i)
        {
            var partSum = 0;
            while (i > 0)
            {
                partSum += i % 10;
                i /= 10;
            }
            return partSum;
        }


        public static bool IsPrime(int number)
        {
            if (number <= 1) return false;
            if (number == 2) return true;
            if (number % 2 == 0) return false;

            var boundary = (int)Math.Floor(Math.Sqrt(number));

            for (int i = 3; i <= boundary; i += 2)
                if (number % i == 0)
                    return false;
            return true;
        }
    }
}
@matslindh

This comment has been minimized.

Copy link

@matslindh matslindh commented Dec 23, 2019

Rett fram python-løsning med numerisk tallsum + parallelisering fordi det er en slik dag

import math
from multiprocessing import Pool

n = 9*9
prime_table = list(range(0, n + 1))
mm = int(math.sqrt(n))

for x in range(3, mm, 2):
    for y in range(x*2, n, x):
        prime_table[y] = 0


def is_prime(n):
    if n == 2:
        return True

    if n % 2 == 0 or n < 2:
        return False

    return bool(prime_table[n])


def sum_digits(n):
    s = 0
    
    while n:
        s += n % 10
        n //= 10
        
    s += n
    
    return s


def is_harshad(n):
    d = sum_digits(n)

    if n % d == 0:
        return is_prime(d)
        
    return False
    

def test_is_harshad():
    assert is_harshad(1729)
    assert not is_harshad(1730)


def is_harshad_range(interval):
    s = 0

    for i in range(interval[0], interval[1]):
        if is_harshad(i):
            s += 1

    return s


if __name__ == '__main__':
    p = Pool(8)
    sums = p.map(is_harshad_range, [
        (1, 12345679),
        (12345679 * 1, 12345679 * 2),
        (12345679 * 2, 12345679 * 3),
        (12345679 * 3, 12345679 * 4),
        (12345679 * 4, 12345679 * 5),
        (12345679 * 5, 12345679 * 6),
        (12345679 * 6, 12345679 * 7),
        (12345679 * 7, 12345679 * 8+1),
    ])
    
    print(sum(sums))
@lassem

This comment has been minimized.

Copy link

@lassem lassem commented Dec 23, 2019

Kjøretid 1m52s... Mao lite optimalisering og kun 1 tråd.

#!/usr/bin/env python3

def isPrime(n):
    if n <= 3:
        return n > 1
    elif n % 2 == 0 or n % 3 == 0:
        return False
    i = 5
    while i * i <= n:
        if n % i == 0 or n % (i + 2) == 0:
            return False
        i = i + 6
    return True

def sumDigits(n):
    s = 0
    while n:
        s += n % 10
        n //= 10
    return s

def isHarshad(n):
    s = sumDigits(n)
    m = n % s
    return (m == 0) and isPrime(s)

r = 0
for n in range(1, 98765433):
    if isHarshad(n):
        r += 1

print(r)
@jonathangjertsen

This comment has been minimized.

Copy link

@jonathangjertsen jonathangjertsen commented Dec 23, 2019

Pythonløsning som tok "rakk å lage kaffe og rydde opp på kjøkkenet" sekunder

primes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, }
ans = 0
for num in range(1, 98765432+1):
    digit_sum = sum(int(d) for d in list(str(num)))
    if digit_sum in primes and num % digit_sum == 0:
        ans += 1
print(ans)
@gunnar2k

This comment has been minimized.

Copy link

@gunnar2k gunnar2k commented Dec 23, 2019

Elixir ❤️

is_prime? = fn n ->
  cond do
    n == 2 || n == 3 ->
      true

    true ->
      floored_sqrt =
        :math.sqrt(n)
        |> Float.floor
        |> round

      !Enum.any?(2..floored_sqrt, &(rem(n, &1) == 0))
  end
end

is_harshadprime? = fn n ->
  sum =  Integer.digits(n) |> Enum.sum()

  Integer.mod(n, sum) == 0 && is_prime?.(sum)
end

1..98_765_432
|> Enum.filter(&is_harshadprime?.(&1))
|> Enum.count()
|> IO.puts()
@Aha43

This comment has been minimized.

Copy link

@Aha43 Aha43 commented Dec 23, 2019

C#

class Luke23 : KnowItJuleKalenderLuke
{
    public Task<string> OpenAsync()
    {
        var result = 0;
        var primes = new Primes(100);
        for (long i = 2; i <= 98765432; i++)
        {
            var sum = i.Digits().Sum();
            if (primes.IsPrime(sum) && (i % sum == 0)) result++;
        }

        return Task.FromResult(result.ToString());
    }
}

class Primes
{
    private bool[] _primes;

    public long N => _primes.Length;

    public Primes(long N)
    {
        _primes = new bool[N];
        for (var i = 2; i < N; i++) _primes[i] = true;
        long lastPrime = 2;
        while (lastPrime != -1)
        {
            MarkNonPrimes(lastPrime);
            lastPrime = NextPrime(lastPrime);
        }
    }

    private long NextPrime(long lastPrime)
    {
        for (long nextPrime = lastPrime + 1; nextPrime < _primes.Length; nextPrime++)
        {
            if (_primes[nextPrime]) return nextPrime;
        }

        return -1;
    }

    private void MarkNonPrimes(long prime)
    {
        for (long i = prime + prime; i < _primes.Length; i += prime)
        {
            _primes[i] = false;
        }
    }

    public bool IsPrime(long n)
    {
        return n < 0 ? false : _primes[n];
    }

}
@Bullhill

This comment has been minimized.

Copy link

@Bullhill Bullhill commented Dec 23, 2019

Dagens må kunne sies å ha vært plankekjøring

$limit = 98765432;
$c = 0;

for($i=1;$i<=strlen($limit)*9;$i++)
{
	if(is_prime($i))
		$primenums[] = $i;
}

for($i = 1; $i <= $limit; $i++)
{
	$nums = str_split($i);
	$sum = 0;
	foreach($nums as $num)
	{
		$sum += $num;
	}
	
	if(in_array($sum, $primenums) && $i % $sum == 0)
		$c++;
	
	if($i % 500000 == 0)
		echo "\r" . $i ;
}
echo "\rSvar: " . $c . "\n";

function is_prime($number)
{
    if ( $number == 1 ) {
        return false;
    }
    if ( $number == 2 ) {
        return true;
    }
    $x = sqrt($number);
    $x = floor($x);
    for ( $i = 2 ; $i <= $x ; ++$i ) {
        if ( $number % $i == 0 ) {
            break;
        }
    }
 
    if( $x == $i-1 ) {
        return true;
    } else {
        return false;
    }
}
@Suppen

This comment has been minimized.

Copy link

@Suppen Suppen commented Dec 23, 2019

C++. Kjøretid på litt over 19s, selv med Sieve of Eratosthenes

#include <iostream>

bool* eratosthenes(int n) {
        bool* sieve = new bool[n];
        sieve[0] = false;
        sieve[1] = false;
        for (int i = 2; i <= n; i++) {
                sieve[i] = true;
        }

        for (int i = 2; i < n; i++) {
                if (!sieve[n]) {
                        continue;
                }

                int prod = i+i;
                while (prod < n) {
                        sieve[prod] = false;
                        prod += i;
                }
        }

        return sieve;
}

int digit_sum(int n) {
        int sum = 0;
        while (n > 0) {
                sum += n % 10;
                n = n / 10;
        }

        return sum;
}

bool is_harshad(bool* primes, int n) {
        int ds = digit_sum(n);
        return primes[ds] && n % ds == 0;
}

int main(int argc, char** argv) {
        int limit = 98765432;
        bool* primes = eratosthenes(limit+1);

        int count = 0;
        for (int n = 0; n <= limit; n++) {
                if (is_harshad(primes, n)) {
                        count++;
                }
        }

        std::cout << count << std::endl;
}

EDIT: Prøvde å bytte ut Sieve of Eratosthenes med @bobcat4 sin is_prime. Kjøretiden ble da redusert til litt over 1.4 sekunder, men jeg skal ikke si at jeg forstår hvordan is_prime funker. Et eller annet fuckery med at et primtall (med unntak av 2 og 3) alltid er ±1 fra et multiplum av 6?

EDIT2: Ved værmere analyse av is_prime ser den ut til å være en variant av primtallsgeneratoren som brukes i Haskellkoden under.

EDIT3: Etter å la hest litt kommentarer her innser jeg jo at å generere alle primtall opp til 98765432 er idiotisk, da siffersummen maksimalt er 71... Det er der mesteparten av tidsinnsparingen kommer fra ved å bruke is_prime

@Suppen

This comment has been minimized.

Copy link

@Suppen Suppen commented Dec 23, 2019

Haskell med kjøretid ca. 23.4s.

import qualified Data.Set as S

-- Note: Digits are in reverse order
type Digits = [Int]

limit :: Int
limit = 98765432

-- Stolen from https://www.reddit.com/r/haskell/comments/35vc31/the_real_way_to_generate_a_list_of_primes_in/
primes :: S.Set Int
primes = S.fromAscList . takeWhile (<maxPrime) $ 2 : 3 : 5 : primes'
    where isPrime (p:ps) n = p*p > n || n `rem` p /= 0 && isPrime ps n
          primes' = 7 : filter (isPrime primes') (scanl (+) 11 $ cycle [2,4,2,4,6,2,6,4])
          maxPrime = (*9) . length . toDigits $ limit

toDigits :: Int -> Digits
toDigits 0 = []
toDigits n = n `mod` 10 : toDigits (n `div` 10)

digitSum :: Int -> Int
digitSum = sum . toDigits

isHarshad :: Int -> Bool
isHarshad n = let ds = digitSum n
              in  n `mod` ds == 0 && S.member ds primes

main =
    print
        . length
        . filter isHarshad
        $ [1..limit]
@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 23, 2019

Denne var nesten for enkel. Men kanskje et (nytt) objekt for optimalisering? ;-)

Edit: Her er det ikke noe CPU-cache-optimalisering å hente, så optimaliseringen må finnes matematisk. Interessant!

Kjøretid "out of the box": 1s (uoptimalisert av mennesker).

C-kode:

#include <stdio.h>
#include <string.h>
#include <stdlib.h>

int is_prime(int n)
{
        if (n == 1) return 0;
        if (n < 4) return 1;
        if (n % 2 == 0) return 0;
        if (n % 3 == 0) return 0;

        int i = 5;
        int w = 2;

        while (i * i <= n) {
                if (n % i == 0) return 0;
                i += w;
                w = 6 - w;
        }
        return 1;
}

int digitSum (int n)
{
        int sum = 0;
        while (n != 0)
        {
                int rem = n % 10;
                sum += rem;
                n /= 10;
        }
        return sum;
}

int main ()
{
        int cnt = 0;
        for (int i = 1 ; i < 98765433 ; i++) {
                int dsum = digitSum(i);
                if ((i % dsum == 0) && is_prime(dsum)) cnt++;
        }
        printf("%d\n", cnt);
        return 0;
}
@jostein-skaar

This comment has been minimized.

Copy link

@jostein-skaar jostein-skaar commented Dec 23, 2019

primes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, }

Ah, det hadde jeg ikke tenkt på. Og jeg som tenkte at nå må jeg jo i hvert fall finne en smart måte å finne primtall på...

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 23, 2019

@jostyposty

Spørs om denne er raskere:
primes = { 0,0,1,1,0,1,0,1,0,0,0,1,0,1,0,0,0,1,0,1,0,0,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,0,0,1,0}

I hvert fall i teorien. "Hardware loves arrays". Og her blir det en array-lookup. primes = { 2, 3, 5, 7, 11, 13,... krever traversering.

Edit: Ingen målbar forskjell på min boks.

@Aha43

This comment has been minimized.

Copy link

@Aha43 Aha43 commented Dec 23, 2019

@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 23, 2019

Rett-frem uskviset løsning i C.

Single thread:

  • Time (mean ± σ): 1.003 s ± 0.007 s

16 thread

  • Time (mean ± σ): 228.6 ms ± 3.2 ms
#define MULTITHREADED

#include <stdio.h>
#include <stdint.h>
#include <stdbool.h>
#include <math.h>
#include <string.h>
#include <assert.h>
#include "../timer.h"

#ifdef MULTITHREADED
#include <pthread.h>
typedef struct
{
  pthread_t thread_ref;
  bool* prime_map;
  uint32_t from;
  uint32_t to;
  int result;
} thread_definition_t;
static const int NUM_THREADS = 16;
#endif

static const int MAX_PRIME = 100;

static inline void generate_primes_sieve(uint64_t until, bool num_map[]) {
    uint64_t max_square_root = sqrt(until)+0.5;
    bool eliminated[until];
    memset(eliminated, 0, sizeof(eliminated));
    
    for (uint64_t i = 2; i <= until; ++i)
    {
        if (!eliminated[i])
        {
            num_map[i] = true;
            if (i <= max_square_root)
            {
                for (uint64_t j = i * i; j <= until; j += i)
                {
                    // printf("j: %llu\n", j);
                    eliminated[j] = true;
                }
            }
        }
    }
}

static inline int sum_of_digits(uint64_t num) {
    int sum = 0;

    while(num > 0) {
        sum += num % 10;
        num /= 10;
    }

    return sum;
}


static inline int solve(uint64_t limit) {
    bool NUM_IS_PRIME_MAP[MAX_PRIME] = {0};
    generate_primes_sieve(MAX_PRIME, NUM_IS_PRIME_MAP);
    int result = 0;
    int sum;
    for (int i=0; i<limit+1; i++) {
        sum = sum_of_digits(i);
        if(NUM_IS_PRIME_MAP[sum] && i % sum == 0) {
            result++;
        }
    }

    return result;
}

#ifdef MULTITHREADED
void* mt_worker(void* arg) {
    thread_definition_t *self = arg;
    int result = 0;

    int sum;
    for (int i=self->from; i<self->to; i++) {
        sum = sum_of_digits(i);
        if(self->prime_map[sum] && i % sum == 0) {
            result++;
        }
    }


    self->result = result;
    return NULL;
}

static inline uint32_t min(uint32_t a, uint32_t b) {
    if(a < b) return a;
    return b;
} 

static inline int mt_solve(uint64_t limit) {
    bool NUM_IS_PRIME_MAP[MAX_PRIME] = {0};
    generate_primes_sieve(MAX_PRIME, NUM_IS_PRIME_MAP);
    thread_definition_t threads[NUM_THREADS];

    uint64_t chunk_size = ceil(1.0*limit/NUM_THREADS);
    for(int i=0; i<NUM_THREADS; i++) {
        threads[i].prime_map = NUM_IS_PRIME_MAP;
        threads[i].from = i*chunk_size;
        threads[i].to = min(limit, (i+1)*chunk_size);
        threads[i].result = 0;
        pthread_create(&threads[i].thread_ref, NULL, mt_worker, &threads[i]);
    }

    int result = 0;
    for(int i=0; i<NUM_THREADS; i++) {
        pthread_join(threads[i].thread_ref, NULL);
        result += threads[i].result;
    }

    return result;
}
#endif

int main() {
    // timer_t timer;
    // timer_start(&timer);
    int solution;
    #ifdef MULTITHREADED
    solution = mt_solve(98765432+1);
    #else
    solution = solve(98765432+1);
    #endif
    printf("Solution: %d\n", solution);
    // assert(solution == 724713);
    // printf("solve finished: %llu\n", timer_lap(&timer));
    return 0;
}
@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 23, 2019

@Bullhill

Dagens må kunne sies å ha vært plankekjøring

Helt klart. Spesielt om man har tatt de foregående dagene da dette da brått bare var en sammensetning+lett tweaking av av biter som man gjerne da har tatt frem.

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 23, 2019

@michaelo

Rett-frem uskviset løsning i C.

Single thread:

  • Time (mean ± σ): 1.003 s ± 0.007 s

16 thread

  • Time (mean ± σ): 228.6 ms ± 3.2 ms

Lukter "false sharing" av denne MT-en (aka dirty cache lines). Du bør kunne komme under 100 ms med 16 tråder.

@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 23, 2019

@bobcat4

Lukter "false sharing" av denne MT-en. Du bør kunne komme under 100 ms med 16 tråder.

Bør helt klart kunne komme lavere enn jeg er pr nå, men ikke så lavt som du postulerer. i7 7700HQ'en min har 4 cores / 8 threads. Har enda ikke fått pinpointet hvorfor den faktisk skreller av ~30 ms når jeg går fra 8 til 16 tråder heller.

@Bullhill

This comment has been minimized.

Copy link

@Bullhill Bullhill commented Dec 23, 2019

primes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, }

Ah, det hadde jeg ikke tenkt på. Og jeg som tenkte at nå må jeg jo i hvert fall finne en smart måte å finne primtall på...

Startet med en kjempeliste av primtall jeg fra luke 20 (største primtall var 1299827), helt til jeg kom på at største primtall tall man kom til å sjekke om var primtall var lengden av største tall * 9 som ble 72. Utgjorde mye på kjøretiden å ikke la den få for mange primtall å jobbe med
73 vil den aldri sjekke etter i denne oppgaven.

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 23, 2019

@Bullhill

Startet med en kjempeliste av primtall jeg fra luke 20 (største primtall var 1299827), helt til jeg kom på at største primtall tall man kom til å sjekke om var primtall var lengden av største tall * 9 som ble 72

Største primtall er 61 og det blir også siste primtall i lookup-tabellen.

@bobcat4

This comment has been minimized.

@Bullhill

This comment has been minimized.

Copy link

@Bullhill Bullhill commented Dec 23, 2019

@bobcat4
89999999 vil vell gi 71?

@stianeikeland

This comment has been minimized.

Copy link

@stianeikeland stianeikeland commented Dec 23, 2019

Parallell clojure utgave, med @jonathangjertsen's primtallstriks, som fungerer rimelig fiffig i clojure, da datastrukturene også er funksjoner.

Rimelig langt unna de parallelle C-løsningene over i hastighet dog...

(ns knowit2019.23a
  (:require [clojure.core.reducers :as r]))

(def prime? #{2, 3, 5, 7, 11, 13, 17, 19,
              23, 29, 31, 37, 41, 43, 47,
              53, 59, 61, 67, 71, 73})

(defn digit-sum [^long n]
  (loop [n n
         acc 0]
    (if (zero? n) acc
        (recur (quot n 10) (+ acc (mod n 10))))))

(defn harshad-prime? [^long n]
  (let [sum (digit-sum n)]
    (and (zero? (mod n sum))
         (prime? sum))))

(->> (range 1 98765433)
     (vec)
     (r/filter harshad-prime?)
     (r/map (fn [_] 1))
     (r/fold +))
=> 724713
@jonathangjertsen

This comment has been minimized.

Copy link

@jonathangjertsen jonathangjertsen commented Dec 23, 2019

@bobcat4

Spørs om denne er raskere:
primes = { 0,0,1,1,0,1,0,1,0,0,0,1,0,1,0,0,0,1,0,1,0,0,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,0,0,1,0}

I hvert fall i teorien. "Hardware loves arrays". Og her blir det en array-lookup. primes = { 2, 3, 5, 7, 11, 13,... krever traversering.

Edit: Ingen målbar forskjell på min boks.

I et Python-program som regner ut list(str(num)) for hvert eneste tall tror jeg det har fint lite å si altså :D Men det er jo et snedig triks, hvis du skal blodtrimme C-kode på denne så tipper jeg det er den raskeste måten å bestemme primtallitet på.

@MarcusTL12

This comment has been minimized.

Copy link

@MarcusTL12 MarcusTL12 commented Dec 23, 2019

Dette var kanskje letteste oppgaven så langt. Endte opp med å bli en oneliner i julia uten å prøve. Kjører på 20 sec ish.

count(n -> ((s = sum(digits(n))) |> isprime) && n % s == 0, 1 : 98765432)

Edit: Om det er en ting jeg har lært av å gjøre de her kalenderoppgavene, så må det være at all form for innebygd omgjøring til siffer er forferdelig tregt. Ingenting slår rett fram modulo og heltalsdivisjon der. Så med å skrive min egen crosssum så gikk kjøretiden ned til ca 3 sec.

Edti2: Med alt snakket her om at det å sjekke primtall tar tid så måtte jeg teste litt om det var mulig å gjøre noe raskere en den generelle isprime funksjonen i Primes.jl (som forøvrig er helt sykt rask) siden vi vet at tverrsummen er 71 eller mindre. Ved å lage et Set av de relevante primtallene og sjekke om de var der så gikk kjøretiden ned til 2.4 sec.

using Primes


function crosssum(n)
    acc = 0
    while n > 0
        acc += n % 10
        n ÷= 10
    end
    acc
end


function main()
    rel_primes = Set(primes(71))
    count(n -> ((s = crosssum(n)) in rel_primes) && n % s == 0, 1 : 98765432)
end
@Sjark

This comment has been minimized.

Copy link

@Sjark Sjark commented Dec 23, 2019

c#

static void Main(string[] args)
{
    var numberOfHarshadprimtall = 0;
    Enumerable.Range(1, 98765432).AsParallel()
        .ForAll(i =>
        {
            var sumOfDigits = i.ToString().Select(a => a - '0').Sum();

            if (i % sumOfDigits == 0 && IsPrime(sumOfDigits))
            {
                Interlocked.Increment(ref numberOfHarshadprimtall);
            }
        });

    Console.WriteLine(numberOfHarshadprimtall);
}

static bool IsPrime(int number)
{
    if (number <= 1)
    {
        return false;
    }

    if (number == 2)
    {
        return true;
    }

    if (number % 2 == 0)
    {
        return false;
    }

    var boundary = (int)Math.Floor(Math.Sqrt(number));

    for (int i = 3; i <= boundary; i += 2)
    {
        if (number % i == 0)
        {
            return false;
        }
    }

    return true;
}
@danieman

This comment has been minimized.

Copy link

@danieman danieman commented Dec 23, 2019

Deilig med en kjapp og enkel nå i julestria. 😄

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 23, 2019

@bobcat4
89999999 vil vell gi 71?

Det har du naturligvis rett i, @Bullhill (morsomt nick som jeg ikke skjønte før ganske nylig) :-)

Det gikk litt fort i svingene. Det jeg mente å si var at 69 er den største tverrsummen et Harshad-tall opp til og med 98765432 kan få. Ikke rare forbedringen, men det er jo litt. Jaja, bare ett primtall mindre men... :-)

@Runebjo

This comment has been minimized.

Copy link

@Runebjo Runebjo commented Dec 23, 2019

JavaScript:

function isPrime(number) {
    if (number === 1) return false;
    for (let i = 2; i < number; i++) {
        if (number % i === 0) return false;
    }
    return true;
}

function isHarshadPrime(number) {
    let harsh = number.toString().split('').map(n => parseInt(n)).reduce((a, b) => a + b);
    return number % harsh === 0 && isPrime(harsh);
}

let counter = 0;
for (let i = 1; i <= 98765432; i++) {
    if (isHarshadPrime(i)) {
        counter++;
    }
}
console.log(counter);
@hakonrossebo

This comment has been minimized.

Copy link

@hakonrossebo hakonrossebo commented Dec 23, 2019

F#, genererte først et primtall sett, men var jo mye greiere med noe ferdig

let sluttAntall = 98765432 
let primtall = [2; 3; 5; 7; 11; 13; 17; 19; 23; 29; 31; 37; 41; 43; 47; 53; 59; 61; 67; 71; 73] |> Set.ofList

let inline sumAvSiffer n =
    let mutable sum = 0
    let mutable i = n
    while i > 0 do
        sum <- sum + i % 10
        i <- i / 10
    sum
    
let inline erPrimtall n = Set.contains n primtall

let finnDagensSvar () =
    seq {
        for n in 1..sluttAntall do
            let sifferSum = sumAvSiffer n
            if n % sifferSum = 0 && erPrimtall sifferSum then
                yield n
    }
    |> Seq.length

finnDagensSvar()

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 23, 2019

Edit: Her er det ikke noe CPU-cache-optimalisering å hente, så optimaliseringen må finnes matematisk. Interessant!

Hvor er @terjew når vi trenger han? ;-)

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 23, 2019

@bobcat4 @michaelo

Edit: Her er det ikke noe CPU-cache-optimalisering å hente, så optimaliseringen må finnes matematisk. Interessant!

Hvor er @terjew når vi trenger han? ;-)

Var opptatt med å prøve å optimalisere C-koden din fra i går 😁

Har laget en optimalisert løsning (C#) basert på kode fra tidligere løsninger, men skal ut og julehandle og så rydde litt før jeg poster koden.

EDIT: Var helt ute og syklet når jeg kom frem til de tallene jeg postet her tidligere. Glem dem, jeg prøver på nytt senere.

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 23, 2019

@terjew, imponerende kjøretid. Litt rart at kjøretiden bare ble halvvert.

Det er kun primtall opp til 69 som brukes, så en hvilken som helst primtallgenerator duger her.

@aude

This comment has been minimized.

Copy link

@aude aude commented Dec 23, 2019

Go (1s)

package main

import (
	"fmt"
	"runtime"
)

// Hvor mange tall fra 1 til og med 98765432 er Harshadprimtall?
const min uint64 = 1
const max uint64 = 98765432

// largest possible prime is less than 9+9+9+9+9+9+9+9=72
const largestPrime = 72

func digitSum(n uint64) uint64 {
	sum := uint64(0)
	for ; n > 0; n /= 10 {
		sum += n % 10
	}
	return sum
}

// false == prime
var sieve [largestPrime + 1]bool

func fillSieve() {
	sieve[0] = true
	sieve[1] = true
	for i := 2; i < len(sieve); i++ {
		// if prime
		if sieve[i] == false {
			// fill all multiples of this prime (except i*1)
			for j := 2; i*j < len(sieve); j++ {
				sieve[i*j] = true
			}
		}
	}
}

func prime(n uint64) bool {
	return !sieve[n]
}

func harshadPrime(n uint64) bool {
	ds := digitSum(n)
	if n%ds == 0 && prime(ds) {
		return true
	}
	return false
}

func main() {
	count := uint64(0)

	fillSieve()

	// spawn workers, one per available thread
	workers := uint64(runtime.GOMAXPROCS(0))
	c := make(chan uint64, workers)
	split := (max - min) / workers
	for i := uint64(0); i < workers; i++ {
		start := min + i*split + i
		end := min + (i+1)*split + i
		if end > max {
			end = max
		}

		go func() {
			var partcount uint64 = 0
			for n := start; n <= end; n++ {
				if harshadPrime(n) {
					partcount++
				}
			}
			c <- partcount
		}()
	}

	// collect results from workers
	for i := uint64(0); i < workers; i++ {
		count += <-c
	}

	fmt.Println(count)
}
@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 23, 2019

Ingen som har prøvd å løse denne på et par tusen Cuda-cores i parallell? :-)

@christianfosli

This comment has been minimized.

Copy link

@christianfosli christianfosli commented Dec 23, 2019

Rust på 4-5 sekunder

use primes::is_prime;

fn is_harshad_prime(n: u32) -> bool {
    let sum_digits = sum_digits(n);
    let is_harshad = n % sum_digits == 0;
    is_harshad && is_prime(sum_digits as u64)
}

fn sum_digits(n: u32) -> u32 {
    let mut n = n;
    let mut sum = 0;
    for _ in 0..((n as f32).log10() as u32 + 1) {
        sum += n % 10;
        n = n / 10;
    }
    sum
}

fn main() {
    let mut count = 0;
    for i in 1..98_765_433 {
        if is_harshad_prime(i) {
            count += 1;
        }
    }
    println!("{}", count);
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_sum_digs() {
        assert_eq!(sum_digits(123), 6);
    }

    #[test]
    fn test_is_harshad_prime_should_be_true() {
        assert_eq!(is_harshad_prime(1729), true);
    }

    #[test]
    fn test_is_harshad_prime_should_be_false() {
        assert_eq!(is_harshad_prime(1730), false);
    }
}
@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 23, 2019

@terjew
Nydelig! Jeg har et par tiltak jeg skal fullføre når treet er pynta og kidsa er i seng.

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 23, 2019

@michaelo @bobcat4
Må dessverre skuffe, det var helt feil tall. Sånn går det når man har dårlig tid! Ser mer på det senere, men tallene jeg postet var off med en faktor på 100

@AndersSteenNilsen

This comment has been minimized.

Copy link

@AndersSteenNilsen AndersSteenNilsen commented Dec 23, 2019

Løste denn først med brute force. Deretter løste den ved å generere alle kombinasjoner av termer av alle primtall opp til 71.
Brute-force ~200sek, "Smart"-løsning ~100sek, så ikke helt optimaliseringen jeg hadde håpet på.
Det er 5994 forskjellige kombinasjoner/måter å addere opp til alle primtallene (gitt maks 8 tall og hvert tall er 1<=x<=9) men antall permisjoner (når rekkefølgen betyr noe) av disse kombinasjonene er 24 786 388...

from itertools import combinations, product, permutations
from sympy.utilities.iterables import multiset_permutations
from time import time

start_brute_dorce = time()
from math import factorial

composite_roof = 98765432
primes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71}


def isHarshadPrime(contestant):
    h = contestant
    digit_sum = 0
    while h > 0:
        digit_sum += h % 10
        h //= 10
    if digit_sum in primes:
        return contestant % digit_sum == 0
    return False


prime_terms = {p: set() for p in primes}


def generate_terms(target, terms, last_term=0, iteration=0):
    if sum(terms) == target:
        prime_terms[target].add(tuple(terms))
        return
    if sum(terms) > target:
        return
    if iteration >= 8:
        return False

    for i in range(last_term, 10):
        generate_terms(target, terms + [i], i, iteration + 1)


def concat_to_number(number_list):
    return_number = 0
    for number in number_list:
        return_number *= 10
        return_number += number
    return return_number


total_perms = 0


def calculate_permutations(numbers, p, number_length=8):
    global total_perms
    length = len(numbers)
    numbers_list = list(numbers) + [0] * (number_length - length)
    harshad_numbers = 0
    for perm in multiset_permutations(numbers_list):
        total_perms += 1
        harshad_numbers += concat_to_number(perm) % p == 0
    return harshad_numbers


def calcAllIsHarshad(lower, upper):
    harshads = 0
    for i in range(lower, upper):
        harshads += isHarshadPrime(i)
    return harshads


# brute_force_harshad_primes = calcAllIsHarshad(1, composite_roof)
# print(brute_force_harshad_primes, time() - start_brute_dorce)
start = time()

out_of_bound_harshad_primes = calcAllIsHarshad(composite_roof + 1, 100000000)
for i in range(composite_roof + 1, 100000000):
    out_of_bound_harshad_primes += isHarshadPrime(i)

for p in primes:
    generate_terms(p, [], 1)

number_of_harshad_primes = 0
for p in primes:
    print(p)
    for term in prime_terms[p]:
        number_of_harshad_primes += calculate_permutations(term, p)

print(number_of_harshad_primes - out_of_bound_harshad_primes, time() - start)
@DarioSucic

This comment has been minimized.

Copy link

@DarioSucic DarioSucic commented Dec 23, 2019

Ingen som har prøvd å løse denne på et par tusen Cuda-cores i parallell? :-)

👋

Nå vil jeg nok være den første til å innrømme at jeg egentlig ikke kan CUDA noe særlig, men jeg måler i alle fall en kjøretid på ~5.4ms for hele beregningen med CUDA events (ikke vist i koden). Det uheldige er at CUDA sin oppstart tar ~200ms, så kjøretiden til programmet blir jo rundt 200ms.

Dette er på en 1080 Ti gjennom en PCIE 3.0 x2 tilkobling (via USB-C / Thunderbolt 3).

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

#define u32 uint32_t
#define u8  uint8_t

__global__ void kernel(u32 *count, u8 *prime, u32 n)
{
    u32 idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx > n)
        return;

    u32 k = idx + 1;

    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;

    atomicAdd(count, (k % digit_sum == 0) & prime[digit_sum]);
}

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

    u8 *prime = sieve(prime_limit);
    u8 *device_prime;

    cudaMalloc(&device_prime, prime_limit);
    cudaMemcpy(device_prime, prime, prime_limit, cudaMemcpyHostToDevice);

    u32 count, *device_count;
    cudaMalloc(&device_count, sizeof(u32));
    kernel<<<(n + 255) / 256, 256>>>(device_count, device_prime, n);
    cudaMemcpy((void *)&count, device_count, sizeof(u32), cudaMemcpyDeviceToHost);
    printf("Count: %d\n", count);

    cudaFree(device_count);
    cudaFree(device_prime);
    free(prime);

    return 0;
}
@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 23, 2019

@DarioSucic 👍

Det uheldige er at CUDA sin oppstart tar ~200ms

Sant, men det blir fort mer morsomt om vi øker "n" her med et par 10-faktorer ;)

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 23, 2019

@DarioSucic

Rått parti!
1080 TI har vel 3584 CUDA cores om jeg ikke husker feil (neida, jeg har googlet) :-)

Jeg mener 5.4ms er reellt og at du kan se bort fra oppstartstid. Det gjør vi med JIT også.

@joakimwinum

This comment has been minimized.

Copy link

@joakimwinum joakimwinum commented Dec 23, 2019

AWK løsning:
awk -f door_23.awk

BEGIN {
	numberOfIterations = 98765432
	generatePrimes(digitSum(maximizeDigits(numberOfIterations)))
	for (i = 1; i <= numberOfIterations; i++) {
		number = digitSum(i)
		if (i % number == 0 && (number in primes)) {
			++count
		}
	}
}

BEGIN {
	print count
}


function digitSum(number,    numberOfDigits, i, sum)
{
	number += 0
	numberOfDigits = length(number)
	for (i = 0; i < numberOfDigits; i++) {
		sum += int(number / (10 ^ i)) % 10
	}
	return sum
}

function generatePrimes(maxInteger,    limit, i, prime)
{
	limit = int(sqrt(maxInteger))
	for (i = 2; i <= maxInteger; i++) {
		primes[i]
	}
	PROCINFO["sorted_in"] = "@ind_num_asc"
	for (prime in primes) {
		prime += 0
		if ((prime in primes) != 1) {
			continue
		}
		if (prime > limit) {
			break
		}
		for (i = prime * 2; i <= maxInteger; i += prime) {
			if (i in primes) {
				delete primes[i]
			}
		}
	}
	PROCINFO["sorted_in"] = ""
}

function maximizeDigits(number,    numberOfDigits)
{
	numberOfDigits = length(number)
	return (10 ^ (numberOfDigits + 1) - 1)
}
@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 23, 2019

@michaelo

@DarioSucic 👍

Det uheldige er at CUDA sin oppstart tar ~200ms

Sant, men det blir fort mer morsomt om vi øker "n" her med et par 10-faktorer ;)

Enig. Og denne skalerer linært. Kastet på en ekstra null til slutt og fikk 6349358 Harshadprimtall med en kjøretid på 10.305s. Enda en null og jeg passerer halvannet minutt. Estimert DarioSucic-kjøretid: 540 ms. Den raskeste boksen jeg har tilgang til har 2 x E5-2620, altså 32 kjerner. Men jeg har ikke sjanse mot 3584 CUDA-kjerner.

@DarioSucic

This comment has been minimized.

Copy link

@DarioSucic DarioSucic commented Dec 23, 2019

@DarioSucic 👍

Det uheldige er at CUDA sin oppstart tar ~200ms

Sant, men det blir fort mer morsomt om vi øker "n" her med et par 10-faktorer ;)

Hmm, det var ikke en dum idé. Ettersom vi opererer med u32 kan vi like godt sette n = 2^31 - 1 - 256.
Her trekker jeg bort 256 til så jeg slipper å tenke på overflow og blokk-størrelsesproblemer, men i praksis
er det grensen for hvor stort jeg kan beregne det uten å bytte til u64.

Uansett gir det en n som er ~21.743 ganger så stor!

Får da:
Kernel execution: 45.830 ms
Count: 6971592

edit: Dette var visst feil. Fikset litt på koden slik at den håndterer 10 siffer (grensen for u32).
Får i stedet:

Kernel execution: 103.267 ms
Count: 13650795
Total runtime: 261.788 ms

Det bekrefter i alle fall at den skalerer rimelig lineært 😄

@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 23, 2019

@DarioSucic: sånt varmer hjertet mitt 👌🏼😊

@bobcat4: nettopp. Da spiller fort 200ms forsvinnende liten rolle uansett.

@amumurst

This comment has been minimized.

Copy link

@amumurst amumurst commented Dec 23, 2019

Dagens lille scala løsning. Kjøretid ~3sek.

object Day23 extends App {
  def isPrime(i: Int): Boolean = if (i == 1 || (i != 2 && i % 2 == 0)) false else !(3.until(i, 2)).exists(i % _ == 0)
  def digitsum(i: Int, acc: Int = 0): Int = if (i == 0) acc else digitsum(i / 10, (i % 10) + acc)
  println((2 to 98765432).count(i => (i % digitsum(i) == 0) && isPrime(digitsum(i))))
}
@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 23, 2019

@DarioSucic

Hmm, det var ikke en dum idé. Ettersom vi opererer med u32 kan vi like godt sette n = 2^31 - 1 - 256.
Her trekker jeg bort 256 til så jeg slipper å tenke på overflow og blokk-størrelsesproblemer, men i praksis
er det grensen for hvor stort jeg kan beregne det uten å bytte til u64.

Uansett gir det en n som er ~21.743 ganger så stor!

Får da:

Kernel execution: 45.830 ms
Count: 6971592

Jeg får en kjøretid på 24s på denne. Og jeg får count=13648618 så vi er ikke helt enige. :-)

Loopen: for (int i = 1 ; i <= 2147483391 ; i++)

Edit: Jeg manger masse primtall i lookup-tabellen min. Duh!

Edit2: Fikset lookup-tabellen, og fikk 13650795 så vi er fremdeles ikke enige.

Edit3: Da ble vi enige ja! 😃

@DarioSucic

This comment has been minimized.

Copy link

@DarioSucic DarioSucic commented Dec 23, 2019

Da er vi nok enige! Se editen min 😉

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 23, 2019

Thumbs up til @DarioSucic som løste "ekstraoppgaven" på 46ms 261.788 ms mens jeg med min ensomme CPU-kjerne måtte bruke 24 lange sekunder.

Jeg foreslår av vi samler krefter til siste luke @terjew, og prøver å glemme denne. 🤣 Selv håper jeg siste oppgave er like lettløst som denne, for kona og jeg skal feire jul sammen med svigers og all erfaring tilsier at det ikke blir mye tid til koding i noe annet språk enn ribbe-steking. Det ble min jobb i år... 😯

Edit: Prøver å se litt mer på gårsdagens oppgave, men Sissel Kyrkjebø på NRK1 prøver etter beste evne å sette en stopper for det. 😉

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 23, 2019

@bobcat4 @michaelo
Beklager PEBCAK i sted, det var rett og slett en kopieringstabbe som gjorde at jeg målte noe helt annet enn jeg trodde jeg målte.

Relle tall er som følger:

(Kjorte enkelttrådet beregning 10 ganger på 13937 ms, snitt 1393,74 ms per kjoring)
Count: 724713
(Kjorte flertrådet beregning 10 ganger på 7464 ms, snitt 746,44 ms per kjoring)
Count: 724713

Optimaliseringene jeg trodde jeg hadde gjort viste seg bare å være tull, så da endte det til slutt med en helt rett fram løsning i C# for min del:

using Common;
using System;
using System.Linq;
using System.Threading;

namespace KnowitJulekalender23
{
    class Program
    {
        private static bool[] notPrime;

        static void Main(string[] args)
        {
            int count = 10;

            int result = 0;
            using (new CodeTimer("enkelttrådet beregning", null, count))
            {
                for (int i = 0; i < count; i++) result = SolveSinglethreaded();
            }
            Console.WriteLine($"Count: {result}");
            using (new CodeTimer("flertrådet beregning", null, count))
            {
                for (int i = 0; i < count; i++) result = SolveMultithreaded();
            }
            Console.WriteLine($"Count: {result}");
        }

        static int SolveSinglethreaded()
        {
            int largestSum = 9 + 9 + 9 + 9 + 9 + 9 + 9 + 9;
            notPrime = GeneratePrimeTruthTable(largestSum);
            int count = 0;
            for (int i = 2; i < 98765432; i++)
            {
                count += IfHarshadtPrime(i);
            }
            return count;
        }

        static int SolveMultithreaded()
        {
            int largestSum = 9 + 9 + 9 + 9 + 9 + 9 + 9 + 9;
            notPrime = GeneratePrimeTruthTable(largestSum);
            int count = Enumerable.Range(2, 98765432)
                    .AsParallel()
                    .Sum(i => IfHarshadtPrime(i));
            return count;
        }

        private static int IfHarshadtPrime(int number)
        {
            int rest = number;
            int digitsum = 0;
            while(rest > 0)
            {
                digitsum += rest % 10;
                rest /= 10;
            }
            if (number % digitsum == 0 && !notPrime[digitsum])
            {
                return 1;
            }
            return 0;
        }

        private static bool[] GeneratePrimeTruthTable(int maxPrime)
        {
            var maxSquareRoot = (int)Math.Sqrt(maxPrime);
            var eliminated = new bool[maxPrime + 1];
            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;
                        }
                    }
                }
            }
            return eliminated;
        }

    }
}
@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 23, 2019

@DarioSucic 👏
Om det var en oppgave som var laget for SIMD eller CUDA så var det vel denne :)

@bobcat4
Mistenker at grunnen til at jeg bare oppnår ca 50% forbedring ved flertrådet er aksessen av primtallstabellen, C# skjønner vel ikke at dette kan trygt gjøres uten låsing. Tipper en variant med 16 eksplisitte tråder med hver sin variant av denne vil gjøre det vesentlig bedre.

EDIT: Mistanken stemte. Ved å opprette eksplisitte tråder og la hver tråd opprette sin egen kopi av primtallstabellen kommer jeg ned i 275 ms per kjøring. Tenker jeg kaller det et greit resultat og tar kvelden :)

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 24, 2019

@terjew

@bobcat4
Mistenker at grunnen til at jeg bare oppnår ca 50% forbedring ved flertrådet er aksessen av primtallstabellen, C# skjønner vel ikke at dette kan trygt gjøres uten låsing. Tipper en variant med 16 eksplisitte tråder med hver sin variant av denne vil gjøre det vesentlig bedre.

EDIT: Mistanken stemte. Ved å opprette eksplisitte tråder og la hver tråd opprette sin egen kopi av primtallstabellen kommer jeg ned i 275 ms per kjøring. Tenker jeg kaller det et greit resultat og tar kvelden :)

Har du prøvd å bruke en read-only primtall-tabell?

Akademisk interesse. Det koster ikke stort å lage en kopi til hver tråd.

@finnmich

This comment has been minimized.

Copy link

@finnmich finnmich commented Dec 24, 2019

En liten oppdatering før jula setter inn for fullt, C# single threaded på ca. 0.5 sek.
Brukt litt inspirasjon fra andre med sieving av primtall, har også innsett at man ikke trenger å regne ut delsummen på nytt hvis tallet ikke er delelig på 10, da kan man bare legge til 1 på forrige delsum.

public class Luke23
{
    private static bool[] primeTable;
    public static long Run()
    {
        primeTable = GeneratePrimeTruthTable(81);

        var count = 0;
        var currentSum = 0;
        for (int i = 1; i < 98765432; i++)
        {
            if (i % 10 != 0)
            {
                currentSum++;
                    
            }
            else
            {
                var n = i;
                currentSum = 0;
                while (n > 0)
                {
                    currentSum += n % 10;
                    n /= 10;
                }
            }

            if (i % currentSum == 0 && !primeTable[currentSum])
            {
                count++;
            }

        }
        return count;
    }

    private static bool[] GeneratePrimeTruthTable(int maxPrime)
    {
        var maxSquareRoot = (int)Math.Sqrt(maxPrime);
        var eliminated = new bool[maxPrime + 1];
        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;
                    }
                }
            }
        }
        return eliminated;
    }
}
@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 24, 2019

@finnmich den optimaliseringen satt jeg og funderte på i bilen fra Oslo til Telemark nå, artig at noen andre hadde implementert den :)

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 26, 2019

Rakk ikke utforske denne helt på lille julaften, så jeg tok den opp igjen nå for å ta en kikk.

Har tenkt å implementere samme optimalisering som @finnmich, men før det ville jeg teste å tune antall tråder, og ser at ved å gå fra 16 til 12 tråder oppnår jeg ganske mye bedre resultat; 230ms. Det var jo en ganske billig optimalisering, og gir mening siden jeg har en CPU med 6 kjerner, hver av dem med Hyperthreading. Det er ingen venting på I/O eller liknende i denne oppgaven, så hver kjerne bør kunne være busy 100% gjennom hele kjøringen, samtidig er det en del aritmetikk der hyperthreadingen sannsynligvis gir avkastning.

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 26, 2019

Ved å i tillegg ta inn samme optimalisering som @finnmich brukte på sin løsning havner jeg på:

(Kjorte flertrådet beregning (eksplisitt tråding) 100 ganger på 7265 ms, snitt 72,65 ms per kjoring)
Count: 724713

Ikke helt CUDA-territorie, men tross alt ikke SÅ langt unna heller...

Edit: Og én enkeltkjøring målt fra shellet havner på ca 140 ms, så litt raskere enn CUDA der faktisk :)
Edit2: Ved å flippe sannhetstabellen fra sievingen endte jeg med en 10% forbedring igjen

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 26, 2019

Prøvde meg på en SIMD-variant vha .Net core Numerics også, men uten en modulo-operator er jeg litt usikker på hvordan jeg skal få det til å spinne særlig effektivt.

Prøvde meg på å implementere modulo selv som A % B = A - (B * (A / B)), men da blir det ganske mange operasjoner for hvert tall, så det ender opp ca 4 ganger så tregt som en naiv variant:

                    Vector<int> vPowers = new Vector<int>(new[] { 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000 }, 0);
                    Vector<int> v10 = new Vector<int>(new[] { 10, 10, 10, 10, 10, 10, 10, 10 }, 0);
                    Vector<int> v1 = new Vector<int>(new[] { 1, 1, 1, 1, 1, 1, 1, 1}, 0);
                    Vector<int> number = new Vector<int>(i);
                    var divided = Vector.Divide(number, vPowers);

                    //calculate modulo 10:
                    var d10 = Vector.Divide(divided, v10);
                    var m10 = Vector.Multiply(d10, v10);
                    var res = Vector.Subtract(divided, m10);
                    //sum of digits:
                    digitSum = Vector.Dot(res, v1);

Noen som har et forslag til hvordan man kan gjøre dette mer effektivt uten modulo?

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 26, 2019

@DarioSucic @michaelo @bobcat4
Prøvde meg på det store datasettet (2^31 - 1 - 256) og fikk følgende:

(Kjorte flertrådet beregning (eksplisitt tråding) 10 ganger på 14379 ms, snitt 1437,96 ms per kjoring)
Count: 13650795

Altså ca 15x kjøretiden til CUDA der også.

Her er (foreløpig?) siste versjon:

        public static int SolveMultithreadedExplicit(int max)
        {
            int numthreads = 16;
            int partitionsize = (max - 2) / numthreads;
            int start = 2;
            Thread[] threads = new Thread[numthreads];
            int[] sums = new int[numthreads];
            for (int threadNo = 0; threadNo < numthreads; threadNo++)
            {
                int threadBegin = start;
                int threadEnd = Math.Min(start + partitionsize, max);
                int currentThread = threadNo;
                var thread = new Thread(new ThreadStart(() =>
                {
                    sums[currentThread] = SolveRangeOptimized(threadBegin, threadEnd);
                }));
                thread.Start();
                threads[threadNo] = thread;
                start += partitionsize;
            }

            int sum = 0;
            for (int threadNo = 0; threadNo < numthreads; threadNo++)
            {
                threads[threadNo].Join();
                sum += sums[threadNo];
            }
            return sum;
        }

        private static int SolveRangeOptimized(int start, int end)
        {
            var largestPossibleSum = (int)(Math.Round(Math.Log10(end) + 1) * 9);
            bool[] prime = GeneratePrimeTruthTable(largestPossibleSum, false);
            int count = 0;
            int digitsum = 0;
            int rest = start;
            while (rest > 0)
            {
                digitsum += rest % 10;
                rest /= 10;
            }
            if (prime[digitsum] && start % digitsum == 0)
            {
                count += 1;
            }

            for (int i = start + 1; i < end; i++)
            {
                if (i % 10 != 0) digitsum++;
                else
                {
                    digitsum = 0;
                    rest = i;
                    while (rest > 0)
                    {
                        digitsum += rest % 10;
                        rest /= 10;
                    }
                }
                if (prime[digitsum] && i % digitsum == 0)
                {
                    count += 1;
                }
            }
            return count;
        }

        
        private static bool[] GeneratePrimeTruthTable(int maxPrime, bool inverted = true)
        {
            var maxSquareRoot = (int)Math.Sqrt(maxPrime);
            var eliminated = new bool[maxPrime + 1];
            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)
            {
                return eliminated;
            }
            else
            {
                return eliminated.Select(b => !b).ToArray();
            }
        }
@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 26, 2019

inb4 @bobcat4 tar C#-koden min og kompilerer den som C++ 🤣

Her er C++-varianten av samme kode.

Lite datasett: 53ms (snitt over 100 kjøringer)
Stort datasett: 1171ms (snitt over 10 kjøringer)

Da er vi vel på ca 10-11 ganger kjøretiden til CUDA.

#include <iostream>
#include <thread>
#include <algorithm>

bool * GeneratePrimeTruthTable(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;
}

int SolveRangeOptimized(int start, int end)
{
  auto largestPossibleSum = (int)(round(log10(end) + 1) * 9);
  bool* prime = GeneratePrimeTruthTable(largestPossibleSum, false);
  int count = 0;
  int digitsum = 0;
  int rest = start;
  while (rest > 0)
  {
    digitsum += rest % 10;
    rest /= 10;
  }
  if (prime[digitsum] && start % digitsum == 0)
  {
    count += 1;
  }

  for (int i = start + 1; i < end; i++)
  {
    if (i % 10 != 0) digitsum++;
    else
    {
      digitsum = 0;
      rest = i;
      while (rest > 0)
      {
        digitsum += rest % 10;
        rest /= 10;
      }
    }
    if (prime[digitsum] && i % digitsum == 0)
    {
      count += 1;
    }
  }
  delete[] prime;
  return count;
}

int SolveMultithreadedExplicit(int max)
{
  int numthreads = 12;
  int partitionsize = (max - 2) / numthreads;
  int start = 2;
  std::thread * threads = new std::thread[numthreads];
  int * sums = new int[numthreads];

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

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

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

int main()
{
  int count = 10;
  int result = 0;
  for (int i = 0; i < count; i++)
    //result = SolveRangeOptimized(2, 98765432);
    //result = SolveMultithreadedExplicit(98765432);
    result = SolveMultithreadedExplicit(INT_MAX - 256);
  std::cout << "result is : " << result << "\n";
}
@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 26, 2019

@bobcat4 @michaelo

Lukter "false sharing" av denne MT-en. Du bør kunne komme under 100 ms med 16 tråder.
Bør helt klart kunne komme lavere enn jeg er pr nå, men ikke så lavt som du postulerer. i7 7700HQ'en min har 4 cores / 8 threads. Har enda ikke fått pinpointet hvorfor den faktisk skreller av ~30 ms når jeg går fra 8 til 16 tråder heller.

@bobcat4 Kan du forklare hva du mener at indikerer "false sharing" i det caset?

Med en CPU med 4 "ekte" kjerner pluss hyperthreading vil jeg gjette at en best-case forbedring ved MT i de fleste tilfeller ligger på rundt 5X.
4X skal være ganske bankers, pluss at man gjerne får litt ekstra fra HT men neppe så mye som 8X. Intel har vel tidligere hintet om "opp til 30%" forbedring med HT kontra uten, i praksis mener jeg det pleier å ligge en del lavere.

@michaelo Jeg ville gjettet at det ideelle antall tråder ligger et sted mellom 4 og 8 siden det ikke er noen venting på IO. Hos meg så jeg best ytelse på 12 tråder (6 kjerner pluss HT), men ikke stor forskjell fra 9 til 12. Om du jevnt over får bedre ytelse på 16 tråder enn 8 synes jeg det er noe pussig på gang...

@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 26, 2019

@terjew

@michaelo Jeg ville gjettet at det ideelle antall tråder ligger et sted mellom 4 og 8 siden det ikke er noen venting på IO. Hos meg så jeg best ytelse på 12 tråder (6 kjerner pluss HT), men ikke stor forskjell fra 9 til 12. Om du jevnt over får bedre ytelse på 16 tråder enn 8 synes jeg det er noe pussig på gang...

Jeg er helt enig med deg. Dette var resultat av fritthåndskjøringer i et hverdagsbelastet system - på ingen måter akademisk tilnærmet. Mente å få tatt mer kontrollerte målinger, men falt bakpå der blant gamle og nye luker og jul. Men nå er det nye muligheter for forsøk!

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 27, 2019

@terjew

@bobcat4 Kan du forklare hva du mener at indikerer "false sharing" i det caset?

Det var vel denne:

threads[i].prime_map = NUM_IS_PRIME_MAP;

Kan vi stole på at compileren er bombesikker på at dette er en const ?

I utgangspunktet skal ikke en identisk kopi av prime_map (som aldri endres) føre til "dirty cache lines", men ettersom ytelsen ikke samsvarte med antall tråder... Men jeg har tatt feil før. :-)

For øvrig: Jeg har oppgradert Qemu til 4.2.0 og også oppgradert "host"-ens kernel så alle gamle målingene mine gjelder ikke dessverre.

@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 27, 2019

@terjew:

Da er vi vel på ca 10-11 ganger kjøretiden til CUDA.

Pent 👍

@bobcat4

threads[i].prime_map = NUM_IS_PRIME_MAP;

Jeg testet en del med varianter av pointer-indirection og const-varianter for å se hva den hang seg opp i, og testet også med å sette denne i fil-scope uten at det gjorde noe nevneverdig utslag i denne oppgaven. Likeledes med lokal full-kopi forøvrig.

Oppdaterte målinger på min gamle løsning, bare for å bekrefte at det var muffens sist. Nå flater det ut ved 8 tråder som (egentlig) forventet:

tråder: tid
1: 1.042s
2:   498ms ± 7.9ms
4:   298ms ± 3.4ms
8:   224ms ± 5.5ms
12:  222ms ± 1.7ms
16:  224ms ± 4.5ms
@DarioSucic

This comment has been minimized.

Copy link

@DarioSucic DarioSucic commented Dec 27, 2019

Økte ytelsen på min maskin med 20% ved å bytte til unsigned int.

int SolveRangeOptimized(int start, int end)
{
  ...
  int count = 0;
  int digitsum = 0;
  int rest = start;
  ...
}

til

int SolveRangeOptimized(unsigned int start, unsigned int end)
{
  ...
  unsigned int count = 0;
  unsigned int digitsum = 0;
  unsigned int rest = start;
  ...
}

Hvis vi ser på den genererte maskinkoden ser vi fort hvorfor (Clang 9.0).
Gitt funksjonene:

int i32mod(int a) {
    return a % 10;
}

unsigned int u32mod(unsigned int a) {
    return a % 10;
}

Får vi følgende maskinkode:

i32mod(int):                             # @i32mod(int)
        movsxd  rax, edi
        imul    rcx, rax, 1717986919
        mov     rdx, rcx
        shr     rdx, 63
        sar     rcx, 34
        add     ecx, edx
        add     ecx, ecx
        lea     ecx, [rcx + 4*rcx]
        sub     eax, ecx
        ret

u32mod(unsigned int):                             # @u32mod(unsigned int)
        mov     eax, edi
        mov     ecx, edi
        mov     edx, 3435973837
        imul    rdx, rcx
        shr     rdx, 35
        add     edx, edx
        lea     ecx, [rdx + 4*rdx]
        sub     eax, ecx
        ret
@DarioSucic

This comment has been minimized.

Copy link

@DarioSucic DarioSucic commented Dec 27, 2019

Lekte litt med en metode for å generere siffer-summene uten å bruke noen delingsinstruksjoner, og endte opp med det her.
Kjøretiden på min laptop ligger på rundt 8s for 10 kjøringer, men varierer noe veldig. (lavest observert ~7.6s)

Kompilert med clang main.cpp -Ofast -march=skylake og kjørt på en 7700HQ.

#include <iostream>
#include <thread>
#include <algorithm>
#include <math.h>

#include <stdint.h>

#define u8 uint8_t
#define u32 uint32_t
#define i32 int32_t

bool *GeneratePrimeTruthTable(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;
}

u32 SolveRangeOptimized(u32 start, u32 end)
{

    u32 nDigits = (u32)(ceil(log10(end)));
    u32 largestPossibleSum = 9 * nDigits;
    bool *prime = GeneratePrimeTruthTable(largestPossibleSum, false);

    u32 count = 0;
    u32 *digits = new u32[nDigits];
    u32 s = 0;
    u32 _start = start;

    for(u32 d = 0; d < nDigits; ++d)
    {
        digits[nDigits - 1 - d] = _start % 10;
        s += digits[nDigits - 1 - d];
        _start /= 10;
    }

    u32 lastDigit = nDigits - 1;
    u32 digitIdx = lastDigit;
    u32 idx = start;

    s -= digits[lastDigit];

    while(true)
    {
        if(digitIdx == lastDigit)
        {
            for(; digits[lastDigit] < 10; digits[lastDigit]++)
            {
                u32 digitSum = s + digits[lastDigit];

                if(prime[digitSum] && idx % digitSum == 0)
                    ++count;

                if(idx == end)
                    goto END;
                
                ++idx;
            }
            digitIdx -= 1;
            digits[lastDigit] = 0;
        }
        else
        {

            if(digits[digitIdx] == 9)
            {
                digits[digitIdx] = -1;
                digitIdx -= 1;
                s -= 10;
            }
            else
            {
                digits[digitIdx] += 1;
                digitIdx += 1;
                s += 1;
            }
        }
    }

    END:

    delete[] prime;
    delete[] digits;
    return count;
}

int SolveMultithreadedExplicit(int max)
{
    int numthreads = 16;
    int partitionsize = (max - 2) / numthreads;
    int start = 2;
    std::thread *threads = new std::thread[numthreads];
    int *sums = new int[numthreads];

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

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

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

int main()
{
    int count = 10;
    int result = 0;
    for (int i = 0; i < count; i++)
        // result = SolveRangeOptimized(2, 98765432);
        //result = SolveMultithreadedExplicit(98765432);
        result = SolveMultithreadedExplicit(2147483647 - 256);
    std::cout << "result is : " << result << "\n";
}
@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 27, 2019

@DarioSucic

Hvis vi ser på den genererte maskinkoden ser vi fort hvorfor (Clang 9.0).
Gitt funksjonene:

int i32mod(int a) {
    return a % 10;
}

unsigned int u32mod(unsigned int a) {
    return a % 10;
}

Får vi følgende maskinkode:

i32mod(int):                             # @i32mod(int)
        movsxd  rax, edi
        imul    rcx, rax, 1717986919
        mov     rdx, rcx
        shr     rdx, 63
        sar     rcx, 34
        add     ecx, edx
        add     ecx, ecx
        lea     ecx, [rcx + 4*rcx]
        sub     eax, ecx
        ret

u32mod(unsigned int):                             # @u32mod(unsigned int)
        mov     eax, edi
        mov     ecx, edi
        mov     edx, 3435973837
        imul    rdx, rcx
        shr     rdx, 35
        add     edx, edx
        lea     ecx, [rdx + 4*rdx]
        sub     eax, ecx
        ret

GCC:

#int32_t i32mod(int32_t a) 
   0:   48 63 c7                movslq %edi,%rax
   3:   89 fa                   mov    %edi,%edx
   5:   48 69 c0 67 66 66 66    imul   $0x66666667,%rax,%rax
   c:   c1 fa 1f                sar    $0x1f,%edx
   f:   48 c1 f8 22             sar    $0x22,%rax
  13:   29 d0                   sub    %edx,%eax
  15:   8d 04 80                lea    (%rax,%rax,4),%eax
  18:   01 c0                   add    %eax,%eax
  1a:   29 c7                   sub    %eax,%edi
  1c:   89 f8                   mov    %edi,%eax
  1e:   c3                      retq
  1f:   90                      nop

0000000000000020 <u32mod>:

# uint32_t u32mod(uint32_t a)
  20:   89 f8                   mov    %edi,%eax
  22:   ba cd cc cc cc          mov    $0xcccccccd,%edx
  27:   48 0f af c2             imul   %rdx,%rax
  2b:   48 c1 e8 23             shr    $0x23,%rax
  2f:   8d 04 80                lea    (%rax,%rax,4),%eax
  32:   01 c0                   add    %eax,%eax
  34:   29 c7                   sub    %eax,%edi
  36:   89 f8                   mov    %edi,%eax
  38:   c3                      retq

Ikke så rart at unsigned er raskere enn signed (selv uten å se assembly-en). Det er ikke gratis å holde styr på fortegn (det vet alle som har lært algebra) ;-)

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 27, 2019

Hmm... Jeg testet uint32_t vs int og int er faktisk raskere på min CPU:

int inline is_Harshadprime (int num)
{
        int n = num;
        int sum = 0;
        while (n != 0)
        {
                int rem = n % 10;
                sum += rem;
                n /= 10;
        }
//      return (num%sum == 0) && primes[sum];
        return  primes[sum] && (num%sum == 0);
}

Jeg fant også ut at return primes[sum] && (num%sum == 0); er raskere enn return (num%sum == 0) && primes[sum];, men det er vel ganske logisk når jeg tenker litt mer på det. Dette er single threaded. Kom ned i 750ms med int og sjekking av primtall før mod.

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 27, 2019

@bobcat4

Jeg fant også ut at return primes[sum] && (num%sum == 0); er raskere enn return (num%sum == 0) && primes[sum];, men det er vel ganske logisk når jeg tenker litt mer på det. Dette er single threaded. Kom ned i 750ms med int og sjekking av primtall før mod.

Jeg fant samme optimalisering i går, du ser det i koden jeg postet i natt. Ganske dramatisk forskjell faktisk, modulo-operatoren er vesentlig tregere enn jeg først hadde antatt.

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 27, 2019

@DarioSucic

Lekte litt med en metode for å generere siffer-summene uten å bruke noen delingsinstruksjoner, og endte opp med det her.
Kjøretiden på min laptop ligger på rundt 8s for 10 kjøringer, men varierer noe veldig. (lavest observert ~7.6s)

👏 👏 👏

Fantastisk! Jeg visste det var mer å hente i den koden jeg postet, herlig å se at andre klarer å finne veien videre når man er litt stuck.

Overrasket over at det var så stor forskjell på uint og int, men det gir jo mening med assemblyen du har postet. Dette må jeg teste på min maskin hhv med C# og C++ (med cl.exe) senere.

Jeg var også inne på tanken å holde på arrayet over sifre og vedlikeholde disse på lik linje som det du gjør, men tenkte at det kan umulig gi noen gevinst når vi egentlig bare trenger å beregne tverrsummen for hver 10. iterasjon... Hva slags ytelse fikk du på det store "datasettet" før den optimaliseringen?

Må gjøre litt betalt jobbing noen timer, men kommer tilbake senere og ser hvor langt dere har kommet! 😄

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 27, 2019

Jeg fant samme optimalisering i går, du ser det i koden jeg postet i natt. Ganske dramatisk forskjell faktisk, modulo-operatoren er vesentlig tregere enn jeg først hadde antatt.

Jeg trodde det ikke før jeg så det. 😄 Magefølelsen sier at det er feil, men det er jo ikke det. Det er faktisk ganske logisk.

modulo før primes-lookup:

            963.95 msec task-clock                #    1.000 CPUs utilized
                 2      context-switches          #    0.002 K/sec
                 0      cpu-migrations            #    0.000 K/sec
                56      page-faults               #    0.058 K/sec
     3,750,971,822      cycles                    #    3.891 GHz
     8,900,744,843      instructions              #    2.37  insn per cycle
       977,230,258      branches                  # 1013.778 M/sec
         6,946,132      branch-misses             #    0.71% of all branches

primes-lookup før modulo:

            760.08 msec task-clock                #    1.000 CPUs utilized
                 6      context-switches          #    0.008 K/sec
                 0      cpu-migrations            #    0.000 K/sec
                57      page-faults               #    0.075 K/sec
     2,959,304,796      cycles                    #    3.893 GHz
     8,903,189,912      instructions              #    3.01  insn per cycle
       977,120,518      branches                  # 1285.556 M/sec
         6,814,130      branch-misses             #    0.70% of all branches

Her blir det rett og slett færre instruksjoner av å sjekke primtall først.
Det er altså større sannsynlighet for at primes[sum] == 0 enn at (num%sum == 0)

Men... Med -fprofile-generate og -fprofile-use så skjønner compiler-en hvor landet ligger og bytter om på sjekkene. Imponerende!

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 27, 2019

@michaelo

threads[i].prime_map = NUM_IS_PRIME_MAP;

Jeg testet en del med varianter av pointer-indirection og const-varianter for å se hva den hang seg opp i, og testet også med å sette denne i fil-scope uten at det gjorde noe nevneverdig utslag i denne oppgaven. Likeledes med lokal full-kopi forøvrig.

Jeg stusset også over at denne skulle gi noe false sharing da den ikke skrives til.

Jeg var et øyeblikk inne på tanken at variabelen du akkumulerer til kanskje lå på samme cachelinje for de ulike trådenes state. Dersom du hadde oppdatert en teller direkte på state-structen kunne man jo se for seg det, men siden du bruker en lokal variabel i trådenes mt_worker-funksjon vil jeg jo anta at denne havner i et register.

Oppdaterte målinger på min gamle løsning, bare for å bekrefte at det var muffens sist. Nå flater det ut ved 8 tråder som (egentlig) forventet:

Da høres det ut som ting er tilbake til normalen ja 😌

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 27, 2019

@bobcat4

Jeg trodde det ikke før jeg så det. 😄 Magefølelsen sier at det er feil, men det er jo ikke det. Det er faktisk ganske logisk.
Her blir det rett og slett færre instruksjoner av å sjekke primtall først.

Marginalt færre, men ja. Men det blir også (marginalt) mindre branching og ~2% færre branch misses, noe som kanskje kan være vel så viktig?

Det er altså større sannsynlighet for at primes[sum] == 0 enn at (num%sum == 0)

Det er mulig at det er riktig, men jeg kan ikke på stående fot si at det kan utledes fra setningen over... Det kan vel også forklares av at primes[sum] == 0 krever vesentlig færre instruksjoner å sjekke enn (num%sum == 0), selv om primes[sum] ofte skulle returnere 1/true?

Men... Med -fprofile-generate og -fprofile-use så skjønner compiler-en hvor landet ligger og bytter om på sjekkene. Imponerende!

Det var faktisk litt gøy å høre :) Spent på om MSVC catcher den med "whole program optimization" (/GL), skal teste senere.

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 27, 2019

@terjew

Det er altså større sannsynlighet for at primes[sum] == 0 enn at (num%sum == 0)

Det er mulig at det er riktig, men jeg kan ikke på stående fot si at det kan utledes fra setningen over...

En matematiker kan nok skrive side opp og side ned om dette, men jeg utledet det fra antall "cycles". 😉 Rett og slett: Mindre jobb å sjekke primes[sum] først.

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 27, 2019

@bobcat4

Det er altså større sannsynlighet for at primes[sum] == 0 enn at (num%sum == 0)

Det er mulig at det er riktig, men jeg kan ikke på stående fot si at det kan utledes fra setningen over...

En matematiker kan nok skrive side opp og side ned om dette, men jeg utledet det fra antall "cycles". 😉 Rett og slett: Mindre jobb å sjekke primes[sum] først.

Men er ikke "cycles" bare tid * klokkefrekvens? Det at de to variantene har ca like mange instruksjoner på ganske ulikt antall klokkesykler ville jeg tro kan skyldes at:

  1. den ene varianten i snitt utfører dyrere instruksjoner. Noen matte-instruksjoner tar jo f.eks flere klokkesykler enn andre, og ved riktig unrolling kan noen ganger prosessoren utføre flere "vanlige" (ikke SIMD) matteinstruksjoner på samme klokkesyklus.

eller

  1. den ene varianten får flere pipeline flusher, noe som vel typisk skjer ved feil branch prediction.

Det er ikke meningen å drive med flisespikking her, men jeg er langt fra noen ekspert på dette temaet og er interessert i å lære om min mentale modell av dette stemmer med virkeligheten.

Konklusjonen er jeg uansett enig i, mindre jobb å sjekke primes først :)

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 27, 2019

En matematiker kan nok skrive side opp og side ned om dette, men jeg utledet det fra antall "cycles". wink Rett og slett: Mindre jobb å sjekke primes[sum] først.

Men er ikke "cycles" bare tid * klokkefrekvens? Det at de to variantene har ca like mange instruksjoner på ganske ulikt antall klokkesykler ville jeg tro kan skyldes at:

  1. den ene varianten i snitt utfører dyrere instruksjoner. Noen matte-instruksjoner tar jo f.eks flere klokkesykler enn andre, og ved riktig unrolling kan noen ganger prosessoren utføre flere "vanlige" (ikke SIMD) matteinstruksjoner på samme klokkesyklus.

eller

  1. den ene varianten får flere pipeline flusher, noe som vel typisk skjer ved feil branch prediction.

Ja, du har et poeng her, @terjew. Det var jeg som tulla med "cycles" som er tid*klokkefrekvens (single threaded). Det VM i Moskva, og hjernen er i sjakk-modus.

Jeg vil tro at alternativ 2 gjelder her. Pipelining (og her snakker vi om en brukbar bunke rør) er en vesentlig faktor, og pipeline flush koster veldig mye. Google min i5-4690K som en CPU med "14- to 19-stage instruction pipeline". Går man tom (på grunn av flush) vil det koste mye.

Det er ikke meningen å drive med flisespikking her, men jeg er langt fra noen ekspert på dette temaet og er interessert i å lære om min mentale modell av dette stemmer med virkeligheten.

Flisespikking er morsomt, og absolutt min gren! 😉

Konklusjonen er jeg uansett enig i, mindre jobb å sjekke primes først :)

Dette er jo engentlig fryktelig interessant! Jeg måtte sjekke naturligvis.

(num%sum == 0) er false 91876607 ganger,
primes[sum] er false 74273000 ganger.

Ja, da tok jeg feil (igjen). Så da er vel mod så kostbart at det lønner seg å sjekke primtall 17 millioner flere ganger... :-)

Edit: Ser ut som om GCC bruker IDIV for mod10 mod sum og det er vel ikke særlig smart... Latency 22-29 på en Haswell-CPU.
Edit2: Ettersom sum (tverrsummen) er en variabel er IDIV beste alternativ.

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 27, 2019

@DarioSucic @bobcat4

Økte ytelsen på min maskin med 20% ved å bytte til unsigned int.

På min maskin med VS2019 (cl.exe 19.24.28314) gir også unsigned int en besparelse, men her er det kun ~10%.

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 27, 2019

Dette er sikker ikke noe nytt...

Tverrsummen ser ut til å være det som koster mest. Så finnes det noen god metode å finne tverrsummen uten å gå igjennom alle tall?

Altså, uten å bruke

    while (n != 0) 
    { 
     sum = sum + n % 10; 
     n = n/10; 
    } 

Hvis vi teller fra 1 til 1000 og skriver ut tverrsummen i rader på 9 tall dukker det opp et mønster:

 1, 2, 3, 4, 5, 6, 7, 8, 9,
 1, 2, 3, 4, 5, 6, 7, 8, 9,
10, 2, 3, 4, 5, 6, 7, 8, 9,
10,11, 3, 4, 5, 6, 7, 8, 9,
10,11,12, 4, 5, 6, 7, 8, 9,
10,11,12,13, 5, 6, 7, 8, 9,
10,11,12,13,14, 6, 7, 8, 9,
10,11,12,13,14,15, 7, 8, 9,
10,11,12,13,14,15,16, 8, 9,
10,11,12,13,14,15,16,17, 9,
10,11,12,13,14,15,16,17,18,
 1, 2, 3, 4, 5, 6, 7, 8, 9,
10, 2, 3, 4, 5, 6, 7, 8, 9,
10,11, 3, 4, 5, 6, 7, 8, 9,
10,11,12, 4, 5, 6, 7, 8, 9,
10,11,12,13, 5, 6, 7, 8, 9,
10,11,12,13,14, 6, 7, 8, 9,
10,11,12,13,14,15, 7, 8, 9,
10,11,12,13,14,15,16, 8, 9,
10,11,12,13,14,15,16,17, 9,
10,11,12,13,14,15,16,17,18,
10,11,12,13,14,15,16,17,18,
19, 2, 3, 4, 5, 6, 7, 8, 9,
10,11, 3, 4, 5, 6, 7, 8, 9,
10,11,12, 4, 5, 6, 7, 8, 9,
10,11,12,13, 5, 6, 7, 8, 9,
10,11,12,13,14, 6, 7, 8, 9,
10,11,12,13,14,15, 7, 8, 9,
10,11,12,13,14,15,16, 8, 9,
10,11,12,13,14,15,16,17, 9,
10,11,12,13,14,15,16,17,18,
10,11,12,13,14,15,16,17,18,
19,11,12,13,14,15,16,17,18,
19,20, 3, 4, 5, 6, 7, 8, 9,
10,11,12, 4, 5, 6, 7, 8, 9,
10,11,12,13, 5, 6, 7, 8, 9,
10,11,12,13,14, 6, 7, 8, 9,
10,11,12,13,14,15, 7, 8, 9,
10,11,12,13,14,15,16, 8, 9,
10,11,12,13,14,15,16,17, 9,
10,11,12,13,14,15,16,17,18,
10,11,12,13,14,15,16,17,18,
19,11,12,13,14,15,16,17,18,
19,20,12,13,14,15,16,17,18,
19,20,21, 4, 5, 6, 7, 8, 9,
10,11,12,13, 5, 6, 7, 8, 9,
10,11,12,13,14, 6, 7, 8, 9,
10,11,12,13,14,15, 7, 8, 9,
10,11,12,13,14,15,16, 8, 9,
10,11,12,13,14,15,16,17, 9,
10,11,12,13,14,15,16,17,18,
10,11,12,13,14,15,16,17,18,
19,11,12,13,14,15,16,17,18,
19,20,12,13,14,15,16,17,18,
19,20,21,13,14,15,16,17,18,
19,20,21,22, 5, 6, 7, 8, 9,
10,11,12,13,14, 6, 7, 8, 9,
10,11,12,13,14,15, 7, 8, 9,
10,11,12,13,14,15,16, 8, 9,
10,11,12,13,14,15,16,17, 9,
10,11,12,13,14,15,16,17,18,
10,11,12,13,14,15,16,17,18,

"Det holder det, sa Gundersen, det klarer seg med det" (Arne Svendsen) 👎

Det ser i hvert fall ut som et mønster, men det er neppe rett matrise. Kan man finne en "billig" formel som gir samme output?

Slik ser det ut i (ref https://oeis.org/A007953
graph

Ligner på en rekke parallellogram. Burde være grei skuring..

Formelen fra Wikipedia virker ikke særlig smart (vi vil ikke ha divisjon)

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 27, 2019

@bobcat4

Tverrsummen ser ut til å være det som koster mest. Så finnes det noen god metode å finne tverrsummen uten å gå igjennom alle tall?

Du har helt rett, det er kostbart å beregne den, så vi bør unngå det så langt det er mulig.

Jeg har ikke kommet opp med noen algoritme eller mønster som kan beregne alle tverrsummer, men med @finnmich sin optimalisering kan vi kutte 90% av beregningene da vi vet at så lenge (i % 10 != 0) så vil alltid tverrsum være lik forrige tverrsum +1. Dermed trenger vi bare beregne tverrsum på nytt for hver 10. iterasjon, noe som ga en stor besparelse (totalt rundt 50%) fra opprinnelig algoritme.

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 27, 2019

@terjew

Jeg har ikke kommet opp med noen algoritme eller mønster som kan beregne alle tverrsummer, men med @finnmich sin optimalisering kan vi kutte 90% av beregningene da vi vet at så lenge (i % 10 != 0) så vil alltid tverrsum være lik forrige tverrsum +1. Dermed trenger vi bare beregne tverrsum på nytt for hver 10. iterasjon, noe som ga en stor besparelse (totalt rundt 50%) fra opprinnelig algoritme.

Ja, jeg så den, men jeg trodde det var mulig å gjøre det bedre. Etter utallige Google-søk har jeg endret mening. Enten er ikke jeg flink nok med Google eller så finnes ikke en sånn løsning. 😃

Edit: med @finnmich sin optimalisering er jeg nede i 224 ms i én tråd.

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 27, 2019

@DarioSucic @bobcat4

Varianten som beregner sifferene uten divisjon kjører hos meg på 10,2 sekunder for 10 kjøringer, altså 1,02 sekunder per kjøring i snitt.

Jeg valgte en annen tilnærming og gikk i stedet løs på loopen og begynte å unrolle.

Først unrollet jeg while-løkken som beregner tverrsum, og har foreløpig valgt å alltid kjøre 10 iterasjoner på denne. Det vil nok være lønnsomt å holde styr på antall siffer i segmentet man jobber med og kjøre ulike varianter avhengig av det, men med 10 iterasjoner funker det i hvert fall for alle tall som passer i en int32. For tall med få siffer er det selvsagt massiv overhead, men for tall med 10 siffer sparer det en del branching. Om jeg gidder skal jeg lage egen variant for 7,8,9 siffer også. Alle tall under 7 siffer kan egentilg håndteres på "gamlemåten", de tallene utgjør så liten prosent av totalen uansett.

I tillegg unrollet jeg selve telle-loopen slik at den kjører 10 og 10 iterasjoner. Da slipper vi teste i % 10 for hver gjennomkjøring, og kan i tillegg bruke informasjon om hvorvidt tverrsummen er partall til å skippe halvparten av testene for primtall.

Med disse endringene kjører programmet nå 10 iterasjoner av "stort datasett" på ~5,9 sekunder på min maskin, altså 590 ms per iterasjon.

EDIT: Tweaket litt på unrollingen og oppdaterte tallene.

inline void step(uint32_t& digitsum, uint32_t& count, bool* prime, int i)
{
  if (prime[digitsum] && i % digitsum == 0)
  {
    count++;
  }
}

inline int sumDigits(uint32_t i)
{
  uint32_t digitsum = 0;
  uint32_t rest = i;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  return digitsum;
}

int SolveRangeOptimized(int start, int end)
{
  auto largestPossibleSum = (int)(round(log10(end) + 1) * 9);
  bool* prime = GeneratePrimeTruthTable(largestPossibleSum, false);
  uint32_t count = 0;
  uint32_t digitsum = 0;
  uint32_t rest = start;

  int blockend = (end / 10) * 10;

  int i = start;
  //head:
  do {
    digitsum = sumDigits(i);
    step(digitsum, count, prime, i++);
  } while (i % 10 != 0);

  for (; i < blockend; i += 10)
  {
    digitsum = sumDigits(i);
    step(digitsum, count, prime, i);
    if (digitsum % 2 != 0) {
      digitsum += 2;
      step(digitsum, count, prime, i + 2);
      digitsum += 2;
      step(digitsum, count, prime, i + 4);
      digitsum += 2;
      step(digitsum, count, prime, i + 6);
      digitsum += 2;
      step(digitsum, count, prime, i + 8);
    }
    else {
      digitsum += 1;
      step(digitsum, count, prime, i + 1);
      digitsum += 2;
      step(digitsum, count, prime, i + 3);
      digitsum += 2;
      step(digitsum, count, prime, i + 5);
      digitsum += 2;
      step(digitsum, count, prime, i + 7);
      digitsum += 2;
      step(digitsum, count, prime, i + 9);
    }
  }

  //tail:
  while (i < end) {
    digitsum = sumDigits(i);
    step(digitsum, count, prime, i++);
  }

  delete[] prime;
  return count;
}

@DarioSucic

This comment has been minimized.

Copy link

@DarioSucic DarioSucic commented Dec 27, 2019

@bobcat4

Tverrsummen ser ut til å være det som koster mest. Så finnes det noen god metode å finne tverrsummen uten å gå igjennom alle tall?

Mønsteret er fler-periodisk, og er i grunnen ganske enkelt å beregne med nøstede for-løkker for et fiksert antall siffer, eller rekursivt for et arbitrært antall siffer.

Her er et eksempel på det sistnevnte:

def digit_sums(n_digits):
    def rec(prev, digit_index):
        if digit_index == n_digits - 1:
            for c in range(10):
                print(prev + c)
        else:
            for c in range(10):
                rec(prev + c, digit_index+1)

    return rec(0, 0)

Og en versjon med for-løkker kunne sett slik ut for 3 siffer:

for a in range(10):
    aa = a * 10**2
    for b in range(10):
        bb = aa + b * 10**1
        for c in range(10):
            cc = bb + c * 10**0
            print(cc)

Det som er kult er at denne også tar hensyn til @finnmich sin optimalisering (synlig i den innerste løkken).
Problemet er at det blir vanskelig å replikere strukturen ovenfor for et arbitrært antall siffer.

Jeg løste det ved å lage en tilstandsautomat (i den siste koden jeg postet) som kun benytter seg av
grener og addisjon/subtraksjon, men jeg lurer på om det ikke er en lurere måte å gjøre det på?

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 27, 2019

@terjew

I tillegg unrollet jeg selve telle-loopen slik at den kjører 10 og 10 iterasjoner. Da slipper vi teste i % 10 for hver gjennomkjøring, og kan i tillegg bruke informasjon om hvorvidt tverrsummen er partall til å skippe halvparten av testene for primtall.

Dette har jeg:

const int primes[] = { 0,0,1,1,0,1,0,1,0,0,0,1,0,1,0,0,0,1,0,1,0,0,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0,1,0,0 };

int lastDigitSum;

int inline is_Harshadprime (int num)
{
        int sum;
        if (num % 10 != 0) sum = (++lastDigitSum); // @finnmich sin optimalisering
        else {
                int n = num;
                sum = 0;
                while (n != 0)
                {
                        int rem = n % 10;
                        sum += rem;
                        n /= 10;
                }
                lastDigitSum = sum;
        }
        return  primes[sum] && (num%sum == 0);
}

Jeg tror ikke if (digitsum % 2 != 0) er raskere enn if (!primes[2]) hvis vi antar at primes ligger i CPU-cache under hele kjøringen.
(Men man skal aldri si aldri).

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 28, 2019

@bobcat4

Jeg tror ikke if (digitsum % 2 != 0) er raskere enn if (!primes[2]) hvis vi antar at primes ligger i CPU-cache under hele kjøringen.
(Men man skal aldri si aldri).

Det tror ikke jeg heller, men siden jeg kjører gjennom 10 og 10 tall om gangen trenger jeg bare sjekke i % 2 != 0 én gang for hver 10 tall. Med én sånn test kan jeg utelukke 5 unødvendige inkrementeringer og primtallsjekker for tverrsummer som jeg vet er partall. Dette gir avkastning:

Teste primtall for alle 10 (skippe sjekken for partall): 660ms på stort datasett
Teste primtall bare for oddetall: 590ms på stort datasett

"stort datasett" her er fra 2 til INT_MAX - 256.

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 28, 2019

@DarioSucic @bobcat4

Jeg laget en liknende variant av det @DarioSucic gjorde for å unngå divisjon og mod ved beregning av tverrsummer. Nå har jeg i stedet et array over siffere, og for hver runde i loopen (hver 10. tall) inkrementerer jeg det nest siste sifferet med 1, og håndterer overflowingen i en while-løkke. Dette sparte ytterligere nesten 50% av kjøretiden, nå er jeg nede i 3024ms på 10 kjøringer av stort datasett, altså 303ms per kjøring.

#include <iostream>
#include <thread>
#include <algorithm>
#include <inttypes.h>

bool * GeneratePrimeTruthTable(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 step(uint32_t& digitsum, uint32_t& count, bool* prime, int i)
{
  if (prime[digitsum] && i % digitsum == 0)
  {
    count++;
  }
}

inline int sumDigits(uint32_t i)
{
  uint32_t digitsum = 0;
  uint32_t rest = i;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  rest /= 10;
  digitsum += rest % 10;
  return digitsum;
}

inline void add10(int* digits)
{
  int i = 1;
  while (true) {
    digits[i]++;
    if (digits[i] == 10) {
      digits[i] = 0;
      i++;
    }
    else {
      break;
    }
  }
}

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

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

int SolveRangeOptimized(int start, int end)
{
  auto largestPossibleSum = (int)(round(log10(end) + 1) * 9);
  bool* prime = GeneratePrimeTruthTable(largestPossibleSum, false);
  uint32_t count = 0;
  uint32_t digitsum = 0;
  uint32_t rest = start;

  int blockend = (end / 10) * 10;

  int digits[10];

  int i = start;
  //head:
  do {
    digitsum = sumDigits(i);
    step(digitsum, count, prime, i++);
  } while (i % 10 != 0);

  loadDigits(digits, i - 1);
  digits[0] = 0;
  for (; i < blockend; i += 10)
  {
    add10(digits);
    sumDigits(digits, digitsum);
    step(digitsum, count, prime, i);
    if (digitsum % 2 != 0) {
      digitsum += 2;
      step(digitsum, count, prime, i + 2);
      digitsum += 2;
      step(digitsum, count, prime, i + 4);
      digitsum += 2;
      step(digitsum, count, prime, i + 6);
      digitsum += 2;
      step(digitsum, count, prime, i + 8);
    }
    else {
      digitsum += 1;
      step(digitsum, count, prime, i + 1);
      digitsum += 2;
      step(digitsum, count, prime, i + 3);
      digitsum += 2;
      step(digitsum, count, prime, i + 5);
      digitsum += 2;
      step(digitsum, count, prime, i + 7);
      digitsum += 2;
      step(digitsum, count, prime, i + 9);
    }
  }

  //tail:
  while (i < end) {
    digitsum = sumDigits(i);
    step(digitsum, count, prime, i++);
  }

  delete[] prime;
  return count;
}

int SolveMultithreadedExplicit(int max)
{
  int numthreads = 12;
  int partitionsize = (max - 2) / numthreads;
  int start = 2;
  std::thread * threads = new std::thread[numthreads];
  int * sums = new int[numthreads];

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

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

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

int main()
{
  int count = 10;
  int result = 0;
  for (int i = 0; i < count; i++)
    //result = SolveRangeOptimized(2, 98765432);
    //result = SolveRangeOptimized(2, 1000);
  //result = SolveMultithreadedExplicit(98765432);
    result = SolveMultithreadedExplicit(INT_MAX - 256);
  std::cout << "result is : " << result << "\n";
}
@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 28, 2019

En liten meta-post her for min del.

@terjew @bobcat4 @DarioSucic: Herlig å se dere skvise videre på denne. Har dessverre vært et par dager uten særlig rom for testing her - men ser ikke akkurat ut til at dere har lidd noen nød av den grunn! 🙌

Forøvrig; denne diskusjonen/tråden her er kanskje den rakeste motsetning av den typiske kommentartråden på Internett: Det vises nøkternhet/beskjedenhet rundt egne bidrag, respekt for de andre, og alle påstander backes før eller siden i tester og kvalitative data (AKA kildekritikk/fakta-verifikasjon). Sånn skal det være!

@bobcat4 @DarioSucic
Flott mønstergjenkjenning

@terjew:
Pen unrolling

(Jeg innser ved litt ettertanke at dette ikke er komplimenter jeg får gitt hver dag - dessverre)

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 28, 2019

@michaelo

@terjew @bobcat4 @DarioSucic: Herlig å se dere skvise videre på denne. Har dessverre vært et par dager uten særlig rom for testing her - men ser ikke akkurat ut til at dere har lidd noen nød av den grunn!

Nei,vi har holdt det gående, men nå trenger vi din hjelp!

Forøvrig; denne diskusjonen/tråden her er kanskje den rakeste motsetning av den typiske kommentartråden på Internett: Det vises nøkternhet/beskjedenhet rundt egne bidrag, respekt for de andre, og alle påstander backes før eller siden i tester og kvalitative data (AKA kildekritikk/fakta-verifikasjon). Sånn skal det være!

Enig, umåtelig hyggelig debattklima i disse trådene. Det har vært skikkelig gøy (og utrolig lærerikt) å optimalisere sammen med dere 😄 👍

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 28, 2019

@terjew og @michaelo

Enig, umåtelig hyggelig debattklima i disse trådene. Det har vært skikkelig gøy (og utrolig lærerikt) å optimalisere sammen med dere 😄 👍

Helt enig! Her er det rom for flisespikking (med et smil), øldrikking og koding (med de problemer det måtte medføre), og veldig, veldig lærerikt.

Jeg skal prøve å implementere @DarioSucic sin rekursive digitSum() som jeg tar tro på! 😄

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 28, 2019

Da har jeg følgende C-kode som kjører i én tråd på 170ms (veldig konsekvent etter flere forsøk) på det originale datasettet.
digit_sum_rec() er portet fra @DarioSucic sin Python-kode.

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

const int primes[] = { 0,0,1,1,0,1,0,1,0,0,0,1,0,1,0,0,0,1,0,1,0,0,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0,1,0,0 };
int i;
int cnt;

void digit_sums_rec(int prev, int dindex, int ndigits) {
        if (dindex == ndigits-1) {
                for (int c = 0; c < 10 ; c++) {
                        if (i == 98765433) {
                                printf("Svar: %d\n", cnt);
                                exit(0);
                        }
                        if (primes[prev+c] && (i % (prev+c) == 0)) cnt++;
                        i++;

                }
        }
                else {
                        for (int c = 0; c < 10 ; c++) digit_sums_rec(prev+c, dindex+1, ndigits);
        }
}

int main()
{
        i = cnt = 0;
        digit_sums_rec(0,0,8);
}

perf sier:

            169.95 msec task-clock                #    0.999 CPUs utilized
                 2      context-switches          #    0.012 K/sec
                 0      cpu-migrations            #    0.000 K/sec
                56      page-faults               #    0.330 K/sec
       659,365,257      cycles                    #    3.880 GHz
     1,203,381,625      instructions              #    1.83  insn per cycle
       276,032,967      branches                  # 1624.240 M/sec
         5,480,502      branch-misses             #    1.99% of all branches

       0.170148353 seconds time elapsed

       0.170192000 seconds user
       0.000000000 seconds sys

Den har en "stygg" exit(0) (ja, det gjør litt vondt), men altså den raskeste jeg har skrevet (med god hjelp) så langt i "vanlig" single threaded C

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 28, 2019

@bobcat4

Da har jeg følgende C-kode som kjører i én tråd på 170ms (veldig konsekvent etter flere forsøk):
digit_sum_rec() er portet fra @DarioSucic sin Python-kode.

Pent! Hva med det store datasettet?

Den har en "stygg" exit(0) (ja, det gjør litt vondt), men altså den raskeste så langt i "vanlig" single threaded C

Strengt tatt har du rett i det, men det er ikke så store endringer som må til for å bygge min siste variant som ren C, så jeg vet ikke hvor lenge det vil vare ;)

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 28, 2019

@bobcat4

Single threaded C.
92ms på originalproblemet,
1950ms på utvidet problem

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

#define true 1
#define false 0
#define bool char
#define DIGIT char

bool * generatePrimeTruthTable(uint32_t maxPrime, bool inverted)
{
  uint32_t maxSquareRoot = (uint32_t)sqrt(maxPrime);
  bool* eliminated = (bool*)calloc(maxPrime + 1LL, sizeof(bool));
  eliminated[0] = 1;
  eliminated[1] = 1;
  for (uint32_t j = 4; j <= maxPrime; j += 2)
  {
    eliminated[j] = 1;
  }

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

inline void step(uint32_t * digitsum, uint32_t * count, bool * prime, int i)
{
  if (prime[*digitsum] && i % *digitsum == 0)
  {
    (*count)++;
  }
}

inline void add(DIGIT* digits, int digitno)
{
  digits[digitno]++;
  if (digits[digitno] == 10)
  {
    digits[digitno] = 0;
    add(digits, digitno + 1);
  }
}

inline void add10(DIGIT * digits)
{
  add(digits, 1);
}

inline void add1(DIGIT* digits)
{
  add(digits, 0);
}

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

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

uint32_t solveRange(uint32_t start, uint32_t end)
{
  auto largestPossibleSum = (int)(ceil(log10(end)) * 9);
  bool * prime = generatePrimeTruthTable(largestPossibleSum, 0);
  uint32_t count = 0;
  uint32_t digitsum = 0;
  uint32_t rest = start;

  uint32_t blockend = (end / 10) * 10;

  uint32_t i = start;
  DIGIT digits[10];
  loadDigits(digits, start);

  //head:
  do {
    sumDigits(digits, &digitsum);
    step(&digitsum, &count, prime, i++);
    add1(digits);
  } while (i % 10 != 0);

  for (; i < blockend; i += 10)
  {
    sumDigits(digits, &digitsum);

    step(&digitsum, &count, prime, i);
    if (digitsum % 2 != 0) {
      digitsum += 2;
      step(&digitsum, &count, prime, i + 2);
      digitsum += 2;
      step(&digitsum, &count, prime, i + 4);
      digitsum += 2;
      step(&digitsum, &count, prime, i + 6);
      digitsum += 2;
      step(&digitsum, &count, prime, i + 8);
    }
    else {
      digitsum += 1;
      step(&digitsum, &count, prime, i + 1);
      digitsum += 2;
      step(&digitsum, &count, prime, i + 3);
      digitsum += 2;
      step(&digitsum, &count, prime, i + 5);
      digitsum += 2;
      step(&digitsum, &count, prime, i + 7);
      digitsum += 2;
      step(&digitsum, &count, prime, i + 9);
    }
    add10(digits);
  }

  loadDigits(digits, i);
  //tail:
  while (i < end) {
    sumDigits(digits, &digitsum);
    step(&digitsum, &count, prime, i++);
    add1(digits);
  }

  free(prime);
  return count;
}

int main()
{
  int count = 10;
  uint32_t result = 0;
  for (int i = 0; i < count; i++)
  {
    //result = solveRange(2, 98765432);
    result = solveRange(2, INT_MAX - 256);
  }
  printf("result is : %d\n", result);
}
@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 29, 2019

@bobcat4

Det pussige er at C++-versjonen kjører raskere enn C-versjonen enkelttrådet hos meg. Kanskje noe med at Microsoft-kompilatoren er mer aktivt vedlikeholdt for C++-kode enn for C?

Single threaded C++:
Opprinnelig datasett: 85ms
Utvidet datasett: 1765ms

@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 29, 2019

@terjew

Det pussige er at C++-versjonen kjører raskere enn C-versjonen enkelttrådet hos meg. Kanskje noe med at Microsoft-kompilatoren er mer aktivt vedlikeholdt for C++-kode enn for C?

Det vil jo være interessant å se hvordan dette er sammenliknet med f.eks. GCC eller clang. Jeg kan ikke se noen av de "typiske" tingene som skal hjelpe C++-kompilatoren med optimalisering vs C's heller.

Har du sett på den genererte assemblyen?

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 29, 2019

@michaelo

Det pussige er at C++-versjonen kjører raskere enn C-versjonen enkelttrådet hos meg. Kanskje noe med at Microsoft-kompilatoren er mer aktivt vedlikeholdt for C++-kode enn for C?

Det vil jo være interessant å se hvordan dette er sammenliknet med f.eks. GCC eller clang. Jeg kan ikke se noen av de "typiske" tingene som skal hjelpe C++-kompilatoren med optimalisering vs C's heller.

Har ikke clang eller gcc på maskinen min så langt, det så ut til å være litt mas å få det opp å kjøre på windows. Du kan gjerne kjøre en bench hos deg med de 2 siste versjonene jeg postet (hhv C og C++).

Har du sett på den genererte assemblyen?

Nei, hender jeg tar en kikk på godbolt om det er noe jeg lurer på, men jeg bruker ikke så mye tid på å inspisere assemblykode. For å være helt ærlig så er det ganske lenge siden jeg skrev noe assembly (20 år minst), og da var det for enklere plattformer. Dermed føler jeg dessverre at jeg får lite ut av å kikke på moderne X86/X64-assembly.

@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 29, 2019

@terjew

Det pussige er at C++-versjonen kjører raskere enn C-versjonen enkelttrådet hos meg. Kanskje noe med at Microsoft-kompilatoren er mer aktivt vedlikeholdt for C++-kode enn for C?

Det vil jo være interessant å se hvordan dette er sammenliknet med f.eks. GCC eller clang. Jeg kan ikke se noen av de "typiske" tingene som skal hjelpe C++-kompilatoren med optimalisering vs C's heller.

Har ikke clang eller gcc på maskinen min så langt, det så ut til å være litt mas å få det opp å kjøre på windows. Du kan gjerne kjøre en bench hos deg med de 2 siste versjonene jeg postet (hhv C og C++).

Jepp! Var mest ledende tanker for min egen del. Blir endelig litt rom i kveld igjen for å leke med slikt.

Har du sett på den genererte assemblyen?

Nei, hender jeg tar en kikk på godbolt om det er noe jeg lurer på, men jeg bruker ikke så mye tid på å inspisere assemblykode. For å være helt ærlig så er det ganske lenge siden jeg skrev noe assembly (20 år minst), og da var det for enklere plattformer. Dermed føler jeg dessverre at jeg får lite ut av å kikke på moderne X86/X64-assembly.

Hehe. Ser den! Er ingen proff på det feltet selv, men funnet det spesielt lærerikt ganske enkelt å diffe ulike varianter for f.eks. å se om den har klart å inline som forventet. Godbolt tjener jo forøvrig fort samme nytte 👍

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 29, 2019

@michaelo @bobcat4 @DarioSucic
Ryddet litt i datatyper og forsøkte med forskjellige størrelser på lagring av hhv tverrsum, siffere, tellere etc. Ser ut som best ytelse er med uint32_t på det meste, med unntak av selve sifferene som pussig nok gir best ytelse for uint64_t. Kanskje noe spes med at int-addering er hakket raskere i 64 bit på en 64-bit cpu?

Det ser ut som jeg er helt på grensen av L1 cache eller registerbruk, for selv de minste, ubetydelige endringer (bytte rekkefølge på variabler etc) kan fort øke kjøretiden betraktelig.

Uansett, med de siste små endringene kommer jeg meg ned på 14.5ms på det opprinnelige problemet og 299ms på det store problemet (riktignok uten å kjøre Brave i bakgrunnen :P). Det skulle vel tilsi ca 3 ganger kjøretiden til CUDA-løsningen?

#include <iostream>
#include <thread>
#include <algorithm>
#include <inttypes.h>

typedef uint64_t digit_t;

bool* GeneratePrimeTruthTable(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)
{
  if (prime[digitsum] && i % digitsum == 0)
  {
    count++;
  }
}

inline void add10(digit_t* digits)
{
  uint32_t i = 1;
  while (true) {
    digits[i]++;
    if (digits[i] == 10) {
      digits[i] = 0;
      i++;
    }
    else {
      break;
    }
  }
}

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


uint32_t SolveRangeOptimized(uint32_t start, uint32_t end)
{
  bool* prime = GeneratePrimeTruthTable((int)(round(log10(end) + 1) * 9), false);

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

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

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

  delete[] prime;
  return count;
}


int SolveMultithreadedExplicit(uint32_t start, uint32_t max)
{
  unsigned concurentThreadsSupported = std::thread::hardware_concurrency();
  uint32_t numthreads = concurentThreadsSupported;
  uint32_t partitionsize = (max - 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, max);
    uint32_t currentThread = threadNo;
    threads[threadNo] = std::thread(
      [threadBegin, threadEnd, currentThread, &sums] {
        sums[currentThread] = SolveRangeOptimized(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;
}

int main()
{
  int count = 10;
  int result = 0;
  for (int i = 0; i < count; i++)
    //result = SolveRangeOptimized(2, 98765432);
    //result = SolveRangeOptimized(2, INT_MAX - 256);
    //result = SolveMultithreadedExplicit(98765432);
    result = SolveMultithreadedExplicit(2, INT_MAX - 256);
  std::cout << result << "\n";
}
@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 29, 2019

Er det forresten flere enn meg som begynner å kjenne litt på likheten til denne XKCD-stripen?

@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 29, 2019

FWIW så er dette målingene (for 10x kjøringer -> justert) på to siste løsningene som postet av @terjew over her (C og C++) på min MBP 2018 Intel i7-7700HQ 2.8GHz med clang v11.0.0:

c-solution (st):
small: 1.899s -> 189.9ms
big: 39.969s -> 3996.9ms

cpp-solution (st):
small: 0.764s -> 76.5ms
big: 15.789s -> 1578.9ms

cpp-solution (mt):
small: 0.200s -> 20.0ms
big: 4.020s -> 402.0ms

Obs, har ikke fått sett på hits/misses osv her enda.

@terjew

Er det forresten flere enn meg som begynner å kjenne litt på likheten til denne XKCD-stripen?

👏 Man kan geeke hardt og langt ut om hva som helst ja. Men jeg vil gjerne argumentere for at forståelse rundt HW og effektivitet kan ha en ikke uvesentlig verdi (og samfunnsverdi sett i det store bildet) all den tid det er så mange som ikke bryr seg om det.

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 29, 2019

@michaelo
Sant det du sier om nytten av dette, jeg hadde nok ikke giddet å dra det så langt om jeg følte det var helt bortkastet 😁

For øvrig, jeg installerte Ubuntu gjennom WSL på Windows 10 fordi grunner.

Det morsomme er at den samme c++-koden som kjører på 299ms i Windows kjører på 279ms i WSL når jeg bygger med g++ og -O3 🤣

wiesener@BOBTOP:~$ g++ -pthread -O3 -o knowit23 knowit23.cpp
wiesener@BOBTOP:~$ time ./knowit23
13650795

real    0m2.785s
user    0m31.813s
sys     0m0.016s
wiesener@BOBTOP:~$                                                         

Og med clang:

wiesener@BOBTOP:~$ clang++ -pthreads -o knowit23-clang -O3 knowit23.cpp
wiesener@BOBTOP:~$ time ./knowit23-clang
13650795

real    0m3.185s
user    0m36.734s
sys     0m0.000s
wiesener@BOBTOP:~$                                                                                                                                                                                                  
@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 29, 2019

@terjew

Jeg var veldig nysgjerrig på ytelsesforskjellen mellom C og C++ og har foreløpig funnet følgende:

Ved å endre type på DIGIT fra char (slik den var i din C-versjon) til uint64_t (slik digit_t er i din C++-versjon) så reduserer jeg kjøretid fra 1.9s til 0.8s. Det gir forsåvidt god mening iom at CPUen ikke har noen kjennskap tiloperasjoner for typer under sin native vidde og dermed må konvertere. Da er vi kun ~40ms/4ms (for 10x/1x kjøring) i forskjell mellom de to variantene.

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 29, 2019

@michaelo

Jeg var veldig nysgjerrig på ytelsesforskjellen mellom C og C++ og har foreløpig funnet følgende:

Ved å endre type på DIGIT fra char (slik den var i din C-versjon) til uint64_t (slik digit_t er i din C++-versjon) så reduserer jeg kjøretid fra 1.9s til 0.8s. Det gir forsåvidt god mening iom at CPUen ikke har noe kjennskap til typer under sin native vidde og dermed må konvertere. Da er vi kun ~40ms/4ms (for 10x/1x kjøring) i forskjell mellom de to variantene.

Ah, se der ja. Det hadde jeg ikke tenkt over engang, til tross for at jeg så at å endre fra 32 til 64-bit i c++ gjorde stort utslag... Bra sett!

Hos meg havner da C-versjonen bygget med gcc på 71ms for lite sett og 1503ms for stort.
Med clang havner jeg på 68ms på lite og 1427ms på stort.

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 29, 2019

@michaelo

Ved å endre type på DIGIT fra char (slik den var i din C-versjon) til uint64_t (slik digit_t er i din C++-versjon) så reduserer jeg kjøretid fra 1.9s til 0.8s. Det gir forsåvidt god mening iom at CPUen ikke har noe kjennskap til typer under sin native vidde og dermed må konvertere. Da er vi kun ~40ms/4ms (for 10x/1x kjøring) i forskjell mellom de to variantene.

Vent litt, redusert til 0.8s på stort datasett enkelttrådet? Det er jo en forbedring på nesten 60%, og langt forbi de andre versjonene... er du sikker på at det tallet stemmer?

@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 29, 2019

@terjew

Vent litt, redusert til 0.8s på stort datasett enkelttrådet? Det er jo en forbedring på nesten 60%, og langt forbi de andre versjonene... er du sikker på at det tallet stemmer?

Dette var målingsforskjellen på lite datasett, med 10x loop.
Tilsvarende for stort sett er hopp fra ~40s til ~16.6s

@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 29, 2019

Så /10 så harmoniserer det godt med det du ser :)

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 29, 2019

@terjew

Vent litt, redusert til 0.8s på stort datasett enkelttrådet? Det er jo en forbedring på nesten 60%, og langt forbi de andre versjonene... er du sikker på at det tallet stemmer?

Dette var målingsforskjellen på lite datasett, med 10x loop.
Tilsvarende for stort sett er hopp fra ~40s til ~16.6s

Hmm, så dramatiske forskjeller observerer jeg ikke her. Noe forverring ved å gå tilbake til char men ikke all verden.

clang 6.0.0:
uint64_t: 14.273s
char: 15.087s

gcc 7.4.0:
uint64_t: 15.030s
char: 16.273s

cl 19.24.28314:
uint64_t: 17.680s
char: 19.631s
@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 30, 2019

@terjew:
Humm. Jeg feil-leste nok posten din lenger opp ja med C-målingene ja. Testene/målingene er veldig gjenskapbare hos meg ihvertfall. Ganske brutal påvirkning! Men - når alt annet er optimalisert så blir jo også de "mindre" kostnadene store relativt sett.

Men du fikk da C-versjonen ned til C++-versjonen?

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 30, 2019

@michaelo

Man kan geeke hardt og langt ut om hva som helst ja. Men jeg vil gjerne argumentere for at forståelse rundt HW og effektivitet kan ha en ikke uvesentlig verdi (og samfunnsverdi sett i det store bildet) all den tid det er så mange som ikke bryr seg om det.

Enig der, sånn rent akademisk er dette gull, selv om dette kanskje begynner å miste litt "reell verdi" ettersom det er snakk om veldig spesifikk HW.

Apropos: Jeg gjetter denne oppgaven løses raskest med FPGA, men da er vi virkelig langt unna "allfarvei", og det blir også en kostbar løsning. 😄

@michaelo

This comment has been minimized.

Copy link

@michaelo michaelo commented Dec 30, 2019

@bobcat4

Enig der, sånn rent akademisk er dette gull, selv om dette kanskje begynner å miste litt "reell verdi" ettersom det er snakk om veldig spesifikk HW.

Det er nok ikke så mange arbeids-/oppdragsgivere vi kunne fakturert for dette nei 😁

Apropos: Jeg gjetter denne oppgaven løses raskest med FPGA, men da er vi virkelig langt unna "allfarvei", og det blir også en kostbar løsning. 😄

Altså. Hvis vi går den veien så kan vi sikkert få laget noe i dedikert HW også. Bli kvitt all den ugreie koden. Kanskje tom få en egen ALU i silikonet for Harshadprimtall? 🤔

Men altså: Det hadde vært gøy med noen VHDL/FPGA-løsninger på oppgavene her også.

@terjew

This comment has been minimized.

Copy link

@terjew terjew commented Dec 30, 2019

@bobcat4 @michaelo
Mye gode ideer her :) usikker på betalingsvilje for hardwareløsning på dette problemet... Men jeg har fortsatt noen ideer som jeg vil teste ut, tror softwareløsningen fortsatt har en del rom for forbedring :)

@bobcat4

This comment has been minimized.

Copy link

@bobcat4 bobcat4 commented Dec 30, 2019

@michaelo

Altså. Hvis vi går den veien så kan vi sikkert få laget noe i dedikert HW også. Bli kvitt all den ugreie koden. Kanskje tom få en egen ALU i silikonet for Harshadprimtall? 🤔

LOL!

Men altså: Det hadde vært gøy med noen VHDL/FPGA-løsninger på oppgavene her også.

Ja, jeg tenkte det samme!

Det er kanskje ikke så veldig mange som skriver VHDL. Jeg har bare s