Поразвлечься :)

/*
 * Filename: automorph.cpp
 * Created:  Thu Oct 20 2011
 * Поиск больших автоморфных чисел.
 * См. "Квант", 1973, № 7, стр. 55
 */

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

int LAST  = 1000;  // Сколько цифирей надо
int PRECISION;     // Сколько выделять int под большое число

// Чтоб не вычислять - берем готовое
const int tens[] = { 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000 };
const int radix  = 100000000;

class Number   // Простой и очень ограниченный класс под нашу задачу
{
public:
    Number(int N = 0)
    {
        x = new int[PRECISION];
        memset(x,0,sizeof(int)*PRECISION);
        x[0] = N;
        n = 0;
    }
    void operator = (const Number& N)
    {
        memcpy(x,N.x,sizeof(int)*PRECISION);
        n = N.n;
    }
    ~Number() { delete[] x; }
private:
    int * x;   // Массив для представления чисел кусками по 10^8
    int   n;   // Последний задействованный элемент массива,
               // все выше - нули.
    // k-я цифра
    friend int  n_th_digit(const Number& N, int k);
    // N1 += N2
    friend void add(Number& N1, const Number& N2);
    // N += digit*10^{eex}
    friend void add(Number& N,  int digit, int eex);
    // N *= digit, digit < 10
    friend void mul(Number& N, int digit); // Не больше 10!!!
    // N *= 10^k
    friend void mul10(Number& N, int k);   // * 10^k
    // Вывод числа на stdout
    friend void out(const Number& N);
};


inline int n_th_digit(const Number& N, int k)
{
    // Ищем элемент массива
    int x = N.x[k/8];
    return (x/tens[k%8])%10;
}

inline void add(Number& N1, const Number& N2)
{
    int n = (N1.n > N2.n) ? N1.n : N2.n;
    int carry = 0;
    for(int i = 0; i <= n; ++i)
    {
        N1.x += N2.x + carry;
        if (N1.x >= radix)
        {
            carry = 1;
            N1.x -= radix;
        }
        else carry = 0;
    }
    if (carry)
    {
        N1.x[++n] = 1;
    }
    N1.n = n;
}

inline void add(Number& N,  int digit, int eex)
{
    int k = eex/8;
    digit *= tens[eex%8];
    N.x[k] += digit;
    if (N.x[k] >= radix)
    {
        N.x[k+1] += N.x[k]/radix;
        N.x[k] %= radix;
        ++k;
    }
    if (k > N.n) N.n = k;
}


inline void mul(Number& N, int d)
{
    int carry = 0;
    for(int i = 0; i <= N.n; ++i)
    {
        N.x *= d;
        N.x += carry;
        if (N.x >= radix)
        {
            carry = N.x/radix;
            N.x %= radix;
        }
        else carry = 0;
    }
    if (carry)
    {
        N.x[++N.n] = carry;
    }
}

inline void mul10(Number& N, int k) // * 10^k
{
    // сперва - сдвиг на ячейки. На сколько?
    int shift = k/8;
    if (shift)
    {
        memmove((void*)(N.x+shift),
                (const void*)N.x,(N.n+1)*sizeof(int));
        memset(&N.x[0],0,sizeof(int)*shift);
        N.n += shift;
    }
    if (k%8)
    {
        k %= 8;
        int carry = 0;
        for(int i = shift; i <= N.n; ++i)
        {
            int x = N.x;
            N.x = (x%tens[8-k])*tens[k] + carry;
            carry = x/tens[8-k];
        }
        if (carry)
        {
            N.x[++N.n] = carry;
        }
    }
}

void out(const Number& N)
{
    printf("%d",N.x[N.n]);
    for(int i = N.n-1; i>=0; --i)
        printf("%08d",N.x);
}


int main(int argc, const char * argv[])
{

    if (argc > 1)
    {
        LAST = atoi(argv[1]);
        if (LAST < 3) LAST = 1000;
    }

    PRECISION = LAST/4+10;  // С небольшим запасом...

    // Начальные значения
    Number x = 5, v = 6, y = 25, w = 36, g;
    int a, tenn = 1;

    // Засекаем время
    clock_t start = clock();
    for(int n = 1; n < LAST; ++n, tenn++)
    {
        a = n_th_digit(y,tenn);
        g = x;
        mul(g,2);
        add(g,a,tenn);
        mul(g,a);
        mul10(g,tenn);
        add(y,g);
        add(x,a,tenn);

        a = n_th_digit(w,tenn);
        if (a != 0) a = 10 - a;
        g = v;
        mul(g,2);
        add(g,a,tenn);
        mul(g,a);
        mul10(g,tenn);
        add(w,g);
        add(v,a,tenn);

        // Чтоб отслеживать, что не висим :)
        if (n%1000 == 0) {
            printf("N = %7d, %8d s\n",n,(clock()-start)/CLOCKS_PER_SEC);
        }
    }
    clock_t stop = clock();
    printf("%d s\n",(stop-start)/CLOCKS_PER_SEC);

    // Выводим числа
    out(x);
    printf("\n");
    printf("\n");
    out(v);
    printf("\n");
    printf("\n");
}