Saltar al contenido

Encontrar el máximo común divisor de una matriz en MATLAB

Ya no tienes que buscar más en otros sitios ya que has llegado al lugar correcto, poseemos la solución que buscas y sin problemas.

Solución:

Como no le gustan los bucles, ¿qué hay de las funciones recursivas?

iif = @(varargin) varargin2 * find([varargin1:2:end], 1, 'first')();
[email protected](v,gcdr) iif(length(v)==1,v, ...
                     v(1)==1,1, ...
                     length(v)==2,@()gcd(v(1),v(2)), ...
                     true,@()gcdr([gcd(v(1),v(2)),v(3:end)],gcdr));
[email protected](v)gcdrec(v(:)',gcdrec);

A=[0,0,0; 2,4,2;-2,0,8];
divisor=mygcd(A);
A=A/divisor;

La primera función iif definirá una construcción condicional en línea. Esto permite definir una función recursiva, gcdrec, para encontrar el máximo común divisor de su array. Esta iif funciona así: prueba si el primer argumento es true, si es así, devuelve el segundo argumento. De lo contrario, prueba el tercer argumento, y si esa estrue, luego devuelve el cuarto, y así sucesivamente. Debe proteger las funciones recursivas y, a veces, otras cantidades que aparecen en su interior con @(), de lo contrario, puede obtener errores.

Utilizando iif la función recursiva gcdrec funciona así:

  • si el vector de entrada es un escalar, lo devuelve
  • de lo contrario, si el primer componente del vector es 1, no hay posibilidad de recuperación, por lo que devuelve 1 (permite un retorno rápido para matrices grandes)
  • de lo contrario, si el vector de entrada es de longitud 2, devuelve el máximo común divisor a través de gcd
  • de lo contrario, se llama a sí mismo con un vector abreviado, en el que los dos primeros elementos se sustituyen por su máximo común divisor.

La función mygcd es solo una interfaz para mayor comodidad.

Debería ser bastante rápido, y supongo que solo la profundidad de recursividad podría ser un problema para problemas muy grandes. Hice una verificación de tiempo rápida para comparar con la versión en bucle de @Adriaan, usando A=randi(100,N,N)-50, con N=100, N=1000 y N=5000 y tic/toc.

  1. N=100:
    • bucle de 0,008 segundos
    • recursivo 0,002 segundos
  2. N=1000:
    • bucle 0,46 segundos
    • recursivo 0,04 segundos
  3. N=5000:
    • bucle 11,8 segundos
    • recursivo 0,6 segundos

Actualizar: Lo interesante es que la única razón por la que no supere el límite de recursividad (que por defecto es 500) es que mis datos no tenían un divisor común. Establecer una matriz aleatoria y doblarla conducirá a alcanzar el límite de recursividad ya por N=100. Entonces, para matrices grandes, esto no funcionará. Por otra parte, para matrices pequeñas, la solución de @ Adriaan está perfectamente bien.

También intenté reescribirlo a la mitad del vector de entrada en cada paso recursivo: esto de hecho resuelve el problema del límite de recursividad, pero es muy lento (2 segundos para N=100, 261 segundos para N=1000). Puede haber un término medio en algún lugar, donde el tamaño de la matriz es grande (más o menos) y el tiempo de ejecución no es tan malo, pero aún no lo he encontrado.

 A = [0,0,0; 2,4,2;-2,0,8];
 B = 1;
 kk = max(abs(A(:))); % start at the end
 while B~=0 && kk>=0
     tmp = mod(A,kk);
     B = sum(tmp(:));
     kk = kk - 1;
 end
 kk = kk+1;

Probablemente esta no sea la forma más rápida, pero lo hará por ahora. Lo que hice aquí fue inicializar algún contador, B, para almacenar la suma de todos los elementos en su matriz después de tomar la mod. los kk es solo un contador que pasa por números enteros. mod(A,kk) calcula el módulo después de la división para cada elemento en A. Por lo tanto, si todos sus elementos son totalmente divisibles por 2, devolverá un 0 para cada elemento. sum(tmp(:)) luego hace una sola columna a partir de la matriz módulo, que se suma para obtener algún número. Si y solo si ese número es 0, hay un divisor común, ya que entonces todos los elementos en A son totalmente divisibles por kk. Tan pronto como eso suceda, su bucle se detendrá y su divisor común es el número en kk. Ya que kk disminuye cada recuento, en realidad es un valor demasiado bajo, por lo que se agrega uno.

Nota: acabo de editar el bucle para ejecutar hacia atrás ya que está buscando el CD más grande, no el más pequeño. Si tuvieras una matriz como [4,8;16,8] se detendría en 2, no 4. Disculpas por eso, esto funciona ahora, aunque las otras dos soluciones aquí son mucho más rápidas.

Finalmente, dividir matrices se puede hacer así:

divided_matrix = A/kk;

De acuerdo, ¡a mí tampoco me gustan los bucles! Vamos a matarlos –

unqA = unique(abs(A(A~=0))).';             %//'
col_extent = [2:max(unqA)]';               %//'
B = repmat(col_extent,1,numel(unqA));
B(bsxfun(@gt,col_extent,unqA)) = 0;
divisor = find(all(bsxfun(@times,bsxfun(@rem,unqA,B)==0,B),2),1,'first');
if isempty(divisor)
    out = A;
else
    out = A/divisor;
end

Ejecuciones de muestra

Caso 1:

A =
     0     0     0
     2     4     2
    -2     0     8
divisor =
     2
out =
     0     0     0
     1     2     1
    -1     0     4

Caso # 2:

A =
     0     3     0
     5     7     6
    -5     0    21
divisor =
     1
out =
     0     3     0
     5     7     6
    -5     0    21

Te invitamos a patrocinar nuestra ocupación mostrando un comentario o valorándolo te damos la bienvenida.

¡Haz clic para puntuar esta entrada!
(Votos: 0 Promedio: 0)



Utiliza Nuestro Buscador

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *