/* GccFactorial.c

 *

 * Created: 31.05.2019

 * Author : EH

 */

// 31.05.2015 version based on original QuickC version,

//           but added "if zd[j] == 0", "if zd[j] == 1" and "if fakul[k] == 0", which speeds up a lot.

//           Unrolling "for (j = 0; j <= 3; j++)" does not make any difference.

 

#define F_CPU 18000000UL   // 18 MHz

#include <avr/io.h>

#include <util/delay.h>

 

#define LED_ON      PORTC |= (1 << PORTC0)

#define LED_OFF     PORTC &= ~(1 << PORTC0)

#define ExpLen      2568         // erwartete Anzahl der Ergebnisstellen von 1000!

 

/**********************************************************************

**      Variablen-Deklarationen

**********************************************************************/

 

int z;                            // z > 1 --> z! ist zu berechnen

uint8_t fakul[ExpLen+3];          // Fakultät = Ergebnis, ExpLen+3 Stellen Genauigkeit

 

uint8_t tp[4][ExpLen+3];          // 4 Teilprodukte

uint8_t zd[4];                    // 4-stellige Zahl in Digit-Darstellung

 

int zi;                           // Zwischenwert

int a, b;                         // Zwischenwerte bei der Fakultät-Berechnung

int i, j, k;                      // Schleifen-Zähler

      

/**********************************************************************

**      Berechnung der Fakultät

***********************************************************************/

 

int factorial(int z)

{

       fakul[1] = 1;                     // fakul = 1 setzen, so geht's los

 

       for (i = 2; i <= ExpLen+2; i++)         // restl. Ziffern von fakul auf Null setzen

       {

             fakul[i] = 0;

       }

       tp[1][1] = 0;                     // tp(1) rechts auf   0 setzen

       tp[2][1] = 0;                     // tp(2) rechts auf  00 setzen

       tp[2][2] = 0;                     // tp(2) rechts auf  00 setzen

       tp[3][1] = 0;                     // tp(3) rechts auf 000 setzen

       tp[3][2] = 0;                     // tp(3) rechts auf 000 setzen

       tp[3][3] = 0;                     // tp(3) rechts auf 000 setzen

 

       tp[0][ExpLen] = 0;                // tp(0) links auf 000 setzen

       tp[0][ExpLen+1] = 0;              // tp(0) links auf 000 setzen

       tp[0][ExpLen+2] = 0;              // tp(0) links auf 000 setzen

       tp[1][ExpLen+1] = 0;              // tp(1) links auf 00  setzen

       tp[1][ExpLen+2] = 0;              // tp(1) links auf 00  setzen

       tp[2][ExpLen+2] = 0;              // tp(2) links auf 0   setzen

      

       // alle anderen tp werden im Programm definiert

      

       // ab hier z!-Berechnung

 

       for (i = 2; i <= z; i++)                       // 1 * 2 * 3 * ... z

       {

             zi = i;                                 // Zwischenwert

             for (j = 0; j <= 3; j++)

             {

                    zd[j] = zi % 10;                  // mod10 = Einerstelle

                    zi = zi / 10;                     // nächste Stelle

             }

             for (j = 0; j <= 3; j++)                // Teilprodukte erstellen

             {

                    if (zd[j] == 0)                   // erspart viele Multiplikationen und viel Zeit

                    {

                           for (k = 1; k<= ExpLen+2; k++)

                           tp[j][k] = 0;              // aber tp dafür löschen

                    }

                    else if (zd[j] == 1)                    // erspart viele Multiplikationen

                    {

                           for (k = 1; k<= ExpLen-3; k++)

                           {

                                  tp[j][k + j] = fakul[k];   // dafür ganze Zeile kopieren

                           }

                           tp[j][ExpLen-2 + j] = 0;

                           tp[j][ExpLen-1 + j] = 0;

                    }

                    else                    {

                           b = 0;                            // Übertrag

                           for (k = 1; k<= ExpLen-3; k++)

                           {

                                  if (fakul[k] == 0)         // auch das bringt viel

                                  {

                                        tp[j][k + j] = b;   // b immer <= 3

                                        b = 0;              // kein weiterer Übertrag

                                  }

                                  else if (fakul[k] == 1)    // diese Abfrage bringt kaum Zeitgewinn

                                  {

                                        a = b + zd[j];

                                        tp[j][k + j] = a % 10;     // Einer

                                        a = a / 10;

                                        b = a % 10;         // Zehner = Übertrag

                                  }

                                  else

                                  {

                                        a = b + fakul[k] * zd[j];

                                        tp[j][k + j] = a % 10;     // Einer

                                        a = a / 10;

                                        b = a % 10;         // Zehner = Übertrag

                                  }

                           }

                           tp[j][ExpLen-2 + j] = b;

                           tp[j][ExpLen-1 + j] = 0;

                    }

             }

             b = 0;                                         // Übertrag

             for (j = 1; j<= ExpLen+1; j++)                 // 4 Teilprodukte addieren

             {

                    a = b + tp[0][j] + tp[1][j] + tp[2][j] + tp[3][j];

                    fakul[j] = a % 10;                      // Einer

                    a = a / 10;

                    b = a % 10;                             // Zehner = Übertrag

             }

             fakul[ExpLen+2] = b;                           // um ganz korrekt zu sein

       }

       return 0;

}

 

int main(void)

{

       DDRC |= (1 << DDC0);              // Port C, Bit 0 ist ein LED-Ausgang

       LED_ON;

       DDRA = 0xff;               // Port A ist ein Ausgang

       DDRD = 0xff;               // Port D ist ein Ausgang

       PORTA = 0x00;

       PORTD = 0x00;

       _delay_ms(500);

      

    while (1)

    {

       LED_OFF;

       factorial(1000);

       LED_ON;

       PORTA = 16*fakul[ExpLen] + fakul[ExpLen-1];    // erste beiden Digits (BCD) von 1000!,

                                                      // sollte 40 sein, dauert 2:16 min

       PORTD = 16*fakul[ExpLen-2] + fakul[ExpLen-3];  // nächste beiden Digits (BCD) von 1000!,

                                                      // sollte 23 sein         

       _delay_ms(500);     // kurzer HIGH-Puls (500 ms) für die LED zum Bestimmen der Laufzeit             

    }

}