/*
 * Program: The Multi-Layer PerceptronTrainer
 * Programmer: Robert John Morton UK-YE572246C
 * Programming Language: C
 * Started:  December 1997

                     pointer to input pattern array
                     |          pointer to 'correct' output pattern array
                     |          |       Shift factor corresponding
                     |          |       |  to weight gain term η.
                     |          |       |      Shift factor corresponding 
                     |          |       |      |  to Momentum Factor α.
                     |          |       |      |  */
void mlptrain(short *pi, short *pt, int h, int a) {
  short
    E1[N1],E2[N2],E3[N2],    // output errors for each layer
    *E[] = {NULL,E1,E2,E3},  // array of ptrs to above arrays

    M10[N0],M11[N0],M12[N0],M13[N0], // Layer 1 input δw
    M20[N1],M21[N1],M22[N1],         // Layer 2 input δw
    M30[N1],M31[N1],                 // Layer 3 input δw
    *M1[] = {M10, M11, M12, M13},    // Access to L1 δw
    *M2[] = {M20, M21, M22},         // Access to L2 δw
    *M3[] = {M30, M31},              // Access to L3 δw
    **M[] = {NULL, M1, M2, M3};      // Access to all δw

  int nl;       // layer number
    L[0] = pi;  // points to start of network inputs array
    h += 15;    // shift factor to multiply by η / R

  for(nl = NAL; nl > 0; nl−−) {  // for each layer of network
    short 
      **ppw = *(W + nl),   // ptr to this layer's weights
      **ppm = *(W + nl),   // ptr to layer's delta-weights
      *pe = *(E + nl),     // ptr to layer's output errors
      *po = *(L + nl);     // ptr to this layer's outputs
    int 
      nn, ni,              // neuron number, input number
      NN = *(N + nl),      // number of neurons in layer
      NI = *(N − 1 + nl);  // number of inputs to layer

    /* If processing the output layer, prime each element 
       of the error array with −(t[j] − o[j]). */
    if(nl == NAL)
      for (nn = 0; nn < NN; nn++)
        *(pe + nn) = *(po + nn) − *(pt + nn);

    pi = *(L + nl − 1);  // pointer to start of layer's inputs

    //Compute the output error for each neuron in this layer.
    for(nn = 0; nn < NN; nn++) {
      short 
        *pw = *(ppw + nn),  // ptr to neuron's first weight 
        *pm = *(ppm + nn);  // ptr to neuron's 1st delta-weight
      long 
        o = *(po + nn),     // this neuron's output signal
        e = (((R + o) * (R − o)) >> 15) * *(pe + nn) >> 13; 

      if(e > R) e = R; if(e < −R) e = −R;
      *(pe + nn) = e;  // ∂F/∂a = ∂o/∂a * last time's summation

      for(ni = 0; ni < NI; ni++)  // adjust each input weight
        *(pw + ni) += (
                        *(pm +ni) = (*(pm + ni) >> a) 
                                  − (((e * *(pi + ni)) / NI) >> h)
                      );
    }
    if(nl > 1) {                    // if not yet reached 1st active layer 
      short *ps = *(E + nl − 1);    // ptr to prev layer's output errors
      for(ni = 0; ni < NI; ni++) {  // for each input weight to this layer
        long Hi = 0, Lo = 0;        // see mlp() to explan the following code
        for(nn = 0; nn < NN; nn++) {
          long P = (long)*(pe + nn) * *(*(ppw + nn) + ni);
          Hi += P >> 16; Lo += P & 0xFFFF;
        }
        *(ps + ni) = ((Hi << 1) + (Lo >> 15)) / NN;
      }
    }   // ... prime the previous layer's error array elements 
  }     //     with this layer's error * weight summations. 
}