Project Euler - Aufgabe 586 (Wettbewerb)

[GPI]

https://projecteuler.net/problem=586

The number 209 can be expressed as a^2+3ab+b^2 in two distinct ways:
209=8^2+3⋅8⋅5+5^2
209=13^2+3⋅13⋅1+1^2

Let f(n,r)
be the number of integers k not exceeding n that can be expressed as k=a^2+3ab+b^2, with a>b>0 integers, in exactly r different ways.

You are given that f(10^5,4)=237 and f(10^8,6)=59517

Find f(10^15,40)

Ich find die Grundidee von der Seite einfach nett. Es sind relative einfache Aufgabenstellungen, die man schön optimieren kann und sich daher auch als eine Art Wettbewerb eignen. Also hinsetzen, eine Lösung schreiben, optimieren und dann hier posten und vergleichen.

Noch ein paar Hinweise hier:
- f(n,r) soll nicht K ausrechnen, sondern zählen, wieviele K es gibt, die die Bedingungen erfüllen. (hatte ich erst falsch verstanden)
- Vielleicht geht es einigen so wie mir, quadratische Gleichungen hab ich mal vor Ewigkeiten nach einer Unbekannten aufgelöst. Das hier kann helfen: http://www.mathebibel.de/pq-formel
- Das Programm sollte unter 32Bit und 64Bit laufen - die Zahlen verlangen dadurch zwingend quad oder besser
- Es wird mit sehr großen Zahlen gearbeitet. Meine aktuelle Lösung braucht zum überprüfen von f(10^8,6)=59517 Stunden hier (32bit PB, E5400 @2,70Ghz - Arbeitspc, mal schaun, was meine Kiste zuhause braucht). 10^15 wird demnach deutlich länger brauchen. Nehmt also Zeit mit.
- Stellt den Debugger aus! er braucht einfach viel zu viel Zeit.
- ^ steht für die pow - funktion. a^2 kann man auch einfach durch a*a ausdrücken

Ich bin mal gespannt auf eure Lösungen. Meine Poste ich, sobald ich alles durchgerechnet habe

p.s.: Keine Ahnung wo man das hätte posten können....

[Helle]

Ja, sowas ist immer wieder interessant! Übrigens, keine Angst vor langen Rechenzeiten: Mein Test für f(10^8,6)=59517 braucht 875ms (i7-7700K@5.0GHz, reiner PB-Code). Mal sehen, wie f(10^15,40) hinzukriegen ist

[NicTheQuick]

Ich bin zu blöd die Aufgabe zu verstehen. Kann sie mir jemand nochmal erklären?

[Helle]

Da Du Code wohl besser verstehen wirst als meine blumigen Erklärungen hier mein Testcode für f(10^8,6)=59517. Für f(10^15,40) geht der so wegen Speichermangel eh nicht, da muss ich mir was anderes ausdenken.

Code: Alles auswählen

;a^2 + 3ab + b^2 = a*a + 3*a*b + b*b
;a>b>0
;a >= 2
;b >= 1
;f(10^8,6)=59517   10^8=Limit  6=Anzahl
Limit = 1e8
Anzahl = 6

Buffer = AllocateMemory(Limit)

a = 2
aa = 2
b = 1
TA = ElapsedMilliseconds()
For i = 1 To Sqr(Limit)
  X = 0
  While X <= Limit
    PokeB(Buffer + X, PeekB(Buffer + X) + 1)  ;schreibt fleißig in Null, wird aber unten nicht berücksichtigt
    X = a*a + 3*a*b + b*b
    a + 1
  Wend
  aa + 1
  a = aa
  b + 1
Next

j = 0
For i = 1 To Limit ;0 nicht, kann auch etwas höher sein (kleinster Startwert)
  If PeekB(Buffer + i) = Anzahl
    j + 1
  EndIf
Next
TE = ElapsedMilliseconds() - TA
MessageRequester("Anzahl", Str(j) + #LFCR$ + Str(TE))

Edit: For...Next braucht nur Sqr(Limit). Zeit damit 525ms

[GPI]

@nic

Die Zahl 209 kann durch die Formel a^2+3ab+b^2 auf zwei Arten ausgerechnet werden:
209=8^2+3⋅8⋅5+5^2
209=13^2+3⋅13⋅1+1^2

f(n,r) beschreibt die Anzahl der K die exakt r Lösungen für diese Formel haben: k=a^2+3ab+b^2, a>b>0, a und b sind Ganzzahlen.

Hier zwei Beispiele: f(10^5,4)=237, f(10^8,6)=59517

Finde f(10^15,40)

[NicTheQuick]

Ich habe es mittlerweile doch verstanden und anstatt mich zu melden lieber mal herum probiert. Für 10^15 wird es echt etwas knapp
Für 10^8 klappt meine Methode aber auch schon ganz gut, braucht aber dennoch wesentlich länger als deine, weil ich mich für die ollen Maps entschieden habe.

Immerhin 20 Sekunden auf einem i3 und 12 Sekunden auf meinem i7.

Code: Alles auswählen

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

n = 100000000
r = 6

; Zeit
time = ElapsedMilliseconds()

; Obergrenze für a
maxA = (Sqr(4 * n + 5) - 3) / 2
Debug maxA

; Häufigkeit von k
NewMap counts.i(10000000)

; Alle Kombinationen von a > b > 0 durchgehen

kA = 1     ; 1² + 3 * 1 * 0 + 0²
For a = 2 To maxA
	kA + a + a - 1
	a3 = 3 * a
	k = kA
	For b = 1 To a - 1
		k + b + b + a3 - 1
		If k > n	;maxA reicht nicht aus um sicherzugehen, dass k <= n bleibt
			Break
		EndIf
		;Debug "try: " + k + "  correct: " + Str(func(a, b)) + "  a:" + a + " b:" + b
		counts(Str(k)) + 1
	Next
Next

; Filtern nach Häufigkeit (r) und Zählen
result = 0
ForEach counts()
	If counts() = r
		result + 1
	EndIf
Next

time = ElapsedMilliseconds() - time

MessageRequester("Ergebnis", "n = " + n + #LF$ + "r = " + r + #LF$ + "f(n, r) = " + result + #LF$ + #LF$ + "Zeit: " + time + " ms")
;Debug result

Ich denke der Trick wird ein anderer sein. Vermutlich soll man alle k von 1 bis n in der äußeren Schleife durchlaufen und für jedes k einzeln prüfen, wie viele unterschiedliche Tupel (a, b) es gibt, damit die Gleichung k ergibt. Dazu könnte man sich aus der bisherigen Gleichung eine Wegfunktion herleiten, die der Höhenkurve folgt und dann nur alle ganzzahligen Tupel (a, b) zählen, die dabei auftreten. Ich weiß grad auch nicht wie ich es besser erklären soll.

Vielleicht gibt es aber auch irgendeine Regelmäßigkeit, die man ausnutzen kann. Müsste man vielleicht mal plotten.

[GPI]

Meine Lösung bis jetzt:

Code: Alles auswählen

EnableExplicit
OpenConsole()

Define a.d,b.d
Define h1.d
Define k.q,n.q,r.i,ll.i,myr.i
Define timer.i,timerx.i


N=Pow(10,8)
r=6

ll=0; hier kommt das ergebnis rein

k=10;niedrigster wert, wenn a=1 und b=2 (11) drunter macht k keinen sinn.

timer=ElapsedMilliseconds()
timerx=ElapsedMilliseconds()+5000

;For-Schleife geht nicht, da Quad
While k<n
  k+1
  myr.i=0 ;wie oft k dargestellt werden kann
  
  ;rechenmethode http://www.mathebibel.de/pq-formel
  a=IntQ((-1.5 + Sqr(1.25+k))   ) ; K ist gegeben, wenn B=1 ist, wie groß kann A dann maximal sein?
  
  ;Die schleife beginnt mit den größt möglichen a und zählt dann langsam runter.
  Repeat    
    h1.d=((3*a)/2.0); entspricht q ;)
    b=-h1 + Sqr(h1*h1-a*a+k) ; achtung b kann hier auch Null oder negativ sein!
    
    If b>0 And b<a And b=IntQ(b);k=a*a+3*a*b+b*b  Nur korrekt, wenn B keine Nachkommastellen hat. INTQ ist deutlich schneller als round()
      myr+1
      
      ; Wir haben breits mehr Lösungen als gesucht, von daher können wir aufhörenen
      If myr>r
        Break
      EndIf   
    ElseIf b>=a
      ; je kleiner a wird, um so größer wird b. b muss aber kleiner als a sein - ergo abbruch
      Break
    EndIf    
    a-1  
  Until a<=0.0
  
  If myr=r 
    ll+1
    If ElapsedMilliseconds()>timerx
      timerx=ElapsedMilliseconds()+5000
      PrintN( ""+k+" ("+StrF(k/n*100,2)+"%) "+StrF((ElapsedMilliseconds()-timer)/60000,2) +"m "+ll)
    EndIf
    
  EndIf
  
Wend

PrintN( "ergebnis:"+ll+"  "+StrF((ElapsedMilliseconds()-timer)/60000)+"m")
Input()

Ist verflucht langsam. Ich versuch gerade eure Methode zu verstehen....

[GPI]

Verdammt - man kann auch umgekehrt rechnen.... Meine Lösung braucht dafür fast keinen Speicher

Und bei der Fragestellung 10^15 und wenn man für jede häufigkeit ein Byte ansetzt (man kann ja aufhören zu addieren, wenn man schon über den Grenzwert ist, ein Byte ist also kein Problem):
1000000000000000 bytes = 976562500000 kb = 953674316,40625 mb = 931322,574615478515625 qb = 909,4947017729282379150390625 TB

Sprich: die schnelle Methode wird wohl nicht funktionieren. Selbst wenn man auf eine Datenbank ausweicht, man würde 909 TB nur für die Daten brauchen ohne Verwaltung.

Hier rechnet er gerade an der zweiten vorgegeben Lösung und komme auf: 71.07% aller möglichen Ks bereits durchgegangen, 129.79 Minuten dafür gebraucht und bisher 41142, die alle Bedingungen erfüllen.
Da der Aufwand steigt, je höher die Zahlen werden, rechne ich mit aktuell nochmal 130minuten bis er das erste Ergebnis ausgerechnet hat

Die 10^15-Aufgabe erledige ich aber lieber zuhause mit meinen I7, der dürfte deutlich schneller sein. Hat jemand schon versucht, es mit Threads zu versuchen, also zwei Threads laufen zu lassen, die dann jeweils k um 2 erhöhen und somit jeder die Hälfte rechnet - müsste dann sich nicht besser auf die CPU-Kerne verteilen?

edit: FERTIG! also das kleine:
ergebnis:59517 217.4097442627 minuten

edit2:
mir ist eventuell eine Lösung eingefallen, die meine beschleunigen könnte - kann man feststellen, ob eine wurzel glatt aufgeht? (edit: eher unglatt, obs mit 0.5 endet oder nicht) damit könnte man die SQR-Funktion deutlich beschleunigen.

Ansonsten ist das einzige was mir grob als Idee einfällt: B interessiert mich eigentlich nicht, ich möchte nur wissen, ob ein B überhaupt exestiert....
hmmm.... mir kommt gerade was. die formel nach a auflösen und in eine formel, die nach b aufgelöst ist einsetzen. dann müsste man doch a ausreichnen können.... wenn die glatt ist, noch überprüfen ob b da ist. Das dürfte deutlich schneller gehen als meine lösung.... grob mal testen.
und muhahahahahaha - wurzeln viele wurzeln. Ne geht nicht. zumal noch ein denkfehler, a muss ja mehrere Lösungen geben, bringt also nichts....

Das zu optimieren ist alles andere als leicht. Eine Gesetzmäßigkeit wird schwer werden.

[GPI]

so ich hab jetzt meine Routinen etwas optimiert und eine Thread-Lösung eingebaut. Windows ist tatsächlich so intelligent die Last auf alle Virtuellen Kerne zu verteilen. Und wenn man alle nutzt, dann ruckelt sogar der Mauszeiger

Code: Alles auswählen

EnableExplicit
OpenConsole()
PrintN(#PB_Compiler_File)


Structure daten
  n.q  
  k.q
  dk.q
  r.i
  ll.i
  done.i
EndStructure


Procedure MyThread(*info.Daten)
  Protected a.d,b.d,bb.i
  Protected k.q,n.q,r.i,ll.i,myr.i
  Protected dk.q
  k=*info\k
  dk=*info\dk
  n=*info\n
  r=*info\r
  
  Repeat ;k+dk<=n
    k+dk
    If k>n
      Break
    EndIf
    
    myr.i=0 
    
    a=IntQ((-1.5 + Sqr(1.25+k))   ) 
    
    Repeat    
      
      b=-1.5*a + Sqr(1.25*a*a+k)    
      If b<a
        bb=Int(b)
        If  b=bb
          myr+1
          If myr>r
            Break
          EndIf  
        EndIf
        
      Else
        Break
      EndIf
      a-1
      
    ForEver
    
    If myr=r
      *info\ll +1
      *info\k=k
    EndIf
    
  ForEver
  
  *info\done=#True
EndProcedure

Define pro,do,i
Define timer,myk.d
Define k.q,n.q,r.i,ll.i
Define text.s

N=Pow(10,8)
r=6
pro=7


ll=0; hier kommt das ergebnis rein
k=10;niedrigster wert, wenn a=1 und b=2 (11)

Dim info.daten(pro)
For i=1 To pro
  info(i)\k=k+i-1-pro
  info(i)\dk=pro
  info(i)\n=n
  info(i)\r=r
  CreateThread(@MyThread(),@info(i))
Next

timer=ElapsedMilliseconds()
Repeat
  Delay(5000)
  do=#True
  ll=0
  myk=0
  For i=1 To pro
    If info(i)\done=#False
      do=#False      
    EndIf
    ll+info(i)\ll
    myk+(info(i)\k/pro)
  Next

  PrintN( ""+IntQ(myk)+" ("+StrF(myk*100/n,2)+"%) "+StrF((ElapsedMilliseconds()-timer)/60000,2) +"m "+ll)  
Until do


PrintN( "ergebnis:"+ll+"  "+StrF((ElapsedMilliseconds()-timer)/60000)+"m")
Input()

die 10^8-Aufgabe hab ich jetzt erledigt:
Ergebnis:59517 12.1759662628 Minuten

Was so ein i7-4770k mit 4+4Kernen ausmacht
Die Variable pro regelt, wieviele Threads erstellt werden sollen, theoretisch gingen bei mir auch 8, ich wollte aber noch etwas fürs OS etc. übrig lassen, damit man noch nebenbei etwas machen kann.

[Imhotheb]

Ich musste jetzt auch mal GPI's Code testen.

trotz meines etwas angestaubten i7-3930k:
ergebnis:59517 8.8388004303m

allerdings mit pro=11 (bei 12 wird es wieder langsamer)
Aber durch die Threads ist es doch spürbar schneller.