Posts Tagged ‘liczby armstronga’

PostHeaderIcon Liczby Armstronga

Liczbą Armstronga (inaczej narcystyczną) nazywamy n-cyfrową liczbę naturalną, która jest sumą swoich cyfr podniesionych do potęgi n. Udowodniono, że takich liczb jest 88, przy czym najwyższa z nich ma 39 cyfr. Oczywiście algorytm do obliczania tak dużych liczb nie może operować na powszechnie znanych typach danych – przykładowo liczba typu int w języku C może być maksymalnie równa w przybliżeniu 4 miliardy. Aby poradzić sobie z tym problemem musimy stworzyć własny typ danych i zdefiniować dla niego podstawowe operacje (dla naszych celów przydadzą się jedynie operacje dodawania i mnożenia).

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>

#define MAX_DIGTS 40
#define MAX_DIGTS2 (MAX_DIGTS+1)*(MAX_DIGTS+1)


uint_fast32_t powtabl[10*MAX_DIGTS2];
uint_fast32_t tmp[MAX_DIGTS];
uint_fast32_t maxnumber[MAX_DIGTS];
uint_fast32_t liczba[MAX_DIGTS];
uint_fast32_t digt;
uint_fast32_t ifrom[10], ilosc[10];
uint_fast32_t minnumber[11*MAX_DIGTS];
uint_fast32_t tmpx;
uint_fast32_t digtscnt[10];
uint_fast32_t flag = 1;
uint_fast32_t iend;
uint_fast32_t found = 0;
uint_fast32_t w;


// aa < bb
inline
//static
void addl(uint_fast32_t* c, const uint_fast32_t* a, const uint_fast32_t* b)
{
    uint_fast32_t i;

    for (i = 1; i <= *a; i++)
        c[i] = a[i] + b[i];

    for (i = *a+1; i <= *b; i++)
        c[i] = b[i];

    for (i = 1; i < *a; i++)
        if (c[i] > 9)
        {
            c[i] += -10;
            ++c[i+1];
        }

    for (i = *a; i < *b; i++)
    {
        if (c[i] > 9)
        {
            c[i] += -10;
            ++c[i+1];
        }
        else
        {
            *c = *b;
            return;
        }
    }

    if (c[*b] > 9)
    {
        c[*b] += -10;
        *c = *b+1;
        c[*c] = 1;
        return;
    }
    *c = *b;
}

inline
//static
void muli(uint_fast32_t *c, const uint_fast32_t *b, const uint_fast32_t a)
{
    uint_fast32_t i, tmpc;

    for (i = 1; i <= *b; i++)
        c[i] = b[i] * a;

    i = 0;
    while (++i < *b)
        if (c[i] > 9)
        {
            tmpc = 0.1*c[i];
            c[i+1] += tmpc;
            c[i] += -(tmpc << 3);
            c[i] += -(tmpc << 1);
        }

    if (c[i] > 9)
    {
        c[i+1] = 0.1*c[i];
        c[i] += -(c[i+1] << 3);
        c[i] += -(c[i+1] << 1);
        if (c[++i] > 9)
        {
            c[i+1] = 0.1*c[i];
            c[i] += -(c[i+1] << 3);
            c[i] += -(c[i+1] << 1);    // dwa zera z przodu
        }
    }

    *c = i;
}


//inline
static
uint_fast32_t *powll(uint_fast32_t a, uint_fast32_t b)
{
    uint_fast32_t tmp[MAX_DIGTS], i, *xx;

    if (powtabl[a*MAX_DIGTS2+b*(MAX_DIGTS)] == 0)
    {
        xx = powll(a,b-1);
        for (i = 0; i < MAX_DIGTS; i++)
            *(tmp+i) = *(xx+i);
           
        while (tmp[*tmp] == 0 && *tmp >= 2)
            (*tmp)--;
           
        muli(tmp, tmp, a);
       
        for (i = 0; i < MAX_DIGTS; i++)
            powtabl[a*MAX_DIGTS2+b*(MAX_DIGTS+1)+i] = tmp[i];
           
        while (powtabl[a*MAX_DIGTS2+b*(MAX_DIGTS+1)+powtabl[a*MAX_DIGTS2+b*(MAX_DIGTS+1)]] == 0 &&
                powtabl[a*MAX_DIGTS2+b*(MAX_DIGTS+1)] >= 2)
            powtabl[a*MAX_DIGTS2+b*(MAX_DIGTS+1)]--;
    }
    return powtabl+a*MAX_DIGTS2+b*(MAX_DIGTS+1);
}


int main()
{
    int_fast32_t i, ito[10], it[10];
    uint_fast32_t n, j;

    memset(powtabl, 0, 10*MAX_DIGTS2*sizeof(uint_fast32_t));

    for (i = 9; i >= 0; i--)
    {
        //liczymy pierwsza i zerowa potege
        powtabl[i*MAX_DIGTS2] = 1;
        powtabl[i*MAX_DIGTS2+1] = 1;
        powtabl[i*MAX_DIGTS2+MAX_DIGTS+1] = 1;
        powtabl[i*MAX_DIGTS2+MAX_DIGTS+2] = i;
        ito[i] = 0;
    }
   
    //liczymy pozostale potegi (1 do 39)
    for (i = 1; i < 10; i++)
        powll(i,MAX_DIGTS-1);

    minnumber[10*(MAX_DIGTS+1)] = 1;
    minnumber[10*(MAX_DIGTS+2)] = 0;

    for (n = 1; n <= MAX_DIGTS; n++)
    {
        digt = 9;
        ifrom[9] = it[9] = n;

        while (digt < 10)
            if (it[digt] >= ito[digt])
            {
                if (digt == 1)
                {
                    liczba[0] = liczba[1] = 1;

                    muli(liczba, liczba, *(it+1));
                    if (*liczba >= *(minnumber+2*(MAX_DIGTS+1)))
                        addl(tmp, minnumber+2*(MAX_DIGTS+1), liczba);
                    else
                        addl(tmp, liczba, minnumber+2*(MAX_DIGTS+1));
                    memcpy(liczba, tmp, (n+1)*sizeof(uint_fast32_t));
                    while (liczba[*liczba] == 0 && *liczba > 1)
                    {
                        (*liczba)--;
                    }

                    if (*liczba == n)
                    {
                        //memset(ilosc, 0, 10*sizeof(int_fast32_t));
                        for (i = 9; i >= 0; i--)
                            *(ilosc+i) = 0;
                        for (i = *liczba; i > 0; i--)
                            ilosc[liczba[i]]++;
                        if (((ilosc[1] - it[1]) | (ilosc[2] - it[2]) | (ilosc[3] - it[3]) | (ilosc[4] - it[4]) | (ilosc[5] - it[5]) | (ilosc[6] - it[6]) | (ilosc[7] - it[7]) | (ilosc[8] - it[8]) | (ilosc[9] - it[9])) == 0)
                        {
                            printf("%d/89: ", ++found);
                            for (i = *liczba; i > 0; i--) printf("%d", liczba[i]);
                            printf(" (%d)\n", n);
                        }
                    }
                    it[digt]--;
                    continue;
                }

                else if (digt != 9)
                {
                    muli(tmp, powtabl+digt*MAX_DIGTS2+n*(MAX_DIGTS+1), it[digt]);
                    if (*(minnumber+(MAX_DIGTS+1)*(digt+1)) >= *tmp)
                        addl(minnumber+(MAX_DIGTS+1)*digt, tmp, minnumber+(MAX_DIGTS+1)*(digt+1));
                    else addl(minnumber+(MAX_DIGTS+1)*digt, minnumber+(MAX_DIGTS+1)*(digt+1), tmp);
                }

                else muli(minnumber+(MAX_DIGTS+1)*9, powtabl+9*MAX_DIGTS2+n*(MAX_DIGTS+1), it[9]);

                ifrom[digt-1] = ifrom[digt]-it[digt];
                muli(maxnumber, powtabl+(digt-1)*MAX_DIGTS2+n*(MAX_DIGTS+1), ifrom[digt-1]);

                tmpx = (MAX_DIGTS+1)*digt;
                if (*maxnumber < *(minnumber+tmpx))
                    addl(maxnumber, maxnumber, minnumber+tmpx);
                else
                    addl(maxnumber, minnumber+tmpx, maxnumber);

                // ilosc najbardziej znaczacych, wspolnych cyfr
                if (*maxnumber < n)
                {
                    --it[++digt];
                }
                else if (minnumber[tmpx] > n)
                    --it[digt];
                else
                {
                    j = n;
                    while (minnumber[tmpx+j] == maxnumber[j] && j > 0)
                    {
                        --j;
                    }
                    iend = minnumber[tmpx];

                    for (i = 9; i >= 0; i--)
                        digtscnt[i] = 0;

                    for (i = n-j; i > 0; i--)
                        ++digtscnt[minnumber[tmpx+i+iend-n+j]];
                       
                    for (i = 9; i >= digt; i--)
                        flag = flag & (digtscnt[i] <= it[i]);

                    for (i = digt-1; i > 0; i--)
                        ito[i] = digtscnt[i];

                    if (flag != 1)
                    {
                        --it[digt];
                        flag = 1;
                    }
                    else
                    {
                        --digt;
                        it[digt] = ifrom[digt];
                    }
                }
            }
            else
            {
                --it[++digt];
            }
    }

    return 0;
}

Kod zapisujemy pod nazwą big-number.c i kompilujemy pod Linuksem poleceniem:

gcc -Wall -O2 -ffast-math -funroll-loops -fomit-frame-pointer big-number.c -o numb

Możemy oczywiście stworzyć plik Makefile, aby zaoszczędzić sobie wpisywania lub pamiętania wszystkich flag przy następnych kompilacjach. Ma on postać:

numb: big-number.c
gcc -Wall -O2 -ffast-math -funroll-loops -fomit-frame-pointer big-number.c -o numb