ANALYSIS of ALGORITHMS Bulletin Board
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Open problem? Diophantine equations.



Hi folks,

I recently came across an interesting open problem. Someone asked for
a solution of the equation a+b+c+...=abc... in 2004 variables.

Newsgroup: es.ciencia.matematias
Thread:  Re: 2004 números
Message-Id:  <m3k729pdd0.fsf@professional.local>,
             <m28yio7hqu.fsf@linuxsexi.nasttg>

I wanted to find some test data in order to discover what the
solutions look like. My first attempt was very naive -- a couple of
nested loops to search the entire space of possible
solutions. Needless to say this only works for a small number of
variables, so I had to come up with something better.

I discovered a simple probabilistic algorithm that yields results
quickly. Start with a solution vector consisting entirely of
constants. Then compute all vectors that can be derived from the
current vector by adding or subtracting one from some element. Score
these vectors according to the formula LHS/RHS-1 and sort them
according to their scores. Then toss a die to select the new vector
from the set of branches, where the vector with the best score has,
say, probability 0.9, the next one 0.1 0.9 etc. This will either
converge to a solution or the algorithm halts if it has not found a
solution after a certain number of steps.

Note that the vectors form a tree where nodes differ in one position
from the parent node.

Open problem: given a diophantine equation that has a solution for
every n e.g. a+b+c...=abc..., compute the expected number of steps
until a solution is found. The analysis seems doable in this case.

I'd be thrilled to get some feedback on this! It would be wonderful to
see a formula for the expectation. BTW Hosam -- is this related to
Polya urns?

I am including a courtesy copy of the second article.

Best regards,

-- 
+------------------------------------------------------------+
| Marko Riedel, EDV Neue Arbeit gGmbH, mriedel@neuearbeit.de |
| http://www.geocities.com/markoriedelde/index.html          |
+------------------------------------------------------------+
From: mriedel@neuearbeit.de (Marko Riedel)
Subject: Re: 2004 números
Newsgroups: es.ciencia.matematicas
Date: 27 Feb 2004 14:51:37 +0100

Marko Riedel <mriedel@neuearbeit.de> writes:

> leonsotelo@wanadoo.es (Leon-Sotelo.) writes:
> 
> > ¿Existen 2004 enteros positivos cuya suma es igual a su producto?
> > 
> > Saludos.
> > León-Sotelo.
> 
> Hola a todos,
> 
> he escrito dos programas para estudiar este problema. Parece que para
> n bastante grande, la unica solucion es (2, n, 1^(n-2)), pero no puedo
> probarlo, tal vez es falso. Mis programas usan un algoritmo
> probabilistico. Comienzo con un vector de n elementos, todos igual a

Hola a todos,

me he dado cuenta que el algoritmo que he presentado puede usarse con
cualquier ecuacion diofantica, lo que es bastante interessante.

Estas son las ecuaciones que he anadido a mi programa:

1. a+b+c+ ... = abc ...
2. a_1^2 = a_2^3 + a_3^2 + a_4^3 + ...
3. a_1^3 = a_2^2 + a_3^3 + a_4^2 + ...
4. a_1^2 = a_2^2 + a_3^2 + a_4^2 + ... + a_n^2
5. a_1^5 = a_2^3 + a_3^3 + a_4^3 + ... + a_n^3

Bueno, pues supongo el proximo paso deberia ser el analisis de una de
estas ecuaciones para poder describir la estructura de las
soluciones. Veo en los datos que el numero tres parece bastante
facil. Alguien habra seguramente quien lograr encontrar una expresion
cerrada.

Aqui vienen los datos:

riedel@compaqtnew:~> ./dioeq.pl 15 sumprod
1 2 3
1 1 2 4
1 1 1 2 5
1 1 1 2 1 6
1 1 1 1 2 1 7
1 1 1 1 1 1 2 8
1 1 1 1 1 1 3 1 5
1 1 1 1 1 1 1 1 2 10
1 1 1 1 1 1 1 1 2 1 11
1 1 1 1 1 1 1 1 1 1 2 12
1 1 1 1 1 1 1 1 1 1 1 2 13
1 1 1 1 1 1 1 1 1 1 1 14 1 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 15
riedel@compaqtnew:~> ./dioeq.pl 15 alt23
14 3 13
6 2 1 3
10 2 8 3 1
15 2 5 4 8 4
18 4 8 4 2 4 8
25 2 12 5 7 5 7 5
29 4 7 5 10 5 9 6 9
34 4 11 5 9 5 12 6 8 6
41 6 14 6 10 6 3 6 14 6 10
53 7 13 7 10 7 13 7 13 7 12 7
57 7 14 7 13 7 14 7 11 7 14 8 12
70 8 1 8 14 8 16 8 14 8 15 8 15 9
75 7 12 8 17 8 15 8 17 9 10 9 16 9 16
riedel@compaqtnew:~> ./dioeq.pl 15 alt32
FAIL: 608 1 608
4 6 3 1
FAIL: 6 1 4 5 5
FAIL: 8 3 6 7 6 4
10 3 6 5 7 8 7
FAIL: 12 4 8 5 8 7 8 10
15 13 9 12 9 11 9 5 9
16 13 8 13 10 5 10 10 10 11
19 7 10 12 11 10 11 11 11 11 11
21 17 11 10 12 15 12 8 12 12 12 14
24 16 12 10 13 15 13 16 13 7 13 15 13
26 7 14 15 14 9 14 14 14 13 14 14 14 14
29 8 14 16 15 16 15 15 15 15 15 12 15 15 15
riedel@compaqtnew:~> ./dioeq.pl 15 p2p2
5 4 3
7 3 2 6
10 5 5 5 5
13 5 6 6 6 6
16 8 3 6 7 7 7
20 2 5 9 9 9 8 8
24 10 10 1 8 7 10 9 9
29 3 11 11 11 5 11 11 9 11
35 10 3 13 13 11 12 12 9 12 12
40 13 13 13 13 13 10 2 13 13 13 12
45 14 14 14 14 9 14 12 12 13 13 13 13
49 15 8 15 15 11 13 14 14 14 14 14 14 14
55 15 10 15 15 15 15 15 15 15 15 15 15 15 15
riedel@compaqtnew:~> ./dioeq.pl 15 p5p3
3 6 3
FAIL: 3 3 2 6
3 3 3 5 4
FAIL: 4 3 4 4 8 7
FAIL: 5 11 1 7 9 1 9
FAIL: 5 4 9 9 9 7 8 3
6 10 14 10 8 10 2 10 8
FAIL: 6 11 9 4 11 7 10 10 10 10
FAIL: 7 16 13 6 10 9 12 12 12 12 12
FAIL: 7 7 13 2 13 12 12 12 12 12 12 12
FAIL: 8 14 14 20 13 14 14 14 14 14 4 14 8
FAIL: 8 15 15 5 13 15 8 15 14 14 14 14 14 14
FAIL: 9 24 17 16 16 16 16 16 8 16 14 16 10 16 15
riedel@compaqtnew:~> ./dioeq.pl 15 p5p3
FAIL: 2 3 2
FAIL: 3 4 6 1
FAIL: 4 7 8 5 5
FAIL: 4 3 5 8 3 7
FAIL: 5 10 8 5 9 9 1
FAIL: 5 7 1 9 9 8 7 8
6 14 10 10 2 10 8 10 8
FAIL: 6 4 9 11 11 10 7 10 10 10
FAIL: 7 16 13 2 8 12 12 11 12 12 12
FAIL: 7 12 13 10 10 8 12 12 12 12 12 12
8 17 18 14 14 12 14 14 10 12 13 13 13
FAIL: 8 1 11 15 15 15 13 14 14 14 14 14 14 14
9 25 14 16 16 14 14 15 15 15 15 15 15 15 15
riedel@compaqtnew:~> exit


Sigue el nuevo programa.

#! /usr/bin/perl -w
#

my $upper = shift || 3;


sub val_sumprod {
  my ($s, $p) = (0, 1);

  foreach my $el (@{$_[0]}){
    $s += $el;
    $p *= $el;
  }

  return abs($s/$p-1);
}

sub val_alt23 {
  my $s = 0;
  my @els = @{$_[0]};
  my $first = shift @els;

  my $flag = 1;
  foreach my $el (@els){
    $s += $el*$el*($flag ? $el : 1);
    $flag = 1-$flag;
  }

  return abs($first*$first/$s-1);
}

sub val_alt32 {
  my $s = 0;
  my @els = @{$_[0]};
  my $first = shift @els;

  my $flag = 0;
  foreach my $el (@els){
    $s += $el*$el*($flag ? $el : 1);
    $flag = 1-$flag;
  }

  return abs($first*$first*$first/$s-1);
}


sub val_p2p2 {
  my $s = 0;
  my @els = @{$_[0]};
  my $first = shift @els;

  foreach my $el (@els){
    $s += $el*$el;
  }

  return abs($first*$first/$s-1);
}

sub val_p5p3 {
  my $s = 0;
  my @els = @{$_[0]};
  my $first = shift @els;

  foreach my $el (@els){
    $s += $el*$el*$el;
  }

  return abs($first**5/$s-1);
}


my %valprocs = 
    (sumprod => \&val_sumprod,
     alt23 => \&val_alt23,
     alt32 => \&val_alt32,
     p2p2 => \&val_p2p2,
     p5p3 => \&val_p5p3);

my $name = shift || 'sumprod';
my $valproc = $valprocs{$name};

if(not defined $valproc){
  print "$name not implemented\n";
  exit -1;
}

MAIN : {
  for(my $n=3; $n<=$upper; $n++){
    my ($count, $max) = (0, 16000);
    my ($vect, $v) = [ ($n) x $n ];

    while($v = &$valproc($vect), $v != 0.0){
      my @branches = ();
      $count++;

      for(my $pos=0; $pos<$n; $pos++){
	my $b =  [ @$vect ];
	$b->[$pos]++;
	push @branches, $b;

	if($vect->[$pos]>1){
	  $b = [ @$vect ];
	  $b->[$pos]--;
	  push @branches, $b;
	}
      }

      my %vals = map {
	$_, &$valproc($_)
      } @branches;

      my @sorted =
	sort { $vals{$a} <=> $vals{$b} }
	  @branches;

      my $pick;
      for($pick = 0; $pick <= $#sorted-1; $pick++){
	last if rand() < 0.9;
      }
    
      $vect = $sorted[$pick];
      last if $count==$max;
    }

    print 'FAIL: ' if $count==$max;
    print join(' ', @$vect) . "\n";
  }
}




Un saludo.

-- 
+------------------------------------------------------------+
| Marko Riedel, EDV Neue Arbeit gGmbH, mriedel@neuearbeit.de |
| http://www.geocities.com/markoriedelde/index.html          |
+------------------------------------------------------------+