Problem z zadaniem: długie wyszukiwanie 'ciekawych liczb'

Witam, próbuję rozwiązać takie zadanie: http://www.main.edu.pl/user.phtml?op=showtask&task=cie&con=ALG
Zgodnie z tym co było w kursie powinienem wykorzystać technikę preprocessingu, czyli znaleźć wszystkie 'ciekawe liczby' i umieścić w gotowych programie (zapewne nie jest ich dużo). Mój problem pojawia się w ich znalezieniu, ponieważ zakres sięga aż do 100 000 000. Program do wyszukiwania, który napisałem liczyłby to jakieś 28h na moim PC. Jak mogę go zoptymalizować?

 
#include <cstdio>
#include <fstream>

using std::ofstream;

int dziel, // ile dzielników
    suma;  // suma liczb pierwszych

void pierwsze(int x, int temp = 2)
{
	for(int i = 2; i <= x/2; i++)
	{
		if(temp > 2) break;
		if(x % i == 0) temp++;
	}
	if(temp == 2) 
		suma += x;
}

void dzielniki(int x)
{
	dziel = 2;
	for(int i = 2; i <= x/2; i++)
		if(x % i == 0) 
		{
			dziel++;
			pierwsze(i); // sprawdzenie czy nie jest liczbą pierwszą
		}
}

int main()
{
	ofstream log;
	log.open("log.txt");
	int x;	// liczba od której zacznie się poszukiwanie
	scanf("%d", &x);
	for(int i = x; i < 100000001; i++)
	{
		dzielniki(i);
		printf("%d\n", i);
		if(dziel % 13 == 0 && dziel/13 == suma) // warunek 'ciekawej liczby'
		{
			printf("%d\n", i);
			log << i << "\n"; // zapisanie znalezionej liczby do pliku
		}
		suma = 0;
	}
	log.close();
	return 0;
}

Po pierwsze, najpierw stwórz sobie tablicę liczb pierwszych mniejszych niż milion. Możesz to zrobić osobnym programem np wykorzystując sito Hornera, zaś do tego programu wstawić już gotową tablicę liczb pierwszych - to się nazywa preprocesing.
Po drugie, sprawdzaj podzielność liczby tylko przez liczby pierwsze, i zliczaj ilość, np. 126000 to 2⁴³25³⁷1 więc masz zbiór (2,2,2,2,3,3,5,5,5,7) - ilość dzielników liczby 126000 to ilość różnych podzbiorów tego zbioru - da się obliczyć wg równania.
Po trzecie zamiast sprawdzać czy liczb jest liczbą ciekawą w nieciekawy sposób:
dziel % 13 == 0 && dziel/13 == suma
wystarczy że sprawdzisz:
suma13==dziel
jedno mnożenie jest o wiele szybsze nawet od jednego dzielenia.

Nie chodziło Ci o Sito Eratostenesa? A liczby pierwsze powyżej miliona liczyć w tym programie? Nie lepiej zrobić już tablicę ze wszystkimi potrzebnymi liczbami pierwszymi (do 50 000 000)?
O równaniu na dzielniki liczby całkiem zapomniałem, zaraz go dodam.
Z tym dzieleniem też głupi błąd, nie pomyślałem, że można go tak prosto zapisać. Dzięki za pomoc

Chciałem zrobić taki program, napisałem:

 #include <iostream>
#include <cmath>

using namespace std;

bool spr(int a)
{
   for(int i=2; i<=sqrt(a); i++)
     if(a%i == 0)
        return false;
        
   return true;
}

int dz(int a)
{
    int ile = 1;
    
    for(int i=1; i <= a/2;  i++)
      if(a%i == 0)
         ile++;
         
    return ile;     
}

int suma (int a)
{
    int suma = a;
    
    for(int i=1; i <= a/2; i++)
     if(a%i == 0)
        if(spr(i))
           suma +=i;
    
    return suma;
}

int main()
{
    
    for(long i=331775; i<=331777; i++)
    {
       
       if(suma(i)*13==dz(i))
         {
           cout  << i << endl;
           break;
         }      
    }
       
    system("PAUSE");
    return EXIT_SUCCESS;
}

Powinno wyświetlić się: 331776, niestety tak się nie dzieje. Nie wiem w czym jest problem :(Proszę o pomoc

Wprowadziłem te wszystkie zmiany, ale program nadal chodzi powolnie. Nie jestem za dobry w optymalizacji :/ Może ktoś ma pomysł jak go zoptymalizować albo jak inaczej znaleźć te wszystkie 'szczęśliwe liczby'?

#include <cstdio>
#include <fstream>
#include <vector>
#include <map>

using std::ofstream;
using std::vector;
using std::map;

int dziel = 1, // ile dzielników
    suma;     // suma liczb pierwszych
const int n = 50000000;				
bool* numbersTable = new bool[n];  // tablica na sito Eratostenesa
vector<int> lp;  // tablica liczb pierwszych
map<const int, int> divisor;  // liczby przy rozkładzie na czynniki pierwsze, first->podstawa, second->wykładnik

void sito()
{
	for (int i = 2; i*i <= n; i++ )  // przeszukuj liczby od 2 do sqrt(n), 0 i 1 nie są liczbami pierwszymi
	{
		if (numbersTable[i] == true)   // jeżeli dana liczb jest już wykreślona
			continue;                   // to przejdź do kolejnej
		for (int j = 2 * i ; j <= n; j += i)   // przejdź od liczby 2 * i do n przesuwając się o i
			numbersTable[j] = true;	   // i każdą z nich usuwaj ze zbioru
	}
	for(int i = 2; i < n; i++)   // wpisywanie do tablicy liczb pierwszych
		if(numbersTable[i] != 1)
			lp.push_back(i);
}

int dzielniki(int x)
{
	if(numbersTable[x] != 1) return 1;   // jeżeli liczba pierwsza to wyjdź
	int l = 0, 
		y = lp[l];
	do   // rozkład liczby na czynniki pierwsze
	{
		if(x % y == 0) 
		{
			x /= y;
			divisor[y]++;
		}
		else { l++; y = lp[l]; }
	} while(x/y != 1);
	divisor[y]++;
	suma = 0;
	dziel = 1;
	for(map<const int, int>::const_iterator it = divisor.begin(); it != divisor.end(); suma += it->first, it++)
		dziel *= (it->second + 1);	// mnożenie kolejnych potęg po rozkładzie liczby
	return 0;
}

int main()
{
	sito();
	ofstream log;
	log.open("log.txt");
	int x = 3;	// liczba od której zacznie się poszukiwanie
	for(int i = x; i < 2*n; i++)
	{
		dzielniki(i);
		if(dziel == 13*suma)	// warunek 'ciekawej liczby'
		{
			printf("%d\n", i);		// wyświetlenie i 
			log << i << "\n";			// zapisanie znalezionej liczby do pliku
		}
	}	
	log.close();
	return 0;
}
 

napisałem sobie program do pisania takich programów, wyjszło mu takie coś:

int sprp2(unsigned long long n) {if(8503056<n)if(42515280<n)if(64000000
<n)return 98537472==n ? -1 : 98537472;else return 64000000==n ? 98537472
 : 64000000;else if(33554432<n)return 42515280==n ? 64000000 : 42515280
;else if(12644352<n)return 33554432==n ? 42515280 : 33554432;else
return 12644352==n ? 33554432 : 12644352;else if(5419008<n)if(7680000
<n)return 8503056==n ? 12644352 : 8503056;else if(6967296<n)return 7680000
==n ? 8503056 : 7680000;else return 6967296==n ? 7680000 : 6967296;
else if(1658880<n)return 5419008==n ? 6967296 : 5419008;else if(331776
<n)return 1658880==n ? 5419008 : 1658880;else return 331776==n ? 1658880
 : 331776;}

trochę przeformatowałem

int odp(int n) {
	if(8503056<n)
		if(42515280<n)
			if(64000000<n)return 98537472==n ? -1 : 98537472;
			else return 64000000==n ? 98537472 : 64000000;
		else 
			if(33554432<n)return 42515280==n ? 64000000 : 42515280;
			else 
				if(12644352<n)return 33554432==n ? 42515280 : 33554432;
				else return 12644352==n ? 33554432 : 12644352;
	else 
		if(5419008<n)
			if(7680000<n)return 8503056==n ? 12644352 : 8503056;
			else 
				if(6967296<n)return 7680000==n ? 8503056 : 7680000;
				else return 6967296==n ? 7680000 : 6967296;
		else 
			if(1658880<n)return 5419008==n ? 6967296 : 5419008;
			else 
				if(331776<n)return 1658880==n ? 5419008 : 1658880;
				else return 331776==n ? 1658880 : 331776;
}

jak dodać

main(){
	int i;
	scanf("%d", &i);
	printf("%d", odp(i));
}

to będzie OK
czyli programy są dwa, pomocniczy znajdzie "ciekawe" liczby, a ten właściwy z nich skorzysta

"ciekawe" liczby znalazłem tak:
PARI/GP> for(i=1,100000000,v=factor(i);l=#v~;if(13*sum(i=1,l,v[i,1])==prod(i=1,l,v[i,2]+1),print(i)))
sum(...v[i,1]) to suma dzielników pierwszych
prod(...v[i,2]+1) to liczba dzielników
gdyby kto chciał to jutro zapiszę to w C (faktoryzację da się chyba dość wcześnie przerwać, jażeli tak to może wystarczy jeden program?)

mój programik potrafi jeszcze tak:

#define         _1      int sprp2(unsigned long long n) {
#define         _2      if(
#define         _3      8503056
#define         _4      <n)
#define         _5      42515280
#define         _6      64000000
#define         _7      return
#define         _8      98537472
#define         _9      ==n ?
#define         _10     -1
#define         _11     :
#define         _12     ;
#define         _13     else
#define         _14     33554432
#define         _15     12644352
#define         _16     5419008
#define         _17     7680000
#define         _18     6967296
#define         _19     1658880
#define         _20     331776
#define         _21     }
#define         _22     //\
#define         _23     _1 _2 _3 _4 _2 _5 _4 _2 _6 _4 _7 _8 _9 _10 _11 _8 _12 _13 _7 _6 \
_9 _8 _11 _6 _12 _13 _2 _14 _4 _7 _5 _9 _6 _11 _5 _12 _13 _2 _15 _4 _7 _14 _9 _5 _11 \
_14 _12 _13 _7 _15 _9 _14 _11 _15 _12 _13 _2 _16 _4 _2 _17 _4 _7 _3 _9 _15 _11 _3 \
_12 _13 _2 _18 _4 _7 _17 _9 _3 _11 _17 _12 _13 _7 _18 _9 _17 _11 _18 _12 _13 _2 _19\
_4 _7 _16 _9 _18 _11 _16 _12 _13 _2 _20 _4 _7 _19 _9 _16 _11 _19 _12 _13 _7 _20 _9\
_19 _11 _20 _12 _21 _22 
_23 3.141592653589793238462643383

Nie wiem jak Ci dziękować za te liczby :) Musiałem tylko troche przeedytować ten programik, żeby spełniał wymagania zadania i wszystko śmiga.
Nie rozumiem tylko jak znalazłeś te liczby. To co napisałeś jest dla mnie troche za pogmatwane.
A to z użyciem define to już jest masakra ;)

napisałem sobie na ideone.com programik, w 10 sekundy dotarłem do 50 mln. (połowa roboty bo czas nie jest kwadratowy:-)
(ideone daje tylko 15 sek. więc nie dotrę do końca)
zacząłem od najprostszego programu rozkładającego na czynniki
przy każdej okazji znalezienia czynnika uzupełniam sumę (*13) i iloczyn wykładników (+1) potrzebne do policzenia liczby dzielników
jak liczba już całkiem mi się wydzieli :-) to sprawdzam SUMA==ILOCZYN i jeżeli tak to mam "ciekawskiego"
warto popatrzeć jak w czasie rozkładu rosną suma i iloczyn

Metoda generowania liczb Xitami (Xitami-ego?) jest OK i pozwoli sobie na wygenerowanie tych liczb (w domu). Ten odczyt liczb to chyba jakiś dziwny żart :D
Najprościej - wrzucasz wszystkie poznane przez siebie liczby ciekawe do tablicy i wyszukujesz w pętli (liniowo; można też wyszukiwaniem binarnym, ale po co się męczyć dla 10 elementów?) odpowiednią liczbę...

Z tablicą byłoby by prościej, ale mniejsza o to. @Xitami jak udało Ci się uzyskać taką szybkość? Zmieniłem kod wg tego co napisałeś, ale nadal się wlecze :/
Możesz wstawić swój program?

 
#include <cstdio>
#include <vector>
#include <iostream>

using std::cout;
using std::vector;

const int n = 50000000;				
bool* numbersTable = new bool[n];	// tablica na sito Eratostenesa
vector<int> lp;		          // tablica liczb pierwszych

void sito()
{
	for (int i = 2; i*i <= n; i++ )			// przeszukuj liczby od 2 do sqrt(n), 0 i 1 nie są liczbami pierwszymi
	{
		if (numbersTable[i] == true)		// jeżeli dana liczb jest już wykreślona
			continue;						// to przejdź do kolejnej
		for (int j = 2 * i ; j <= n; j += i)	// przejdź od liczby 2 * i do n przesuwając się o i
			numbersTable[j] = true;				// i każdą z nich usuwaj ze zbioru
	}
	for(int i = 2; i < n; i++)
		if(numbersTable[i] != 1)
			lp.push_back(i);
}

bool happy(int x)
{
	if(numbersTable[x] != 1) return false;	// jeżeli liczba pierwsza to wyjdź
	int dziel = 1,	// ile dzielników
		suma = 0,	// suma liczb pierwszych
		l = -1, 
		y = 0;
	while(x != 1)
	{
		l++;
		y = lp[l];
		if(x % y == 0)
			suma += y;
		else continue;
		
		int div = 0;
		while(x % y == 0)
		{
			x /= y;
			div++;
		}
		dziel *= div + 1;
	}
	suma *= 13;
	if(dziel == suma) return true;
	else return false;
}

int main()
{
	sito();
	int x = 5419000;
	for(int i = x; i < 2*n; i++)
	{
		if(happy(i))
			printf("%d\n", i);

	}
	system("pause");
	return 0;
}

To może ułatwienie: jeśli liczba n = p^{x * q}y * r^{z * ..., gdzie p,q,r,.. to liczby pierwsze, to liczba dzielników tej liczby to (x+1)(y+1)(z+1)...} . Na przykład: 5625 = 2</sup>1 * 3^{2 * 5}3, więc liczba dzielników tej liczby to (1+1)(2+1)(3+1)=24.
Inne zauważenie jest takie, że każda liczba pierwsza dzieląca n jest nie większa niż pierwiastek kwadratowy z n (nietrudno to udowodnić). Przykład ulepszonego algorytmu w pseudokodzie:

liczbadzielnikow := 1;
pierwiastekzn := sqrt(n);
dla każdej liczby pierwszej p <= pierwiastekzn:
    ilerazydzieli := 0;
    dopóki n jest podzielne przez p:
        n := n / p;
        ilerazydzieli++;
    liczbadzielnikow *= (ilerazydzieli+1);

zastosowałem pewien chwyt
nie czekam na całkowite "wydzielenie", przerywam wcześniej (z negatywną odpowiedzią)
popatrzcie na rozkład na czynniki liczb szczęśliwych

o mam jeszcze jednego pomysła (a nawet dwa), muszę to sprawdzić
-------------------------------------------------------------------[ dodane
jeden pomysł to zastosowanie sita, a przy okazji zbieranie w nim sumy i liczby dzielników
ale to do d.... bo czas kwadratowy
wcześniej nie miałem tablic

drugi pomysł to podejście do sprawy od tyłu
manipulując tablicą z wykładnikami da się otrzymać różne liczby m.in. szczęśliwe
jeszcze nad tym pomyślę

ale najważniejsze, liczba dzielników liczb <100mln nie jest większa niż 4032
Highly composite numbers

@mnbvcX Ułatwienie z (x+1)(y+1)(z+1)... już miałem zastosowane tylko może trochę źle zaimplementowane(już poprawiłem). Po dodaniu drugiego ułatwienia program znacznie przyspieszył, ale to nadal ok. 20k/1s. Tylko nie rozumiem dlaczego "każda liczba pierwsza dzieląca n jest nie większa niż pierwiastek kwadratowy z n". Np. biorąc liczbę 18922, jej dzielnikiem jest liczba pierwsza 9461. Może coś źle rozumiem, dlatego pytam.

@Xitami Zaraz sprawdze ten chwyt, ale jak to policzyć szybko bez niego? Bo na początku nie znamy tych liczb, więc nie można zauważyć żadnych cech wspólnych. A info z max liczbą dzielników może się przydać :) Z tą manipulacją też ciekawy pomysł.

dzielniki zawsze występują w parach, w takiej parze jeden jest mniejszy(=) od pierwiastka, a drugi większy np. (1, n), (a, n/a), (b, n/b), nie warto sprawdzać podzielności przez liczbę większą od pierwiastka bo wcześniej sprawdziliśmy
już jej mniejszego partnera, czyli,
gdy nie znalazłem podzielnika <= pierwiastkowi z liczby to nie warto sprawdzać dalej bo nic tam już nie znajdę.

jaka potrzebna nam jest największa liczba pierwsza?

OK, już rozumiem. Największa to pierwiastek z 1mld (największej możliwej szczęśliwej liczby).

po sprawdzeniu tylko 3569 liczb, otrzymałem 11 liczb szczęśliwych w czasie grubo poniżej sekundy!!!
liczba 1959874560 jest szczęśliwa, podobnie jak 19163331719168
mam ich już ponad tysiąc
tylko muszę przemyśleć jeszcze jak ten program zatrzymać :)

I rzeczywiście sprawa polega n rozkładzie na czynniki, ale liczby są bardzo małe.
W sumie to dałoby się to zadanie rozwiązać z pomocą kalkulatora w czasie picia kawy.
Każdą liczbę można jednoznacznie przedstawić tak: n=a^x * b^y * c^z * ...
Gdzie a, b, c, ... to kolejne liczby pierwsze, natomiast x,y,z to liczby całkowite (z zerem)
liczbę dzielników liczby obliczamy ze wzoru (x+1)(y+1)(z+1)...

Weżmy liczby postaci a^x*b^y, suma dzielników pierwszych tej liczby to a+b
liczba ta ma (x+1)(y+1) dzielników naturalnych
szukamy więc rozwiązania róównania 13(a+b) = (x+1)(y+1)
bardzo łatwo znaleźć najmniejszą szczęśliwą, weźmy liczby 2 i 3
13(2+3)=(x+1)(y+1), akurat 2+3 jest liczbą pierwszą ale nawet gdy tak nie jest to żaden problem, 135 = (x+1)(y+1)
naszymi wykładnikami będzie para 12 i 4,
i tak otrzymujemy dwie liczby szczęśliwe: 2^12*3^4 oraz 2^4*3^12
z pary liczb pierwszych można uzyskać różną liczbę szczęściarzy, zależy to od tego ile dzielników ma liczba 13(a+b)
jeżeli zaczniemy od 2 i 11 otrzymamy tylko jedną liczbę, bo 2+12=13, a 1313=1313
zato dla 13, 17 uzyskamy 14 liczb, a to dla tego, że 13+17=30, a 30 to 235, do tego dochodzi jeszcze 13 czyli mamy 2222=16 dzielników, dwie możliwości oczywiście odpadają, wykładnik danej liczby nie może być zerem bo wtedy suma dzielników nie będzie równa a+b.
Ale... co z samotną trzynastką? 13(13)=(x+1), 169-1=x, no to mamy wyjątkowo szczęśliwą liczbę 13^168, liczba ta ma 188 cyfr, oto ona
13882995648971427331025903931877308650968871447
34723560521135815725214852166231801067608072301
97591538215048913800351981208465401843465386281
61904387790998345046721749399293824585785449121
każda liczba postaci p^(13p-1) jest szczęśliwa gdy p jest pierwsze
np.
ile cyfr ma liczba mówiąca o tym ile cyfr ma liczba określająca liczbę cyfr tej liczy? :-)

więc nie szukajmy tych liczb tylko je konstruujmy.

najmniesza liczba w rozkładzie której pojawia się 7 jest poza zakresem zadania, czyli interesują nas tylo liczby o postaci
2^x*3^y*5^z gdzie x, y, z to liczby naturalne (z zerem!)
prawda ale tylko w przypadku iloczynu 2 liczb pierwszy

Genialne! :) Szkoda, że sam na to nie wpadłem. Tylko trochę ciężko to zapisać w c++ żeby uwzględnić wszystkie możliwe rozkłady. Jeśli dobrze rozumiem, trzeba podstawiać liczby do (x+1)(y+1)... jak się strony zrównają to mamy szczęściarza.. a i liczba 98 537 472 ma w rozkładzie 11, a też sie mieści w zakresie.

dla dwu pierwszych szukamy rozwiązania 13(a+b)=(x+1)(y+1)
można zgadywać x i y np.
for x=2..??
_ for y=x+1..??
_ _ if x*y == ab13 then happy2(x-1, y-1)
ale jak daleko szukać?

a można lepiej poszukać dzielników liczby 13(a+b) ?
for(k=2; k*k<=ab13 ; k++)
_ if ab13 % k == 0 happy2(k-1, ab13/k-1)

dla trzech pierwszych podobnie, w pętli szukamy dzielnika 13(a+b)
on (-1) będzie pierwszym wykładnikiem, oznaczmy go d
teraz szukamy dzielników liczby 13(a+b)/d
no i mamy trzy wykładniki, mamy albo nie bo po drodze można trafić na liczby pierwsze i szczęściarza nie znajdziemy

Niby rozumiem o co chodzi, ale przy trzech pierwszych się gubię. Poza tym trzeba uwzględnić wszystkie możliwości np. 2, 5 albo 2, 11...

Coś mi cały czas umykało, wróciłem do programu sprawdzającego po kolei zakres 1..100mln
czas na ideone.com 0,85 sekundy

już bez sztuczek, mam dowód poprawności :-)
łatwo odsiewam większość
pozostałe 24602 rozkładam na czynniki pierwsze
(wcześniej przygotowałem sobie tablicę liczb pierwszych od 2 do 157

Możesz wkleić kod programu, żebym mógł przeanalizować?

Funkcja testująca wygląda jakoś tak

int test(int n){
	.....
	.....tu rozkład i inne takie
	.....
	return ( 13*s == w );     // no właśnie !!! oczywista oczywistość
}

! 13s==w

Czytasz dalej, czy już wiesz? :-)

to samo inaczej s==w/13

Żeby tak mogło być jeden z czynników musi występować w potędze
która dzieli się przez 12
czyli liczba musi dzielić się przez 2^{12==4096 lub 3}12==531441,
a przez 5^12 już nie musi bo to więcej niż 100 mln

Trzeba pamiętać o potęgach liczb pierwszych o postaci p^(13p-1)
czyli 2^{(132-1)==33554432, 3}(133-1) jest większe od 100 mln

no to mamy trzy liczby przez które MUSI dzielić się liczba z podanego zakresu
dodajmy jeden wiersz do funkcji testującej if(n%4096!=0 && n%531441!=0 && n!=33554432) return 0;To sprawdzenie przechodzi już tylko 24'602 kandydatów.

Zakres liczb to 100e6, ale skoro liczba musi dzielić się prze którąś z tych trzech
to iloraz nie może być większy od 100e6/4096=24414
Iloraz jest albo liczbą pierwszą, albo będzie miał dzielniki pierwsze <= sqrt(24414)=156
czyli sprawdzamy podzielność testowanej przez mniej niż 40 liczb pierwszych

include <stdio.h>
int l=0;
int test(int n){
	static const int sp[37]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,
	                73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157};
	if(n%4096!=0 && n%531441!=0 && n!=33554432) return 0;
	l++;
	int s=0, q=1, w, i;
	for(i=0; i<37; i++) 
		if(n%sp[i]==0) {
			s+=sp[i];
			w=2;
			n/=sp[i];
			while(n%sp[i]==0){
				w++; n/=sp[i]; }
			q*=w;}
	if(n>1) {
		s+=n;
		q*=2;}
	return 13*s == q; }

int main(void) {
	int i;
	for(i=1;i<100000001;i++)
		if(test(i))
			printf("%d\n",i);	 
	printf("\n\n%d\n",l);
	return 0;}

Xitami napisał(a)

Żeby tak mogło być jeden z czynników musi występować w potędze
która dzieli się przez 12

A nie powinno być w 'potędze + 1', które dzieli się przez 13? np. (25+1) % 13 = 0
Wtedy zamiast n != 33554432 można by dodać do warunki n % 2^25 != 0, dalej już jest poza zakresem.
A i wielkie dzięki za pomoc przy tym zadaniu :)