Метод, подходящий для n до тысячи (на Python секунд 15).

Базируется на признаке делимости на 7 на основе знакопеременной суммы триад (здесь). Делим число на грани по 3 цифры, складываем их, если результат делится на 7, то и исходное число тоже.

77744744: 77-744+744 = 77 делится на 7

Для младшей триады всего два варианта, они заданы сразу, перебираем со второй триады. Получаем счётчики количества возможных сумм для всех комбинаций, выбираем те суммы, что делятся на 7.

from collections import Counter

def num74(n):
    if n < 5: return 0
    triads = [[444,447,474,477,744,747,774,777],[4,7],[44,47,74,77]]
    sums = Counter({444:1, 744:1})
    mul = 1
    for k in range(n//3-1):
        mul = - mul
        temp = Counter()
        for v in triads[0]:
            for ky in sums:
                temp[ky + v*mul] += sums[ky]
        sums = temp
    if n%3:
        mul = - mul
        temp = Counter()
        for v in triads[n%3]:
            for ky in sums:
                temp[ky + v*mul] += sums[ky]
        sums = temp

    res = sum([sums[ky] for ky in sums if (ky%7) == 0])
    return(res)

На С++:

int main()
{
    for(int n = 1; n <= 50; n++)
        cout << setw(2) << n << "  "
            << setw(20) << ((n < 2) ? 0 : (unsigned long long)((1llu<<(n-2))/7))
            << endl;
}

На Питоне, как я понимаю,

for n in range(1,1001):
    print(n,"  ",(1<<(n-2)) // 7)

https://oeis.org/A033138

Я бы использовал такой алгоритм:

parts = [[444, 447, 474, 477, 744, 747, 774, 777], [4], [44, 47, 74, 77]]
ends = [[444, 744], [4], [44]]


def calculate(count, size, pos, value):
    res = 0

    if pos < count - 1:
        for part in parts[size]:
            res += calculate(count, 0, pos + 1, value + (part if pos % 2 else -part))
    else:
        for part in ends[size]:
            if (value + (part if pos % 2 else -part)) % 7 == 0:
                return 1

    return res


n = 10
print(calculate(n // 3 + (n % 3 != 0), n % 3, 0, 0))

1. основные принципы - перебор идет по тройкам цифр
2. последняя тройка выбирается по критерию делимости на 4 - (2 последних цифры должны делиться на 4)
3. само число проверяется на делимость на 7 по критерию что сумма четных и разность нечетных троек (тысяч) числа должна делиться на 7

У меня не такое красивое решение, как угаданное @Harry, зато

• оно имеет сложность O(log n), где n - количество цифр в числе,
• оно объяснимо, доказуемо и обобщаемо.

Сначала заметим, что все числа должны иметь вид a = <b44>, то есть заканчиваться на две четверки. В противном случае a не будет делиться на 4.

Найдем остаток a % 7: (b*100 + 44) % 7 = (b*100)%7 + 44%7 = (b%7)*(100%7) + 2 = (b%7)*2 + 2

Итак, a делится на 7 и на 4 тогда и только тогда, когда a заканчивается на 44, а число b, полученное из a отбрасыванием двух последних цифр, должно делиться на 7 с остатком 6. Почему 6? Потому, что 6*2 +2 = 14 == 0 (mod 7)

Чтобы найти, сколько таких чисел длиной m = n-2 нужно немного обобщить задачу.

Пусть мы знаем гистограмму остатков для чисел, составленных из 4 и 7, длиной m: r_m = [ r0_m, r1_m, ..., r6_m], где r0_m - количество чисел длиной m, делящихся на 7 с остатком 0, r1_m - c остатоком 1 и т.д.

Как связаны r_m и r_{m+1}?

Пусть b - некоторое число длиной m. Припишем к нему справа цифры 4 и 7. Получим два числа b1 = <b4>, b2 = <b7>. Чему равны их остатки от деления на 7?

b1 % 7 = ((10*b) + 4)%7 = (b%7)*(10%7) + 4 = 3*(b%7) + 4

b2 % 7 = ((10*b) + 7)%7 = (b%7)*(10%7) + 0 = 3*(b%7)

То есть каждый остаток r порождает два остатка 3r и 3r+4. Выпишем их в явном виде. Слева остаток для числа длиной m, справа два остатка для чисел длиной m+1

0 -> 0, 4
1 -> 3, 0
2 -> 6, 3
3 -> 2, 6
4 -> 5, 2
5 -> 1, 5
6 -> 4, 1

В матричном виде это выглядит так:

Для m=0 существует только одно число, 0. Поэтому начальный вектор остатков r_0 = [1,0,0,0,0,0,0]

r_m = (A^m)@r_0

Здесь A - матрица перехода от r_m к r_{m+1}

Для ответа на исходный вопрос нужно вычислить r_{n-2} и взять у него шестой компонент (нумерация с 0), соответствующий количеству чисел длиной n-2, дающих при делении на 7 остаток 6

import numpy as np

A = np.array([[1,1,0,0,0,0,0],
              [0,0,0,0,0,1,1],
              [0,0,0,1,1,0,0],
              [0,1,1,0,0,0,0],
              [1,0,0,0,0,0,1],
              [0,0,0,0,1,1,0],
              [0,0,1,1,0,0,0]
             ], dtype=np.object)

r0 = np.array([1,0,0,0,0,0,0])

for i in range(3, 101):
    B = np.linalg.matrix_power(A, i-2)
    r = B@r0
    print(i, r[6])

В конструкторе A явно задаётся тип object, чтобы элементы были пайтоновскими целыми произвольной длины. Если тип не указать, то numpy соптимизирует и будет хранить целые как 32-х битные.

Для возведения в степень используется numpy.linalg.matrix_power, так как встроенная функция pow возводит в степень поэлементно.

Результат

3 0
4 0
5 1
6 2
7 4
8 9
9 18
10 36
11 73
12 146
13 292
14 585
15 1170
16 2340
17 4681
18 9362
19 18724
20 37449
21 74898
22 149796
23 299593
24 599186
25 1198372
26 2396745
27 4793490
28 9586980
29 19173961
30 38347922
31 76695844
32 153391689
33 306783378
34 613566756
35 1227133513
36 2454267026
37 4908534052
38 9817068105
39 19634136210
40 39268272420
41 78536544841
42 157073089682
43 314146179364
44 628292358729
45 1256584717458
46 2513169434916
47 5026338869833
48 10052677739666
49 20105355479332
50 40210710958665
51 80421421917330
52 160842843834660
53 321685687669321
54 643371375338642
55 1286742750677284
56 2573485501354569
57 5146971002709138
58 10293942005418276
59 20587884010836553
60 41175768021673106
61 82351536043346212
62 164703072086692425
63 329406144173384850
64 658812288346769700
65 1317624576693539401
66 2635249153387078802
67 5270498306774157604
68 10540996613548315209
69 21081993227096630418
70 42163986454193260836
71 84327972908386521673
72 168655945816773043346
73 337311891633546086692
74 674623783267092173385
75 1349247566534184346770
76 2698495133068368693540
77 5396990266136737387081
78 10793980532273474774162
79 21587961064546949548324
80 43175922129093899096649
81 86351844258187798193298
82 172703688516375596386596
83 345407377032751192773193
84 690814754065502385546386
85 1381629508131004771092772
86 2763259016262009542185545
87 5526518032524019084371090
88 11053036065048038168742180
89 22106072130096076337484361
90 44212144260192152674968722
91 88424288520384305349937444
92 176848577040768610699874889
93 353697154081537221399749778
94 707394308163074442799499556
95 1414788616326148885598999113
96 2829577232652297771197998226
97 5659154465304595542395996452
98 11318308930609191084791992905
99 22636617861218382169583985810
100 45273235722436764339167971620

Так как возведение в степень m использует алгоритм со сложностью O(log m), то сложность для n равна O(log(n-2)), то есть O(log n)

Алгоритм можно использовать и для других условий. Например, для чисел, составленных из цифр 5 и 9, которые делятся на 5 и 9. Будет другая матрица и другое условие для остатоков, но схема решения останется той же.

Можно было даже обойтись без оптимизации с двумя четверками . В этом случае матрица размером 28 на 28 .

A28 = np.zeros((28,28), dtype=np.object)
# подготовка матрицы переходов
for r in range(28):
    i1 = (r*10+4)%28
    i2 = (r*10+7)%28
    A28[i1, r] = 1
    A28[i2, r] = 1

# начальный вектор остатков
r028 = np.zeros(28, dtype=np.object)
r028[0] = 1

for i in range(1, 101):
    B = np.linalg.matrix_power(A28, i)
    r = B@r028
    print(i, r[0])

результат тот же.

Решение от @Harry

Вывести точную формулу я не возьмусь, но обосновать 2^(n-2)//7 легко. Плюс минус единица )

Подходящих чисел длины n всего 2^(n-2), так как последние две цифры фиксированы. Остатки от деления на 7 равномерно распеределены, отличаясь не более чем на 1. Поэтому количество чисел с фиксированным остатоком - в нашем случае остаток 6, - равно 2^(n-2) / 7 плюс минус 1.

Как строго доказать тезис о равномерном распределении остатков, я не представляю. Из общих соображений за это говорит структура матрицы A - в каждой строке ровно две единицы, строки перемешаны, поэтому матрица A должна хорошо выравнивать гистограмму остатков, перемешивая максимумы с минимумами.

Самое удивительное в том, что формула работает не только для больших чисел, но и для маленьких, где флуктуации плюс минус 1 вполне соизмеримы с значениями функции