Algoritmo de Schönhage-Strassen

El Algoritmo de Schonhage-Strassen es un algoritmo que consiste en la multiplicación de matrices. Es asintóticamente más rápido que el algoritmo de multiplicación de matrices estándar, pero más lento que el algoritmo más rápido conocido, y es útil en la práctica para matrices grandes.

Historia

El Algoritmo de Schonhage-Strassen fue el primer algoritmo en bajar en tiempo de ejecución de ε=1 , este fue descubierto en 1969 , por Volker Strassen de aquí el nombre del algoritmo. Este algoritmo logró alcanzar ε=0,808 en su tiempo de ejecución .

Desde ese momento muchos matemáticos han intentado hacer a este algoritmo más rápido para su mejora a nivel matemático y para ganar fama . Un gran salto fue obtenido en 1981 por Schönhage que logró ε=0,548 ,; Desde entonces la evolución de este es debida principalmente a Schönhage y a Strassen , ya que con sus investigaciones durante muchos años fueron reduciendo este tiempo y dando un gran paso para las matemáticas.

Análisis

El método clásico requiere n²(2n-1) operaciones aritméticas para multiplicar dos matrices (n x n). Strassen (1969) publicó un método que solo requería  operaciones. Se ha realizado mucho trabajo para intentar reducir el exponente 2,81, actualmente el mejor tiempo conocido es  operaciones (Virginia Vassilovska Williams). Hay todavía un largo camino por recorrer, pues la mejor cota inferior conocida es .

En lo que es la multiplicación es una de las operaciones que debemos repetir varias veces en enteros grandes, por lo que no está de más comprobar que la multiplicación por el método estándar de dos enteros binarios de longitudes k y precisa un máximo k' suma de un bit(Operaciones bit), por lo tanto el producto de dos números enteros con a lo sumo r cifras decimales requiere a lo más log² r operaciones. Es decir, es de orden O(log² r) y por lo tanto es una operación eficiente. Actualmente hay algoritmos de multiplicación más rápidos, el más conocido es el método de Karatsuba (cuya descripción es completamente accesible) que permite multiplicar dos enteros de longitud binaria menor o igual que k en O(k¹·59) operaciones bit. El método más rápido se debe a Schonhage y Strassen y su complejidad temporal es O(k log log log k).

Algoritmo  

Utilizando el algoritmo de Strassen definimos las siguientes matrices:

M1=(A1,2 - A2,2).(B2,1 + B2,2)

M2=(A1,1 + A2,2).(B1,1 + B2,2)

M3=(A1,1 - A2,1).(B1,1 + B2,1)

M4=(A1,1 + A1,2).B2,2

M5=A1,1 .(B1,2 -B2,2)

M6=A2,2.(B2,1 - B1,1)

M7=(A2,1 + A2,2).B1,1

Algoritmo en C

Ya que todos estos procesos son de una complejidad alta , habría que realizarlos en un computador. Por lo tanto habría que hacerlo con un lenguaje de programación como puede ser C, Java, Python... En este caso nos centramos en C++ por lo que el algoritmo es el siguiente: 

void  strassen ( vector <  vector < int >  >  Y Un ,  
              vector <  vector < int >  >  y B ,  
              vector <  vector < int >  >  y C ,  sin firmar int tam ); 
unsigned  int  nextPowerOfTwo ( int  n ); 
void  strassenR ( vector <  vector < int >  >  Y Un ,  
              vector <  vector < int >  >  y B ,  
              vector <  vector < int >  >  y C ,  
              int  tam ); 
void  suma ( vectorial <  vector < int >  >  Y Un ,  
         vector <  vector < int >  >  y B ,  
         vector <  vector < int >  >  y C ,  int  tam ); 
vacío  de resta ( vector <  vector < int >  >  Y Un ,  
              vector <  vector < int >  >  y B ,  
              vector <  vector < int >  >  y C ,  int  tam );

void  printMatrix ( vector <  vector < int >  >  matriz ,  int  n ); 
vacío  de lectura ( cadena de nombre de archivo ,  vector <  vector < int >  >  Y Un ,  vector <  vector < int >  >  y B );

void  ikjalgorithm ( vector <  vector < int >  >  A ,  
                                   vector <  vector < int >  >  B , 
                                   vector <  vector < int >  >  y C ,  int  n )  { 
    a  ( int  i  =  0 ;  i  <  n ;  i ++ )  { 
        a  ( int  k  =  0 ;  k  <  n ;  k ++ )  { 
            a  ( int  j  =  0 ;  j  <  n ;  j ++ )  { 
                C [ i ] [ j ]  + =  a [ i ] [ k ]  *  B [ k ] [ j ]; 
            } 
        } 
    } 
}

void  strassenR ( vector <  vector < int >  >  Y Un ,  
              vector <  vector < int >  >  y B ,  
              vector <  vector < int >  >  y C ,  int  tam )  { 
    si  ( tam  <=  leafsize )  { 
        ikjalgorithm ( A ,  B ,  C ,  tam ); 
        volver ; 
    }

    // Otros casos se tratan aquí: 
    otra cosa  { 
        int  newTam  =  TAM / 2 ; 
        vector < int >  interior  ( newTam ); 
        vector <  vector < int >  >  
            A11 ( newTam , interior ),  A12 ( newTam , interior ),  A21 ( newTam , interior ),  A22 ( newTam , interior ), 
            B11 ( newTam , interior ),  B12 ( newTam , interior ),  b21 ( newTam , interior ),  B22 ( newTam , interior ), 
              C11 ( newTam , interior ),  C12 ( newTam , interior ),  C21 ( newTam , interior ),  C22 ( newTam , interior ), 
            P1 ( newTam , interior ),  P2 ( newTam , interior ),  p3 ( newTam , interior ),  p4 ( newTam , interior ),  
            p5 ( newTam , interior ),  P6 ( newTam , interior ),  p7 ( newTam , interior ), 
            aResult ( newTam , interior ),  bResult ( newTam , interior );

        int  i ,  j ;

        // dividiendo las matrices en 4 submatrices: 
        para  ( i  =  0 ;  i  <  newTam ;  i ++ )  { 
            a  ( j  =  0 ;  j  <  newTam ;  j ++ )  { 
                a11 [ i ] [ j ]  =  A [ i ] [ j ]; 
                A12 [ i ] [ j ]  =  A [ i ] [ j  +  newTam ]; 
                a21 [ i ] [ j ]  =  A [ i  +  newTam ] [ j ]; 
                a22 [ i ] [ j ]  =  A [ i  +  newTam ] [ j  +  newTam ];

                b11 [ i ] [ j ]  =  B [ i ] [ j ]; 
                b12 [ i ] [ j ]  =  B [ i ] [ j  +  newTam ]; 
                b21 [ i ] [ j ]  =  B [ i  +  newTam ] [ j ]; 
                b22 [ i ] [ j ]  =  B [ i  +  newTam ] [ j  +  newTam ]; 
            } 
        }

        // Cálculo de P1 a P7:

        suma ( a11 ,  a22 ,  aResult ,  newTam );  // A11 + a22 
        suma ( b11 ,  b22 ,  bResult ,  newTam );  // B11 + B22 
        strassenR ( aResult ,  bResult ,  p1 ,  newTam );  // P1 = (a11 + a22) * (B11 + B22)

        suma ( a21 ,  a22 ,  aResult ,  newTam );  // A21 + a22 
        strassenR ( aResult ,  b11 ,  p2 ,  newTam );  // P2 = (a21 + a22) * (b11)

        restar ( b12 ,  b22 ,  bResult ,  newTam );  // B12 - B22 
        strassenR ( A11 ,  bResult ,  p3 ,  newTam );  // P3 = (a11) * (b12 - b22)

        restar ( b21 ,  b11 ,  bResult ,  newTam );  // B21 - b11 
        strassenR ( a22 ,  bResult ,  p4 ,  newTam );  // P4 = (a22) * (b21 - b11)

        suma ( a11 ,  a12 ,  aResult ,  newTam );  // A11 + a12 
        strassenR ( aResult ,  B22 ,  p5 ,  newTam );  // P5 = (a11 + a12) * (b22)   

        restar ( a21 ,  a11 ,  aResult ,  newTam );  // A21 - a11 
        suma ( b11 ,  b12 ,  bResult ,  newTam );  // + B11 b12 
        strassenR ( aResult ,  bResult ,  P6 ,  newTam );  // P6 = (A21-A11) * (B11 + B12)

        restar ( a12 ,  a22 ,  aResult ,  newTam );  // A12 - a22 
        suma ( b21 ,  b22 ,  bResult ,  newTam );  // B21 + B22 
        strassenR ( aResult ,  bResult ,  p7 ,  newTam );  // P7 = (A12-A22) * (B21 + B22)

        // Cálculo de c21, c21, c11 c22 e:

        sum ( p3 ,  p5 ,  c12 ,  newTam );  // C12 = p3 + p5 
        suma ( p2 ,  p4 ,  c21 ,  newTam );  // C21 = p2 + p4

        sum ( P1 ,  P4 ,  aResult ,  newTam );  // P1 + p4 
        suma ( aResult ,  p7 ,  bResult ,  newTam );  // P1 + P4 + p7 
        restar ( bResult ,  p5 ,  c11 ,  newTam );  // C11 = P1 + P4 - P5 + p7

        sum ( p1 ,  p3 ,  aResult ,  newTam );  // P1 + p3 
        suma ( aResult ,  P6 ,  bResult ,  newTam );  // P1 + p3 + P6 
        restar ( bResult ,  p2 ,  c22 ,  newTam );  // C22 = p1 + p3 - p2 + p6

        // La agrupación de los resultados obtenidos en una sola matriz: 
        para  ( i  =  0 ;  i  <  newTam  ;  i ++ )  { 
            a  ( j  =  0  ;  j  <  newTam  ;  j ++ )  { 
                C [ i ] [ j ]  =  c11 [ i ] [ j ]; 
                C [ i ] [ j  +  newTam ]  =  c12 [ i ] [ j ]; 
                C [ i  +  newTam ] [ j ]  =  c21 [ i ] [ j ]; 
                C [ i  +  newTam ] [ j  +  newTam ]  =  c22 [ i ] [ j ]; 
            } 
        } 
    } 
}

unsigned  int  nextPowerOfTwo ( int  n )  { 
    retorno  pow ( 2 ,  int ( ceil ( log2 ( n )))); 
}

void  strassen ( vector <  vector < int >  >  Y Un ,  
              vector <  vector < int >  >  y B ,  
              vector <  vector < int >  >  y C ,  sin firmar int n )  { 
    // unsigned int n = tam; 
    unsigned  int  m  =  nextPowerOfTwo ( n ); 
    vector < int >  interior ( m ); 
    vector <  vector < int >  >  APrep ( m ,  interior ),  BPrep ( m ,  interior ),  CPrep ( m ,  interior );

    para ( sin firmar int i = 0 ;  i < n ;  i ++ )  { 
        a  ( unsigned  int  j = 0 ;  j < n ;  j ++ )  { 
            APrep [ i ] [ j ]  =  A [ i ] [ j ]; 
            BPrep [ i ] [ j ]  =  B [ i ] [ j ]; 
        } 
    }

    strassenR ( APrep ,  BPrep ,  CPrep ,  m ); 
    para ( sin firmar int i = 0 ;  i < n ;  i ++ )  { 
        a  ( unsigned  int  j = 0 ;  j < n ;  j ++ )  { 
            C [ i ] [ j ]  =  CPrep [ i ] [ j ]; 
        } 
    } 
}

void  suma ( vectorial <  vector < int >  >  Y Un ,  
         vector <  vector < int >  >  y B ,  
         vector <  vector < int >  >  y C ,  int  tam )  { 
    int  i ,  j ;

    para  ( i  =  0 ;  i  <  tam ;  i ++ )  { 
        a  ( j  =  0 ;  j  <  tam ;  j ++ )  { 
            C [ i ] [ j ]  =  A [ i ] [ j ]  +  B [ i ] [ j ]; 
        } 
    } 
}

vacío  de resta ( vector <  vector < int >  >  Y Un ,  
              vector <  vector < int >  >  y B ,  
              vector <  vector < int >  >  y C ,  int  tam )  { 
    int  i ,  j ;

    para  ( i  =  0 ;  i  <  tam ;  i ++ )  { 
        a  ( j  =  0 ;  j  <  tam ;  j ++ )  { 
            C [ i ] [ j ]  =  A [ i ] [ j ]  -  B [ i ] [ j ]; 
        } 
    }    
}

int  getMatrixSize ( string  nombre )  { 
    string  línea ; 
    ifstream  archivoentrada ; 
    archivoentrada . abierta  ( nombre de archivo . c_str ()); 
    getline ( archivoentrada ,  línea ); 
    volver  recuento ( line . iniciar (),  la línea . final (),  '\ t' )  +  1 ; 
}

vacío  de lectura ( cadena de nombre de archivo ,  vector <  vector < int >  >  Y Un ,  vector <  vector < int >  >  y B )  { 
    string  línea ; 
    ARCHIVO *  matrixfile  =  freopen ( nombre de archivo . C_str (),  "r" ,  la entrada estándar );

    si  ( matrixfile  ==  0 )  { 
        cerr  <<  "No se pudo leer el archivo"  <<  nombre de fichero  <<  endl ; 
        volver ; 
    }

    int  i  =  0 ,  j ,  una ; 
    mientras que  ( getline ( cin ,  línea )  &&  ! line . vacía ())  { 
        istringstream  ISS ( línea ); 
        j  =  0 ; 
        mientras que  ( ISS  >>  una )  { 
            A [ i ] [ j ]  =  una ; 
            j ++ ; 
        } 
        I ++ ; 
    }

    i  =  0 ; 
    mientras que  ( getline ( cin ,  línea ))  { 
        istringstream  ISS ( línea ); 
        j  =  0 ; 
        mientras que  ( ISS  >>  una )  { 
            B [ i ] [ j ]  =  una ; 
            j ++ ; 
        } 
        I ++ ; 
    }

    fclose  ( matrixfile ); 
}

anulará  printMatrix ( vector <  vector < int >  >  matriz ,  int  n )  { 
    a  ( int  i = 0 ;  i  <  n ;  i ++ )  { 
        a  ( int  j = 0 ;  j  <  n ;  j ++ )  { 
            si  ( j  ! =  0 )  { 
                cout  <<  " \ t " ; 
            } 
            Cout  <<  matriz [ i ] [ j ]; 
        } 
        Cout  <<  endl ; 
    } 
}

int  principal  ( int  argc ,  char *  argv [])  { 
    string  nombre ; 
    si  ( argc  <  3 )  { 
        nombre de archivo  =  "2000.in" ; 
    }  Demás  { 
        nombre de archivo  =  argv [ 2 ]; 
    }

    si  ( argc  <  5 )  { 
        leafsize  =  16 ; 
    }  Demás  { 
        leafsize  =  atoi ( argv [ 4 ]); 
    }

    int  n  =  getMatrixSize ( nombre de archivo ); 
    vector < int >  interior  ( n ); 
    vector <  vector < int >  >  A ( n ,  interior ),  B ( n ,  interior ),  C ( n ,  interior ); 
    leer  ( nombre de archivo ,  A ,  B ); 
    strassen ( A ,  B ,  C ,  N ); 
    printMatrix ( C ,  n ); 
    volver  0 ; 
}

Bibliografía 
 http://www.algoritmia.net/articles.php?id=34 http://es.wikipedia.org/wiki/Algoritmo_de_Strassen
 http://decsai.ugr.es/~castro/CA/node5.html  
 https://web.archive.org/web/20100921151810/http://jair.lab.fi.uva.es/~fdelgado/cripto/resumen.pdf 
https://martin-thoma.com/strassen-algorithm-in-python-java-cpp/

Algoritmos de precisión arbitraria