21.

Mathematica:
    AmicableQ[x_] := (
     d1 = Apply[Plus, Take[Divisors[x], {1, -2}]];
     If[d1 == 0, d1 = 1];
     d2 = Apply[Plus, Take[Divisors[d1], {1, -2}]];
     If[d1 == d2, Return[{-1}]];
     If[x == d2, Return[Sort[{d1, d2}]], Return[{-1}]];
    );
    numbers = {};
    For[i = 1, i < 10000, i++, {
      result = AmicableQ[i];
      If[Length[result] == 2, AppendTo[numbers, i]]
     }
    ];
    Print[Apply[Plus, numbers]];

24.

Mathematica:
    Permutations[{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}][[1000000]]

56.

Mathematica:
    sum = 0;
    For[a = 1, a < 100, a++,
     {
      For[b = 1, b < 100, b++,
       {
        tmp = Apply[Plus, ToExpression[Characters[ToString[a^b]]]],
        If[sum < tmp/b>, sum = tmp]
       }
      ]
     }
    ]
    Print[sum];

58.

Analyse:
给原题目中的加上一层，我们可以得到：

65 64 63 62 61 60 59 58 57
66 37 36 35 34 33 32 31 56
67 38 17 16 15 14 13 30 55
68 39 18  5  4  3 12 29 54
69 40 19  6  1  2 11 28 53
70 41 20  7  8  9 10 27 52
71 42 21 22 23 24 25 26 51
72 43 44 45 46 47 48 49 50
73 74 75 76 77 78 79 80 81

现设x是正副对角线上总的数字个数，不难得出
x = floor(layer / 2) * 4 + 1

观察upper right对角线，
1,  3, 13, 31, 57
2, 10, 18, 26
8,  8,  8

可以看出这是一个二阶差分
UR[x] = Ax^2+Bx+C
UR[1] = 1 =  A +   B + C
UR[2] = 3 = 4A +  2B + C
UR[3] = 13= 9A +  3B + C

解出：
A = 4
B = -10
C = 7

则：
UR[x_] := 4x^2 - 10x + 7

upper left对角线，同理可得：
1,  5, 17, 37, 65
4, 12, 20, 28
8,  8,  8
UL[x] = Ax^2+Bx+C
UL[1] = 1 =  A +   B + C
UL[2] = 5 = 4A +  2B + C
UL[3] = 17= 9A +  3B + C
A = 4
B = -8
C = 5
UL[x_] := 4x^2 - 8x + 5

down left对角线，同理可得：
1,  7, 21, 43, 73
6, 14, 22, 30
8,  8,  8
DL[x] = Ax^2+Bx+C
DL[1] = 1 =  A +   B + C
DL[2] = 7 = 4A +  2B + C
DL[3] = 21= 9A +  3B + C
A = 4
B = -6
C = 5
DL[x_] := 4x^2 - 6x + 3

down right对角线最简单，
DR[x_] := x^2

然后就是Mathematica的工作了。

Mathematica:
    UR[x_] := 4 x^2 - 10 x + 7;
    UL[x_] := 4 x^2 - 8 x + 5;
    DL[x_] := 4 x^2 - 6 x + 3;
    DR[x_] := x^2;

    PrimeQI[x_] := (Return[Length[DeleteCases[PrimeQ[x], False]]]);
    Layer[layer_] := (
        index = Floor[layer/2];
        total = index*4 + 1;
        primes = 0;
        Return[(
             PrimeQI[UR[Range[1, index + 1]]] +
             PrimeQI[UL[Range[1, index + 1]]] +
             PrimeQI[DL[Range[1, index + 1]]]
             )/total];
        );
     For[i = 26001, i < 29001, i += 2,
     {
         If[N[Layer[i]] < 0.1, {Print[i], Break[]}]
     }
    ];

63.

Mathematica:
    powers = Table[Range[1, 100]^x, {x, Range[1, 100]}];
    count = 0;
    For[power = 1, power <= 100, power++, {
      digitsOfPower =
       IntegerDigits[Take[powers[[power]], {1, -1}]],
       For[index = 1, indexdigitsOfPower], indexdigitsOfPower[[indexpower, count++]
       }
      ]
     }
    ]
    Print[count];