[Vari] Contest 9: Divisori dei numeri triangolari

[VICIUS]

Ieri sera mi sono imbattuto in questo problema su un altro sito e l'ho trovato veramente carino. Non è un contest ai livelli di quelli di gugoXX quindi dovreste riuscire ad arrivare alla soluzione finale in poco tempo. Usatelo per staccare un po' e divertirvi

Definizioni
L' n-esimo numero triangolare si ottiene sommando tutti i numeri naturali da 1 fino ad n. Quindi il primo numero triangolare è 1, il secondo è 3 e il terzo è 6.

Per trovare l'ennesimo numero triangolare senza dover fare la sommatoria ogni volta si può usare la formula di Gauss.

T(n) = (n * (n + 1)) / 2

Per divisore in questo caso si intende un numero intero positivo d che divide il numero triangolare senza produrre resto.

Problema
Il problema consiste nello scrivere un programma trovi il primo numero triangolare che ha più di x divisori.

 n  |  Tn | Divisori 
--------------------------
 1  |  1  | 1
 2  |  3  | 1,3
 3  |  6  | 1,2,3,6
 4  |  10 | 1,2,5,10
 5  |  15 | 1,3,5,15
 6  |  21 | 1,3,7,21
 7  |  28 | 1,2,4,7,14,28

Osservando la tabella qui sopra e sapendo che x è 4 allora la risposta giusta è 7 perché il settimo numero è il primo con più di 4 divisori.

Brute force

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

int count_divisors(int n)
{
	int i, count;
	
	for (i = 1, count = 0; i <= n; ++i)
		if ((n % i) == 0) ++count;
	
	return count;
}

int find_nth_triangular_number(int i)
{
	return (i * (i + 1)) / 2;
}

int main(int argc, char** argv)
{
	int d, n = 1, i = 1;
	
	if (argc != 2) return -1;
	if ((d = atoi(argv[1])) < 1) return -1;
	
	printf("Inizio la ricerca del primo numero triangolare che ha più di %d divisori\n", d);
	while (count_divisors(n) < d) {
		n = find_nth_triangular_number(++i);
	}
	printf("Si tratta del numero triangolare %d che vale %d\n", i, n);
	
	return 0;
}

La forza bruta è semplice da implementare e funziona anche piuttosto velocemente per valori piccoli di x ma se si comincia a salire sopra 100 allora la storia cambia. Sul mio computer impiega circa 7 secondi per trovare il numero con più di 200 divisori.

Io ho risolto il problema in due modi completamente diversi. Quello più veloce, per fare un paragone con la versione brute force, impiega circa 16 millesimi di secondi per trovare il numero con più di 1000 divisori.

Vi allego un po' di output del programma in C così potete fare un po' di confronti per capire se la vostra soluzione è corretta.

Inizio la ricerca del primo numero triangolare che ha più di 10 divisori
Si tratta del numero triangolare 15 che vale 120

Inizio la ricerca del primo numero triangolare che ha più di 50 divisori
Si tratta del numero triangolare 224 che vale 25200

Inizio la ricerca del primo numero triangolare che ha più di 100 divisori
Si tratta del numero triangolare 384 che vale 73920

Inizio la ricerca del primo numero triangolare che ha più di 200 divisori
Si tratta del numero triangolare 2015 che vale 2031120

Divertitevi

[m.distrutti]

perdona la domanda VICIUS ma la soluzione l'hai fatta te ihih, il nostro compito qual'e'? fornirne altre più efficienti ?

[VICIUS]

Originariamente inviato da m.distrutti

perdona la domanda VICIUS ma la soluzione l'hai fatta te ihih, il nostro compito qual'e'? fornirne altre più efficienti ?

Come per tutti gli altri contest l'obiettivo ultimo è divertirsi. E poi mica è detto che la mia sia la soluzione migliore.

[m.distrutti]

Originariamente inviato da VICIUS

Come per tutti gli altri contest l'obiettivo ultimo è divertirsi. E poi mica è detto che la mia sia la soluzione migliore.

io veramente ad una prima analisi avrei postato la tua ahah, è il primo contest che vedo nel forum perciò scusa la domanda un pò idiota

[VICIUS]

Originariamente inviato da m.distrutti

io veramente ad una prima analisi avrei postato la tua ahah, è il primo contest che vedo nel forum perciò scusa la domanda un pò idiota

Intendi quello scritto in C? Se si quello è solo per confronto, non è di certo la soluzione migliore.

[shinya]

Ok, fuori concorso... la mia versione brute force (molto "brute" purtroppo...) in factor:

USING: kernel math.ranges sequences locals ;

:: divisors ( n -- seq ) n dup [1,b] [ n swap mod zero? ] filter nip ;
: nth-triangular ( x -- n ) dup 1 + * 2 / ;
: count-divisors ( n -- m ) nth-triangular divisors length ;
:: find-number ( n x -- m ) n count-divisors x <= [ n 1+ x find-number ] [ n ] if ;
: contest-9 ( x -- ) 2 swap find-number . ;

Output di esempio:

( scratchpad ) 4 contest-9
7
( scratchpad ) 100 contest-9
384

[DanieleC88]

Originariamente inviato da VICIUS

 n  |  Tn | Divisori 
--------------------------
 1  |  1  | 1
 2  |  3  | 1,2
 3  |  6  | 1,2,3,6
 4  |  10 | 1,2,5,10
 5  |  15 | 1,3,5,15
 6  |  21 | 1,3,7,21
 7  |  28 | 1,2,4,7,14,28

Intendevi 1, 3?

[DanieleC88]

Originariamente inviato da VICIUS

Io ho risolto il problema in due modi completamente diversi. Quello più veloce, per fare un paragone con la versione brute force, impiega circa 16 millesimi di secondi per trovare il numero con più di 1000 divisori.

Tanto per curiosità, qual è il numero triangolare che viene trovato? Voglio verificare che la mia soluzione sia corretta.

[VICIUS]

Originariamente inviato da DanieleC88

Intendevi 1, 3?

Si certo. Ora Correggo.

Originariamente inviato da DanieleC88

Tanto per curiosità, qual è il numero triangolare che viene trovato? Voglio verificare che la mia soluzione sia corretta.

Si tratta del numero triangolare 41040 che vale 842161320

[DanieleC88]

Propongo questa prima soluzione, è poco più di una brute force. Il tempo impiegato sulla mia macchina per la ricerca del numero triangolare con più di 1000 divisori è spropositato rispetto a quello di VICIUS (ben 26 secondi), quindi proverò a spremerlo il possibile ma temo che dovrò escogitare un metodo più "furbo".

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

#define gauss(n) (unsigned)(((n) * ((n) + 1)) / 2)

unsigned divisors(unsigned n)
{
	unsigned count = 2;
	unsigned i, last;

	if (n <= 1)
	{
		return n;
	}

	last = (unsigned) sqrt(n);
	for (i = 2; i <= last; ++i)
	{
		if ((n % i) == 0)
		{
			count += 2;
		}
	}
	return count;
}

unsigned find(unsigned count, unsigned *index, unsigned *number)
{
	unsigned n = count;
	unsigned sum, d;

	sum = gauss(n);
	d = divisors(sum);

	while (d <= count)
	{
		sum += (++n);
		d = divisors(sum);
	}

	(*index) = n;
	(*number) = sum;
	return d;
}

int main(int argc, char *argv[])
{
	time_t t1, t2;
	unsigned count;
	unsigned index;
	unsigned number;
	unsigned divisors;

	--argc;
	++argv;

	if (argc != 1)
	{
		fprintf(stderr, "Numero errato di argomenti.\n");
		return EXIT_FAILURE;
	}

	count = atoi(argv[0]);

	t1 = clock();
	divisors = find(count, &index, &number);
	t2 = clock();

	printf("Ricerca completata in %lf secondi.\n", (double)(t2 - t1) / (double) CLOCKS_PER_SEC);
	printf("Il numero triangolare cercato sembra essere T%u = %u (con %u divisori).\n", index, number, divisors);
	return EXIT_SUCCESS;
}

[U-Boat]

Per una volta mi lancio anch'io in questi contest.

Ho usato Java, ecco l'implentazione bruce force

public class TriangularNaive {

	private static int count_divisors(int n) {
		int i, count;
		
		for (i = 1, count = 0; i <= n; ++i)
			if ((n % i) == 0) ++count;
		
		return count;
	}
	
	private static int find_nth_triangular_number(int i){
		return (i * (i + 1)) / 2;
	}
	
	public static void main(String[] args) {
		int d = 200;
		int n = 1;
		int i = 1;
		
		long startTime = System.currentTimeMillis();
		
		System.out.println("Inizio la ricerca del primo numero triangolare che ha più di " + d + " divisori\n");
		while (count_divisors(n) < d) {
			n = find_nth_triangular_number(++i);
		}
		System.out.println("Si tratta del numero triangolare " + i + " che vale " + n + "\n");
		System.out.println("ci ho messo " + (System.currentTimeMillis() - startTime) + " millisecondi");

	}
}

Giusto per dare un riferimento sul pc che uso qui in ufficio per trovare il numero con più di 200 divisori ci mette circa 11 secondi.

La strada che ho seguito è stata quella di cercare un modo più veloce per sapere quanti divisori ha un numero intero positivo e una possibile risposta si trova questo indirizzo (ultimo punto in fondo).
L'idea è che se abbiamo la scomposizione in fattori primi, basta moltiplicare gli esponenti incrementati di uno e otteniamo il numero di divisori.

L'implementazione è la seguente (essendo stata fatta di corsissima lo stile lascia un po' a desiderare...)

public class TriangularSingleThread {
	
	private static int countDivisors(int n) {
		DivisorCounter divisorCounter = new DivisorCounter(n);
		
		divisorCounter.decompose();
		
		return divisorCounter.getDivisorCount();
	}
	
	private static int findNthTriangularNumber(int i){
		return (i * (i + 1)) / 2;
	}
	
	public static void main(String[] args) {
		int d = 1000;
		int n = 1;
		int i = 1;
		
		long startTime = System.currentTimeMillis();
		
		System.out.println("Inizio la ricerca del primo numero triangolare che ha più di " + d + " divisori\n");
		while (countDivisors(n) < d) {
			n = findNthTriangularNumber(++i);
		}
		System.out.println("Si tratta del numero triangolare " + i + " che vale " + n + "\n");
		System.out.println("ci ho messo " + (System.currentTimeMillis() - startTime) + " millisecondi");
		


	}
	
	private static class DivisorCounter {
		private int n;
		private int divisorCount;
		
		public DivisorCounter(int n) {
			this.n = n;			
			this.divisorCount = 1;
		}
		
		private void reduceUsingDivisor(int d) {
			if (n % d == 0) {
				int count = 0;
				
				while (n % d == 0) {
					count++;
					n = n / d;
				}
				
				divisorCount = divisorCount * (count + 1);
			}
		}
		
		public void decompose() {
			
			int d = 2;
			
			reduceUsingDivisor(d);
			
			d = 3;
			
			while (n != 1) {
				reduceUsingDivisor(d);
				d = d + 2;
			}
			
		}
		
		public int getDivisorCount() {		
			return divisorCount;
		}
		
	}
}

Con questa implementazione impiego 15 - 16 millesimi per trovare il numero con almeno 200 divisori e circa 1600-1700 per quello con 1000 (non è possibile andare oltre altrimenti si ha un overflow e bisogna passare ai long).

Ho cercato di spingermi oltre: ho provato istanziando 2 DivisorCounter e passando a ognuno solo uno dei due fattori della formula di Gauss, ma non ho ottenuto miglioramenti di sorta (proverò a lavorarci su un po')

[DanieleC88]

Originariamente inviato da U-Boat

La strada che ho seguito è stata quella di cercare un modo più veloce per sapere quanti divisori ha un numero intero positivo e una possibile risposta si trova questo indirizzo (ultimo punto in fondo).
L'idea è che se abbiamo la scomposizione in fattori primi, basta moltiplicare gli esponenti incrementati di uno e otteniamo il numero di divisori.

Ho provato anche io questa strada, ma cercando il primo numero triangolare con più di 1000 divisori non riesco a scendere comunque sotto i 12 secondi. L'algoritmo di VICIUS dev'essere particolarmente scaltro... ed ovviamente ci tiene sulle spine fino alla fine per riverlarcelo.

P.S.: VICIUS, ma l'agoritmo più efficiente con cui hai misurato i tempi è scritto in C o l'hai fatto in Ruby?

[VICIUS]

Originariamente inviato da U-Boat

La strada che ho seguito è stata quella di cercare un modo più veloce per sapere quanti divisori ha un numero intero positivo e una possibile risposta si trova questo indirizzo (ultimo punto in fondo).
L'idea è che se abbiamo la scomposizione in fattori primi, basta moltiplicare gli esponenti incrementati di uno e otteniamo il numero di divisori.

Abbiamo trovato la stessa pagina web
Questa è la mia implementazione di quel algoritmo.

class Numeric
  def multiple_of?(n)
    (self % n) == 0
  end
  
  def find_factors
    return [1] if self == 1
    
    n = self
    factors = Array.new
    
    while n.multiple_of? 2
      factors.push 2
      n /= 2
    end
    
    factor = 3
    while n > 1
      if n.multiple_of? factor
        while n.multiple_of? factor
          factors.push factor
          n /= factor
        end
      end
      factor += 2
    end
    factors
  end
  
  def count_divisors
    return 1 if self == 1
    
    factors = self.find_factors
    exponents = factors.uniq
    
    for i in 0...exponents.size
      exponents[i] = (factors.select { |f| f == exponents[i] }).size
    end
    
    exponents.map! { |e| e+1 }
    exponents.inject(1) { |e, tot| tot *= e }
  end
end

def contest_9(d)
  n = i = 1
  while n.count_divisors < d
    i += 1
    n += i
  end
  "Si tratta del numero triangolare #{i} che vale #{n}"
end

d = $*[0].to_i
puts "Inizio la ricerca del primo numero triangolare che ha piu di #{d} divisori"
puts contest_9(d)

[VICIUS]

Originariamente inviato da DanieleC88

Ho provato anche io questa strada, ma cercando il primo numero triangolare con più di 1000 divisori non riesco a scendere comunque sotto i 12 secondi. L'algoritmo di VICIUS dev'essere particolarmente scaltro... ed ovviamente ci tiene sulle spine fino alla fine per riverlarcelo.

P.S.: VICIUS, ma l'agoritmo più efficiente con cui hai misurato i tempi è scritto in C o l'hai fatto in Ruby?

È scritto in ruby e non è niente di speciale. Anzi è talmente semplice che è praticamente stupido

[DanieleC88]

Originariamente inviato da VICIUS

È scritto in ruby e non è niente di speciale. Anzi è talmente semplice che è praticamente stupido

Sarebbe questa versione che hai appena postato? Possibile che da te finisca in pochi millisecondi per un numero con oltre 1000 divisori?

[shinya]

Ho riscritto la mia versione in factor in un modo più idiomatico ed usando il metodo dei fattori primi che avete proposto...

! contest-9.factor - versione 2
USING: kernel sequences math.primes.factors ;

: count-divisors ( n -- x )
    1 swap group-factors [ second 1+ * ] each ;

: nth-triangular ( x -- n ) dup 1 + * 2 / ;

: solve-contest ( x n -- )
    dup nth-triangular count-divisors 
    pick > [ nip . ] [ 1+ solve-contest ] if ;

: contest-9 ( x -- ) 2 solve-contest ;

[gugoXX]

Variazione sul tema, ma sempre li' siamo
(C#)

Time to calculate all the Primes<10000000 : 650ms
First Triangle with more than 300 divisors:
2079th = 2162160 with 320 divisors, in 8ms

First Triangle with more than 1000 divisors:
41040th = 842161320 with 1024 divisors, in 276ms

[VICIUS]

Originariamente inviato da DanieleC88

Sarebbe questa versione che hai appena postato? Possibile che da te finisca in pochi millisecondi per un numero con oltre 1000 divisori?

No è un'altra. Molto più semplice.

[gugoXX]

I tempi vanno su...

Time to calculate all the Primes<10000000 : 695ms
First Triangle with more than 300 divisors:
2079th = 2162160 with 320 divisors, in 8ms

First Triangle with more than 1000 divisors:
41040th = 842161320 with 1024 divisors, in 245ms

First Triangle with more than 2000 divisors:
395199th = 78091322400 with 2880 divisors, in 1416ms

First Triangle with more than 10000 divisors:
73361375th = 2690945707626000 with 12800 divisors, in 41241ms

[gugoXX]

E poi scendono giu' (poco)
Basta, mi fermo qui

Time to calculate all the Primes<10000000 : 681ms
First Triangle with more than 300 divisors:
2079th = 2162160 with 320 divisors, in 8ms

First Triangle with more than 1000 divisors:
41040th = 842161320 with 1024 divisors, in 195ms

First Triangle with more than 2000 divisors:
395199th = 78091322400 with 2880 divisors, in 835ms

First Triangle with more than 10000 divisors:
73361375th = 2690945707626000 with 12800 divisors, in 29774ms

class Program
{
    static void Main(string[] args)
    {
        int MaxPrim = 10000000;
        Stopwatch sw = new Stopwatch();
        sw.Start();
        Decompos dc = new Decompos(MaxPrim);
        sw.Stop();
        Console.WriteLine("Time to calculate all the Primes<{0} : {1}ms", MaxPrim, sw.ElapsedMilliseconds);
        sw.Reset();
                    
        int ndiv;
        
        ndiv=300;
        Console.WriteLine("First Triangle with more than {0} divisors: ", ndiv);
        sw.Start();
        var res = dc.Striang(ndiv);
        sw.Stop();           
        Console.WriteLine("{0}th = {1} with {2} divisors, in {3}ms", res.TriangleOrder, res.Value, res.NDivisors, sw.ElapsedMilliseconds);
        Console.WriteLine();
        sw.Reset();

        ndiv = 1000;
        Console.WriteLine("First Triangle with more than {0} divisors: ", ndiv);
        sw.Start();
        res = dc.Striang(ndiv);
        sw.Stop();
        Console.WriteLine("{0}th = {1} with {2} divisors, in {3}ms", res.TriangleOrder, res.Value, res.NDivisors, sw.ElapsedMilliseconds);
        Console.WriteLine();
        sw.Reset();

        ndiv = 2000;
        Console.WriteLine("First Triangle with more than {0} divisors: ", ndiv);
        sw.Start();
        res = dc.Striang(ndiv);
        sw.Stop();
        Console.WriteLine("{0}th = {1} with {2} divisors, in {3}ms", res.TriangleOrder, res.Value, res.NDivisors, sw.ElapsedMilliseconds);
        Console.WriteLine();
        sw.Reset();

        ndiv = 10000;
        Console.WriteLine("First Triangle with more than {0} divisors: ", ndiv);
        sw.Start();
        res = dc.Striang(ndiv);
        sw.Stop();
        Console.WriteLine("{0}th = {1} with {2} divisors, in {3}ms", res.TriangleOrder, res.Value, res.NDivisors, sw.ElapsedMilliseconds);
        Console.WriteLine();
        sw.Reset();

        Console.ReadKey();
    }
}

public class Decompos
{
    public Decompos(int maxprim)
    {
        Primes = Eratostene.GetPrimes(maxprim).ToArray();
    }

    private int[] Primes;

    public class Result
    {
        public int TriangleOrder;
        public long Value;
        public int NDivisors;
        public Dictionary<long, int> Factors;
    }

    public Result Striang(int ndiv)
    {
        if (ndiv < 700) return StriangSingleCore(ndiv);
        return StriangDualCore(ndiv);
    }

    public Result StriangSingleCore(int ndiv)
    {
        long ctriang = 0;
        for (int t = 1; ; t++)
        {                
            ctriang += (long)t;
            var Fact = Factorize(ctriang, ndiv);
            if (Fact != null)
            {
                Fact.TriangleOrder = t;
                return Fact;
            }
        }
    }

    private Result Stepper(int findex, int lindex, int ndiv)
    {
        for (int u = findex; u < lindex; u++)
        {
            long longu = (long)u;
            long ctriang = (longu * (longu + 1)) / 2;
            var Fact = Factorize(ctriang, ndiv);
            if (Fact != null)
            {
                Fact.TriangleOrder = u;
                return Fact;
            }
        }
        return null;
    }

    private Result StriangDualCore(int ndiv)
    {
        int splitter = 5000;
        int bisplitter = splitter + splitter;
        for (int t = 1; ; t += bisplitter)
        {
            var r1 = Future<Result>.Create(() => Stepper(t, t + splitter, ndiv));
            var r2 = Future<Result>.Create(() => Stepper(t + splitter + 1, t + bisplitter, ndiv));

            Result r1res = r1.Value;
            if (r1res != null) return r1res;
            Result r2res = r2.Value;
            if (r2res != null) return r2res;                
        }
    }

    private Result Factorize(long input, int MinDivisors)
    {
        int PrimesLength = Primes.Length;
        long rem=input;
        Dictionary<long, int> ret = new Dictionary<long, int>();
        int IndexPrimo=0;
        int CurrentNDivisors = 1;
        do
        {
            if (IndexPrimo >= PrimesLength) throw new Exception("Primes Table too small");
            long Primo = Primes[IndexPrimo];

            int expected = (int)(rem / Primo) + 1;
            if ((CurrentNDivisors * expected ) < MinDivisors) return null;
            
            int Esponente = 0;
            while ((rem % Primo) == 0)
            {
                Esponente++;
                rem /= Primo;
            }

            if (Esponente != 0)
            {
                ret[Primo] = Esponente;
                CurrentNDivisors = CurrentNDivisors * (Esponente + 1);
            }
            IndexPrimo++;
        } while (rem != 1);

        if (CurrentNDivisors <= MinDivisors)
            return null;
                
        Result res = new Result();
        res.Value = input;
        res.NDivisors = CurrentNDivisors;
        res.Factors = ret;
        return res;
    }
}

public static class Eratostene
{
    public static IEnumerable<int> GetPrimes(int maxvalue)
    {            
        int[] Buff = new int[maxvalue];

        int mx = (int)Math.Sqrt((double)maxvalue)+1;
        int curprimo=2;
        do
        {
            yield return curprimo;
            for (int t = curprimo; t < maxvalue; t += curprimo)
            {
                Buff[t] = 1;
            }
            do
            {
                curprimo++;
            } while ((Buff[curprimo] != 0) && (curprimo < mx));
        } while (curprimo < mx);
        for (int t = curprimo; t < maxvalue; t++)
        {
            if (Buff[t] == 0) yield return t;
        }            
    }
}