求解所有的水仙花数（回归数）----大数字[40]位数量级别的运算算法

 
#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
 
#define base  100000000
#define ds 5
 
int powD[10][ds];
int sum[ds];
int num[10];
int N;
struct timeval start, end;
 
void init(void)
{
  int i, j, k;
 
  for (i = 0; i <= 9; i ++)
    for (j = 0; j < ds; j ++)
      powD[i][j] = 0;
 
  for (i = 0; i <= 9; i ++)
    num[i] = 0;
 
  for (i = 1; i <= 9; i ++)
    powD[i][0] = i;
}
 
void nextPower(void)
{
   int i, j, t;
   for (i = 2; i <= 9; i ++)
   {
       t = 0;
       for (j = 0; j < ds; j ++)
       {
           t = powD[i][j] * i + t;
           if (t >= base)
           {
               powD[i][j] = t % base;
               t /= base;
           }
           else
           {
               powD[i][j] = t;
               t = 0;
           }
       }
   }
}
 
int circle(int n, int b)
{
  int i, j, k;
  double tm;
 
  if (n == 0)
  {
    if (sort())
    {
      k = ds - 1;
      while ((sum[k] == 0) && (k >= 0))
        k --;
      gettimeofday(&end, NULL);
      tm = (end.tv_sec - start.tv_sec) * 1000000.0 + (end.tv_usec - start.tv_usec);
      tm /= 1000000.0;
      printf("%4.4f  ", tm);
      printf("%2d: ", N);
      printf("%d", sum[k]);
      if (k > 0)
        for (i = k - 1; i >= 0; i --)
          printf("%08d", sum[i]);
      printf("\n");
    }
  }
  else
     for (i = 0; i <= b; i ++)
     {
       num[i] ++;
           circle(n-1, i);
       num[i] --;
     }
}
 
 
int sort(void)
{
  int i, j, t;
  int temp[ds], d[10];
 
  for (i = 0; i < ds; i ++)
    temp[i] = sum[i] = 0;
 
  for (i = 0; i <= 9; i ++)
    d[i] = 0;
 
  for (i = 1; i <= 9; i ++)
    if (num[i])
    {
      t = 0;
      for (j = 0; j < ds; j ++)
      {
        t = powD[i][j] * num[i] + t;
        if (t >= base)
        {
            temp[j] = t % base;
            t /= base;
        }
        else
        {
            temp[j] = t;
            t = 0;
        }
      }
 
      t = 0;
      for (j = 0; j < ds; j++)
      {
        t = temp[j] + sum[j] + t;
        if (t >= base)
        {
            sum[j] = t - 100000000;
            t = 1;
        }
        else
        {
            sum[j] = t;
            t = 0;
        }
      }
    }
 
    for (i = ds - 1; i >= 0; i --)
    {
        t = sum[i];
        if (t > 0)
          while (t > 0)
          {
              d[t % 10] ++;
              t /= 10;
          }
    }
 
    for (i = 1; i <= 9; i ++)
      if (num[i] != d[i])  return 0;
    return 1;
}
 
int main(void)
{
  int i, j;
 
  gettimeofday(&start, NULL);
  init();
  nextPower(); //是平方
  for (i = 2; i <= 40; i ++)
  {
    N = i;
    printf("Search %d...\n", i);
    circle(N, 9);
    nextPower();
  }
  return 0;
}