Project Euler - Aufgabe 586 (Wettbewerb)

[GPI]

oh verdammt. ja auf r=40.....

[NicTheQuick]

Ich hab mal so nebenbei Statistik geführt für alle n = 1.000.000.

Code: Alles auswählen

f(n, 1) = 87505
f(n, 2) = 48727
f(n, 3) =  2149
f(n, 4) =  5306
f(n, 5) =     6
f(n, 6) =   230
f(n, 7) =     0
f(n, 8) =    55

Für r > 8 war f(n, r) bisher immer Null.
Die Statistik hat jetzt 3 Minuten und 39,571 Sekunden gedauert, weil ich wieder so ran gegangen bin wie GPI am Anfang. Ich hab es nur etwas optimiert.

Ich hab natürlich nicht nur die Ergebnisse gezählt, sondern auch für jedes r die zugehörigen ks gespeichert, weil ich das ganze in einen Graph plotten wollte um vielleicht etwas erkennen zu können, was weiter hilft. Leider sieht man nicht viel. Ich hab's mal hier in meine Owncloud gestellt: Euler

Den Code kann man eigentlich sehr gut in Threads auslagern, aber wenn ich für 1 Mio Tests schon 3 Minuten brauche, dann dauert es für 10^15 wohl etwas zu lang.

[NicTheQuick]

Okay, n = 1.000.000 schaffe ich jetzt auch in 2 Sekunden.
Und n = 10.000.000 dauert 77 Sekunden. Alles nur ein Thread. Ich muss da mal Threads einbauen. Außerdem hab ich an einer Stelle immer noch eine Wurzel-Funktion. Das ist noch etwas, was ich unbedingt ändern muss. Ich hab aber noch andere Ideen.

Hier mal die Stats:

Code: Alles auswählen

f(10000000, 1) = 724319
f(10000000, 2) = 496337
f(10000000, 3) = 18750
f(10000000, 4) = 81624
f(10000000, 5) = 108
f(10000000, 6) = 4472
f(10000000, 7) = 1
f(10000000, 8) = 2676
f(10000000, 9) = 53
f(10000000, 10) = 2
f(10000000, 12) = 59
f(10000000, 16) = 1

Diese Version der Herangehensweise verbraucht keinen Speicherplatz und könnte im Grunde durchgehend komplett im Cache der CPU ausgeführt werden, wenn das Betriebssystem nicht mitzureden hätte. Die Laufzeitanalyse gestaltet sich etwas schwierig, ist aber sicherlich größer als O(n). Dennoch skaliert es schöner als bei euren letzten Ansätzen.

Hier mal mein Code:

Code: Alles auswählen

EnableExplicit

Macro func(a, b)
	((a) * (a) + 3 * (a) * (b) + (b) * (b))
EndMacro

; Wenn #True, werden für alle 0 < r <= #maxR alle zutreffenden Zahlen gespeichert und am Ende gelistet
; VORSICHT: Für hohe n könnte sehr viel Speicher beansprucht werden!
#doFullStatistics = #False

; Wenn #True, wird der Algorithmus automatisch für alle 0 < r < #maxR ausgeführt
#doSimpleStatistics = #True

; Obergrenze für r für die Statistik
#maxR = 100

; Millisekunden, die vergehen müssen, bevor ein Statusupdate ausgegeben wird.
; 0 zum Deaktivieren
#progress = 1000

CompilerIf #doFullStatistics
	Structure Sequence
		List number.i()
	EndStructure
	Dim Sequences.Sequence(#maxR)
CompilerEndIf

CompilerIf #doSimpleStatistics
	Define innerLoops.q = 0
	Dim rStat.i(#maxR)
CompilerEndIf

OpenConsole()

Define n.q, r.i

n = 10000000
r = 6

; Zeit
Define time.i, nextPrintTime.i
time = ElapsedMilliseconds()
nextPrintTime = time + #progress

; Ergebnis für f(n, r)
Define.i result = 0

; interne Variablen
Define k.q, maxA.q, maxAIncAt.q, minA.q, minAIncAt.q
k = 11
maxA = 2                         ; Obergrenze für a für aktuelles k
maxAIncAt = func(maxA, 1) + 1    ; das k, ab dem maxA inkrementiert wird
minA = 2                         ; Untergrenz für a für aktuelles k
minAIncAt = func(minA + 1, minA) ; das k, ab dem minA inkrementiert wird

Define matches.i, a.q, a3.q, b.q, minB.q, kA.q
While k <= n
	matches = 0
	
	; Falls a² + 3a + 2 = k mit a = maxA, erhöhe maxA
	If maxAIncAt = k
		maxA + 1
		; Anstatt maxAIncAt immer komplett neu auszurechnen, addiere einfach die Differenz
		maxAIncAt + maxA + maxA + 2
	EndIf
	
	; Falls (a+1)² + 3(a+1)a + a² = k mit a = minA, erhöhe minA
	If minAIncAt = k
		; Anstatt minAIncAt immer komplett neu auszurechnen, addiere einfach die Differenz
		minAIncAt + 10 * minA
		minA + 1
	EndIf
	
	;PrintN("k=" + k + ", a=" + minA + ".." + maxA)
	a = maxA
	a3 = 3 * a
	
	; Berechne das kleinste b, sodass a² + 3ab + b² <= k ist
	;TODO Weg damit!!!
	minB = IntQ(Sqr(5 * a * a + 4 * k) - 3 * a) / 2
	b = minB
	
	; Berechne den Initialwert einmalig pro k
	;TODO Ist nicht so schlimm wie das Sqr() bzw. allgemein Floats
	kA = func(a, b)
	
	Repeat
		CompilerIf #doSimpleStatistics
			innerLoops + 1
		CompilerEndIf
		
		; Haben wir einen Treffen, zähle ihn, und fahre fort mit dem nächst kleineren a
		If kA = k
			matches + 1
			a - 1
			If a < minA
				Break
			EndIf
			a3 - 3
			kA - (a + a + 1 + 3 * b)
		
		; Ist a² + 3ab + b² zu groß geworden, dekrementiere a
		ElseIf kA > k
			;
			a - 1
			If a < minA
				Break
			EndIf
			a3 - 3
			kA - (a + a + 1 + 3 * b)
		
		; Ist a² + 3ab + b² zu klein geworden, inkrementiere b
		ElseIf kA < k
			kA + a3 + b + b + 1
			b + 1
		EndIf
		
		; Sind a = b geworden, gilt die Regeln a > b > 0 nicht mehr, also beenden wir die Schleife
	Until a = b
	
	CompilerIf #doFullStatistics Or #doSimpleStatistics
		If matches > 0 And matches <= #maxR
			CompilerIf #doFullStatistics
				If AddElement(Sequences(matches)\number())
					Sequences(matches)\number() = k
				EndIf
			CompilerEndIf
			CompilerIf #doSimpleStatistics
				rStat(matches) + 1
			CompilerEndIf
		EndIf
	CompilerEndIf
	
	If matches = r
		result + 1
	EndIf
	
	CompilerIf #progress > 0
		If ElapsedMilliseconds() > nextPrintTime
			nextPrintTime + #progress
			PrintN("k=" + k + ", a=" + minA + ".." + maxA + ", f(n, r) >= " + result)
		EndIf
	CompilerEndIf
	
	k + 1
Wend

time = ElapsedMilliseconds() - time

PrintN(~"\nResult:\nn = " + n + ~"\nr = " + r + ~"\nf(n, r) = " + result + ~"\nTime: " + time + ~" ms\n")

CompilerIf #doSimpleStatistics
	Print("Enter for simple statistics...")
	Input()
	PrintN("Inner loops: " + innerLoops + " (avg. per k: " + Str(innerLoops / n) + ")")
	For r = 1 To #maxR
		If rStat(r)
			PrintN("f(" + n + ", " + r + ") = " + rStat(r))
		EndIf
	Next
CompilerEndIf

CompilerIf #doFullStatistics	
	Print("Enter for full statistics...")
	Input()
	For r = 1 To #maxR
		With Sequences(r)
			If ListSize(\number()) > 0
				Print("r=" + r + ":")
				ForEach \number()
					Print(" " + \number())
				Next
				PrintN("")
			EndIf
		EndWith
	Next
CompilerEndIf

Print(~"\nEnter to quit...")
Input()
CloseConsole()

Update:
Hm... Mit 8 Threads schaffe ich es in 22 Sekunden. Ich dachte da kommt mehr bei rum. Aber optimal ist es noch nicht. Im Grunde lasse ich jeden Thread bei einem anderen k starten und inkrementiere das k dann immer um die Anzahl der Threads. So kann es jetzt natürlich passieren, dass am Ende ein Thread fertig ist und die andere noch eine Weile rödeln. Ein Threadpool wäre hier praktischer, aber komplexer zu realisieren. Den Code kann ich vielleicht später posten. Hab ihn grad noch verschlimmbessert.

Update:

Code: Alles auswählen

EnableExplicit

Macro func(a, b)
	((a) * (a) + 3 * (a) * (b) + (b) * (b))
EndMacro

; Wenn #True, wird der Algorithmus automatisch für alle 0 < r < #maxR ausgeführt
#doSimpleStatistics = #True

; Obergrenze für r für die Statistik
#maxR = 100

; Millisekunden, die vergehen müssen, bevor ein Statusupdate ausgegeben wird.
; 0 zum Deaktivieren
#progress = 1000

CompilerIf #doSimpleStatistics
	Dim rStat.i(#maxR)
CompilerEndIf

CompilerIf #progress > 0
	Define lastInnerLoops.q = 0
CompilerEndIf

CompilerIf #progress > 0 Or #doSimpleStatistics
	Define innerLoops.q = 0
CompilerEndIf

OpenConsole()

Define n.q, r.i

#threads = 8

n = 10000000
r = 6

; Zeit
Define startTime.i, nextPrintTime.i
startTime = ElapsedMilliseconds()
nextPrintTime = startTime + #progress


Structure ThreadData
	result.i
	kOffset.q
	kEnd.q
	skip.i
	r.i
	k.q
	
	mutex.i
	refreshStat.i
	started.i
	
	
	CompilerIf #doSimpleStatistics
		rStat.i[#maxR + 1]
	CompilerEndIf
	
	CompilerIf #progress > 0 Or #doSimpleStatistics
		innerLoops.q
	CompilerEndIf
EndStructure

; Ergebnis für f(n, r)
Define.i result = 0

Procedure Thread(*threadData.ThreadData)

	With *threadData
		LockMutex(\mutex)
		PrintN("Thread started. kOffset=" + \kOffset)
		UnlockMutex(\mutex)
		
		SignalSemaphore(\started)
		
		\result = 0
		
		CompilerIf #doSimpleStatistics
			Protected i.i
			For i = 0 To #maxR
				\rStat[i] = 0
			Next
		CompilerEndIf
		
		CompilerIf #progress > 0 Or #doSimpleStatistics
			\innerLoops = 0
		CompilerEndIf
	
		; interne Variablen
		Protected k.q, maxA.q, maxAIncAt.q, minA.q, minAIncAt.q
		k = 11 + \kOffset
		
		;TODO maxA und minA berechnen
		maxA = 2                         ; Obergrenze für a für aktuelles k
		maxAIncAt = func(maxA, 1) + 1    ; das k, ab dem maxA inkrementiert wird
		minA = 2                         ; Untergrenz für a für aktuelles k
		minAIncAt = func(minA + 1, minA) ; das k, ab dem minA inkrementiert wird
		
		Protected matches.i, a.q, a3.q, b.q, minB.q, kA.q
		While k <= \kEnd
			\k = k
			matches = 0
			
			; Falls a² + 3a + 2 = k mit a = maxA, erhöhe maxA
			While maxAIncAt <= k
				maxA + 1
				; Anstatt maxAIncAt immer komplett neu auszurechnen, addiere einfach die Differenz
				maxAIncAt + maxA + maxA + 2
			Wend
			
			; Falls (a+1)² + 3(a+1)a + a² = k mit a = minA, erhöhe minA
			While minAIncAt <= k
				; Anstatt minAIncAt immer komplett neu auszurechnen, addiere einfach die Differenz
				minAIncAt + 10 * minA
				minA + 1
			Wend
			
			;PrintN("k=" + k + ", a=" + minA + ".." + maxA)
			a = maxA
			a3 = 3 * a
			
			; Berechne das kleinste b, sodass a² + 3ab + b² <= k ist
			;TODO Weg damit!!!
			minB = IntQ(Sqr(5 * a * a + 4 * k) - 3 * a) / 2
			b = minB
			
			; Berechne den Initialwert einmalig pro k
			;TODO Ist nicht so schlimm wie das Sqr() bzw. allgemein Floats
			kA = func(a, b)
			
			Repeat
				CompilerIf #doSimpleStatistics Or #progress > 0
					\innerLoops + 1
				CompilerEndIf
				
				; Haben wir einen Treffen, zähle ihn, und fahre fort mit dem nächst kleineren a
				If kA = k
					matches + 1
					a - 1
					If a < minA
						Break
					EndIf
					a3 - 3
					kA - (a + a + 1 + 3 * b)
				
				; Ist a² + 3ab + b² zu groß geworden, dekrementiere a
				ElseIf kA > k
					;
					a - 1
					If a < minA
						Break
					EndIf
					a3 - 3
					kA - (a + a + 1 + 3 * b)
				
				; Ist a² + 3ab + b² zu klein geworden, inkrementiere b
				ElseIf kA < k
					kA + a3 + b + b + 1
					b + 1
				EndIf
				
				; Sind a = b geworden, gilt die Regeln a > b > 0 nicht mehr, also beenden wir die Schleife
			Until a = b
			
			CompilerIf #doSimpleStatistics
				If matches > 0 And matches <= #maxR
					\rStat[matches] + 1
				EndIf
			CompilerEndIf
			
			If matches = \r
				\result + 1
			EndIf
			
			k + \skip
			
		Wend
	
		LockMutex(\mutex)
		PrintN("Thread done. kOffset=" + \kOffset)
		UnlockMutex(\mutex)
	
	EndWith
EndProcedure

Structure Thread
	handle.i
	threadData.ThreadData
EndStructure

Dim threads.Thread(#threads - 1)
Define mutex.i = CreateMutex()

; Initialisiere ThreadData
Define i.i
For i = 0 To #threads - 1
	With threads(i)\threadData
		\kOffset = i
		\kEnd = n
		\skip = #threads
		\r = r
		\mutex = mutex
		\refreshStat = #False
		\started = CreateSemaphore(0)
	EndWith
Next

; Starte Threads
For i = 0 To #threads - 1
	With threads(i)
		\handle = CreateThread(@Thread(), @\threadData)
	EndWith
Next

; Warte bis alle laufen
For i = 0 To #threads - 1
	With threads(i)\threadData
		WaitSemaphore(\started)
	EndWith
Next

CompilerIf #progress > 0
	Define kLast.q = 0, threadsRunning.i
	
	Repeat
		Define currentTime = ElapsedMilliseconds()
		threadsRunning.i = #threads
		If currentTime > nextPrintTime
			nextPrintTime = currentTime
			
			result = 0
			Define.i kMin = n
			Define.i kMax = 0
			Define.q innerLoops = 0
			
			For i = 0 To #threads - 1
				With threads(i)\threadData
					result + \result
					innerLoops + \innerLoops
					If \k < kMin
						kMin = \k
					EndIf
					If \k > kMax
						kMax = \k
					EndIf
				EndWith
				If Not IsThread(threads(i)\handle)
					threadsRunning - 1
				EndIf
			Next
			
			PrintN("k=" + kMin + ".." + kMax +
			       ", f(n, r) >= " + result +
			       ", k/sec=" + StrD((kMin - kLast) * 1000 / #progress, 2) +
			       ", loops/sec=" + StrD((innerLoops - lastInnerLoops) * 1000 / #progress, 2) +
			       ", time=" + Str(currentTime - startTime))
			kLast = kMin
			lastInnerLoops = innerLoops
		EndIf
		Delay(#progress)
	Until threadsRunning = 0
CompilerElse

	; Warte auf Thread-Ende
	For i = 0 To #threads - 1
		With threads(i)
			WaitThread(\handle)
		EndWith
	Next
CompilerEndIf

; Akkumuliere Statistik
CompilerIf #doSimpleStatistics
	Define j.i	
	For i = 0 To #threads - 1
		With threads(i)\threadData
			For j = 1 To #maxR
				rStat(j) + \rStat[j]
			Next
		EndWith
	Next
CompilerEndIf

result = 0
For i = 0 To #threads - 1
	With threads(i)\threadData
		result + \result
	EndWith
Next

startTime = ElapsedMilliseconds() - startTime

PrintN(~"\nResult:\nn = " + n + ~"\nr = " + r + ~"\nf(n, r) = " + result + ~"\nTime: " + startTime + ~" ms\n")

CompilerIf #doSimpleStatistics
	Print("Enter for simple statistics...")
	Input()
	PrintN("Inner loops: " + innerLoops + " (avg. per k: " + Str(innerLoops / n) + ")")
	For r = 1 To #maxR
		If rStat(r)
			PrintN("f(" + n + ", " + r + ") = " + rStat(r))
		EndIf
	Next
CompilerEndIf

Print(~"\nEnter to quit...")
Input()
CloseConsole()

Ergebnis:

Code: Alles auswählen

Result:
n = 10000000
r = 6
f(n, r) = 4472
Time: 24207 ms

Enter for simple statistics...
Inner loops: 21066879501 (avg. per k: 2106)
f(10000000, 1) = 724319
f(10000000, 2) = 496337
f(10000000, 3) = 18750
f(10000000, 4) = 81624
f(10000000, 5) = 108
f(10000000, 6) = 4472
f(10000000, 7) = 1
f(10000000, 8) = 2676
f(10000000, 9) = 53
f(10000000, 10) = 2
f(10000000, 12) = 59
f(10000000, 16) = 1

[GPI]

Ich frag mich wirklich, was der Trick ist....

Meine letzte Idee war, irgendwie die Ks so zu berechnen, das sie Sortiert rauskommen.... Hat nicht wirklich geklappt, in niedrigen Bereich ist es noch einfach, aber sobald es in die Tausende geht, keine Chance...

Ich hab meine letzte Version ein bischen angepasst, das die Vorhersage der Zeit besser wird. Früher hab ich einfach mit a geschätzt. Das ist extrem ungenau, da ja für jedes a a-1 b gibt. Das wächst also gewaltig. Jetzt werte ich die tatsächliche Ks die berechnet werden. Ist ziemlich genau jetzt. Ergebnis: ca. 1300 Stunden!

Es muss definitv ein Trick geben. Die Aufgabe ist von 19.01 - die haben keine 50 Tage dran rechnen können oder haben zugriff auf einen größeren Server (glaub ich ehrlich gesagt weniger).

A*B zu berechnen ist definitv eine Sackgasse. Das sind 215.204.447.601.369 Berechnungen -> ohne Auswertung!

[helpy]

Für eine ganz andere Lösung braucht man wohl Zahlentheorie!

==> https://de.wikipedia.org/wiki/Bin%C3%A4 ... zer_Zahlen
==> https://de.wikipedia.org/wiki/Bin%C3%A4 ... orm#Nutzen

Soweit ich das verstanden hat, braucht man zur Lösung "diophantische Gleichungen" ... und da habe ich nicht die Zeit mich da hineinzudenken!

Aber vielleicht bin ich da auf einer ganz falschen Spur.

[NicTheQuick]

Ich hab mir schon gedacht, dass da ein bisschen mehr Zahlentheorie dahinter stecken müsste, dennoch versucht man bei einer Programmieraufgabe erst mal programmiertechnisch ran zu gehen. Wenn ich wieder Zeit und Lust habe, gehe ich mal näher in die Zahlentheorie.

[GPI]

Mein Ansatz ist mittlerweile folgender:
In einer Zeile ist a immer konstant und b wird einfach nur größerer. Innerhalb einer Zeile kann es keine doppelten Ergebnisse geben, die Formel enthält nur Additionen
Es gilt als mehrere Zeilen zu untersuchen, ob sie sich "schneiden".
Hinzu kommt, das wir nur die Ks überprüfen müssen, die zwischen endzeile minimum-k und startzeile maximum-k.
Damit könnte man gewaltig viele Berechnungen ignorieren.

Ich schau gerade, wie ich das ganze in einem Programm umsetze. Entweder hab ich ein Denkfehler in meiner Theorie oben oder ich programmiere gerade Mist

[NicTheQuick]

So mache ich es eigentlich schon, wenn ich dich nicht ganz missverstanden habe.

Sei g(a, b) = a² + 3ab + b², dann suche ich aMin <= a <= aMax, sodass g(a, 1) <= k <= g(a, a - 1).

Ich fange bei k = 11 an. Hier ist aMin = 2 und aMax = 2, denn dann kann a nur 2 sein, und in dem Fall ist g(a, 1) = g(a, a - 1) = g(2, 1).
Im nächsten Schritt ist k = 12. Jetzt muss ich aMax = 3 setzen, damit die rote Formel oben noch zutrifft. Irgendwann kann ich auch aMin erhöhen, da dann die rote Formel immer noch zutrifft.

Ich teste also für jedes k dann alle a im Intervall [aMin, aMax]. Und da dann in der Gleichung g(a, b) = k nur noch eine Unbekannte, nämlich b, über bleibt, berechne ich das nach der guten alten Schulmathematik. Das Problem ist, dass sich das Intervall mit wachsendem k auch vergrößert, zwar zum Glück unter der Wurzel, aber immerhin. So sieht der Wachstum der Intervallgröße aus: plot 0.5*(sqrt(4*k+5)-3) and 0.1*(sqrt(20*k+5)+5) for k=0..1000

[NicTheQuick]

Ich hab meinen Code mal mit 8 Threads auf einem richtigen Server laufen lassen mit n = 100.000.000. Der geht ganz schön ab. Ist ein Intel(R) Xeon(R) CPU E3-1230 v3 @ 3.30GHz.
Er hat gerade mal 9 Minuten und 40 Sekunden gebraucht.

Code: Alles auswählen

Inner loops: 666516755726 (avg. per k: 6665)
f(100000000, 1) = 6164497
f(100000000, 2) = 4877564
f(100000000, 3) = 155038
f(100000000, 4) = 1051883
f(100000000, 5) = 1074
f(100000000, 6) = 59517
f(100000000, 7) = 25
f(100000000, 8) = 63124
f(100000000, 9) = 970
f(100000000, 10) = 117
f(100000000, 12) = 2936
f(100000000, 13) = 4
f(100000000, 16) = 492
f(100000000, 18) = 13
f(100000000, 24) = 1

Mein virtueller Server mit 6 VCores läuft schon eine Stunde und ist damit noch nicht fertig.
Update: Er hat 69 Minuten und 52 Sekunden gebraucht. Verrückt.

Cool, jetzt hab ich einen Benchmark.