Вычисление среднего арифметического в списке чисел

a = [40, 43, 49]
b = sum(a) / len(a)
while b < 80:
    a.append([80])
    b = sum(a) / len(a)
print(b)
  • Элемент списка

Нужно посчитать среднее в списке до тех пор пока оно не станет равно 80. Если среднее в списке меньше 80, надо добавить в него 80, пока среднее не станет равно 80.


Ответы (2 шт):

Автор решения: Stanislav Volodarskiy

"Торжественно клянусь, что замышляю только шалость"

Интересно, что этот код после правки будет работать и рано или поздно завершится. Надо только узнать когда.

Правка

Правка состоит в удалении квадратных скобок:

a = [40, 43, 49]
b = sum(a) / len(a)
while b < 80:
    a.append(80)  # было a.append([80])
    b = sum(a) / len(a)
print(b)

Запускайте и ждите. А я пока сосчитаю сколько вам придётся ждать.

Длина списка

В списке хранятся целые числа, которые Питон суммирует точно и без переполнений. Как правильно заметили в комментариях к вопросу, среднее списка у которого почти все элементы равны восьмидесяти, а некоторые меньше, не может стать равным восьмидесяти. Математика не позволяет. Но деление sum(a) / len(a) делается в машинной арифметике ограниченной точности. То есть, результат округляется и, в конце концов, окажется равным восьмидесяти.

Код ниже отыскивает сколько надо добавить чисел чтобы программа остановилась:

def bin_search(lo, hi, pred):
    while lo + 1 < hi:
        mid = (lo + hi) // 2
        if pred(mid):
            hi = mid
        else:
            lo = mid
    return hi


a = [40, 43, 49]


def is_equal(k):
    s = sum(a) + 80 * k
    l = len(a) + k
    return s / l == 80


p = 0
while not is_equal(k := 2 ** p):
    p += 1
print(bin_search(k // 2, k, is_equal))

Требуется в районе 1.5·1016 чисел:

$ python find-k.py
15199648742375421

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

Время работы

Время работы программы квадратично зависит от количества добавленных чисел. Это из-за вызова sum(a) в цикле, который сам по себе линейный.

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

import threading


def stop():
    print(len(a))
    a.append(1000) # смешной способ остановить программу
    

t = threading.Timer(float(input()), stop)
t.start()

a = [40, 43, 49]
b = sum(a) / len(a)
while b < 80:
    a.append(80)  # было a.append([80])
    b = sum(a) / len(a)
print(b)
$ echo 1 | python benchmark.py
24768
80.03278159063383

$ echo 4 | python benchmark.py
49496
80.01640470322033

$ echo 16 | python benchmark.py
98777
80.00821903942507

Верхнее число – длина списка, нижнее – мусор, не обращайте внимания. Видно что при увеличении времени в четыре раза, длина списка увеличивается в два раза с хорошей точностью:

98777/49496 ≈ 1.996
49496/24768 ≈ 1.998

Что подтверждает теорию о квадратичной сложности. Время работы квадратичной программы будет T = Ck2. Тогда C = T/k2. Считаем для всех трёх тестов:

1/247682 ≈ 1.630·10-9
4/494962 ≈ 1.633·10-9
16/987772 ≈ 1.640·10-9

Возьмём среднее C = 1.634·10-9 и подставим k = 15199648742375421 из предыдущей главки.

T = Ck2 = 1.634·10-9·151996487423754212 ≈ 3.775·1023 секунд ≈ 6.292·1021 минут ≈ 1.048·1020 часов ≈ 4.369·1018 дней ≈ 1.197·1016 лет

Короче, это значительно (примерно в миллион раз) больше чем возраст Вселенной.

Память

В комментарии сказали что список в память не поместится. Проверим. У меня 64-битная система. На ней Питон может адресовать списки длиной до 263 - 1, значение которое выдаёт вызов sys.maxsize. Я проверил, что 15199648742375421 < 263 - 1, то есть контейнер такой длины создать можно.

Но хватит ли памяти под такой список? Список в Питоне реализуется как массив указателей на элементы. В нашем случае элементы памяти почти не расходуют, потому что все элементы списка ссылаются на одни и те же четыре значения. Питон кэширует маленькие целые, и следит за тем чтобы число 80 было в программе в единственном экземпляре. Сам массив списка требует 8k + 56 байт памяти:

@>>> sys.getsizeof([80] * 10000)
80056
@>>> sys.getsizeof([80] * 100000)
800056
@>>> sys.getsizeof([80] * 1000000)
8000056

Общая память под окончательный список помещается в память 64-битной системы: 8·15199648742375421 + 56 < 264. Причём со стократным запасом:

(8·15199648742375421 + 56) / 264 ≈ 0.0066

Запас нужен так как при наращивании списка иногда может понадобится в три раза больше памяти. Но мы помещаемся. Надо только купить память:

3(8·15199648742375421 + 56) ≈ 356241767399424.1 KiB ≈ 347892350976 MiB ≈ 339738624 GiB ≈ 331776 TiB ≈ 324 PiB

Модуль памяти SODIMM (предпочитаю работать на ноутбуке) 64GiB весит не более 10 грамм. 339738624 GiB тогда будут весить 339738624/64·10 г = 53084160 г ≈ 53084 кг ≈ 53 тонны.

53 тонны это 53 тонны. Домой я такой ноутбук принести не смогу, поставлю во дворе.

Аналогичные оценки для 32-битной версии Питона оставлю в качестве упражнения заинтересованному читателю.

"Шалость удалась"

→ Ссылка
Автор решения: mrgervant

Попробую со своей стороны изложить самую главную часть ответа на сущность вопроса:

Как из списка, в котором три числа меньше 80 и остальные равны 80, вычислить среднее арифметическое значение равное 80?

В ответе @stanislav-volodarskiy объяснение этого момента было опущено:

Но деление sum(a) / len(a) делается в машинной арифметике ограниченной точности. То есть, результат округляется и, в конце концов, окажется равным восьмидесяти.


Задача с точки зрения математики

Интуитивно ощущается, что невозможно прийти к среднему арифметическому 80, если исходные числа меньше 80, а после добавляются только числа равные 80 - всегда будет не хватать значений больше 80, чтобы "выйти в плюс".

Но как это можно доказать?

Для начала можно построить график функции, который будет отображать динамику среднего арифметического списка (y) при увеличении числа 80-ок (x):

# Исходная функция в приложении
y = sum(a) / len(a)

# Вычленяем исходные числа от частей списка только с 80
# где a_80 - это список `а` без начальных чисел неравных 80
y = sum(a_80) + (40 + 43 + 49) / len(a_80) + 3

# Нахождение суммы списка заменяем на умножение на `x` - количество 80-ок
# т.к. в списке теперь только одно и тоже повторяющееся число
y = (80 * x + (40 + 43 + 49)) / (x + 3)

Упрощаем и получаем итоговое уравнение графика:

y = (80 * x + 132) / (x + 3)

Если его построить, то выйдет такой график: введите сюда описание изображения

Код построения графика:

import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, 500, 500)
y = (80 * x + 132) / (x + 3)

plt.plot(x, y)
plt.xlabel('x')
plt.ylabel('y')
plt.title('График функции y = (80*x + 132) / (x + 3)')
plt.grid()

plt.show()

По данному графику видно, что он стремится к 80, но не может данного предела достичь. И это тоже можно доказать - через нахождение асимптоты.

Для начала - мы легко можем сказать о наличии вертикальной асимптоты, т.к. в вышеуказанном уравнении есть операция деления. Соответственно, нельзя допустить, чтобы при определенном значении x в знаменателе получался 0:

x + 3 = 0
x = -3  # уравнение вертикальной асимптоты

Но нас она не интересует - в списке 80-ок не может быть отрицательное число этих самых 80-ок: в контексте задачи минимальное значение x равно 0. Поэтому перейдем к обозначенному нахождению горизонтальной асимптоты.

Наклонная асимптота - это прямая вида y = kx + b. А горизонтальной она будет при k = 0 - проверим это для нашего уравнения:

# `k` ищется по формуле
k = lim(x -> ∞) f(x) / x

# Вместо `f(x)` подставляем наш `y`
k = ((80 * x + 132) / (x + 3)) / x

# Деление на `x` равно умножению на `1 / x`
k = (80 * x + 132) / ((x + 3) * x)

# Перемножаем полученный знаменатель
k = (80 * x + 132) / (x**2 + 3 * x)

# Если мы подставим сейчас `∞`, то получим неопределенность: ∞ / ∞
k = (80 * ∞ + 132) / (∞**2 + 3 * ∞)
# Числитель
80 * ∞ = ∞
∞ + 132 = ∞
# Знаменатель
∞ * ∞ = ∞
3 * ∞ = ∞
∞ + ∞ = ∞
# Результат
k = ∞ / ∞

# Это неопределённость, так как результат деления может быть равен любому числу
# Поэтому используем метод деления числителя и знаменателя на старшую степень переменной x
# https://ru.wikipedia.org/wiki/Раскрытие_неопределённостей
# В нашем случае это x**2 из знаменателя
k = (((80 * x)/x**2) + (132/x**2)) / ((x**2/x**2) + ((3 * x)/x**2))

# Упрощаем что можем
k = ((80/x) + (132/x**2)) / (1 + (3/x))

# Подставляем `∞` вместо `x`
k = ((∞/x) + (132/∞**2)) / (1 + (3/∞))
# Числитель
80 / ∞ = 0
∞ * ∞ = ∞
132 / ∞ = 0
0 + 0 = 0
# Знаменатель
3 / ∞ = 0
1 + 0 = 1
# Результат
k = 0 / 1
k = 0

В результате действительно получаем k = 0, что подходит для горизонтальной асимптоты. Но нужно ещё найти коэффициент b:

# `b` ищется по формуле
b = lim(x -> ∞) f(x) - (k * x)

# Но при `k = 0` это сводится просто к `f(x)` - то есть нашему `y`
b = (80 * x + 132) / (x + 3)

# Аналогично делим на старшую степень переменной - это просто `x`
b = ((80 * x / x) + 132/x) / (x/x + 3/x)

# Упрощаем что можем
b = (80 + 132/x) / (1 + 3/x)

# Подставляем `∞` вместо `x`
b = (80 + (132/∞)) / (1 + (3/∞))
# Числитель
132 / ∞ = 0
80 + 0 = 80
# Знаменатель
3 / ∞ = 0
1 + 0 = 1
# Результат
b = 80 / 1
b = 80

Таким образом, мы получаем для функции y = (80 * x + 132) / (x + 3) предел (в виде горизонтальной асимптоты):

b = 80

Иначе говоря, мы доказали, что y (среднее арифметическое значений списка) не может достигнуть значения 80 ни при каких значениях x (число 80-ок).


Задача с точки зрения информатики

Вместе с тем, в доказательстве выше мы регулярно оперировали понятием бесконечности, что мало применимо к работе реальных компьютерных систем - человечество не имеет бесконечных вычислительных мощностей / ресурсов памяти и т.д.

Следовательно, реализация математических концепций будет иметь ограничения или компромиссы. И одним из таких является стандарт IEEE 754 - стандарт, описывающий формат представления чисел с плавающей точкой (он необходим для единообразной реализации работы с данными числами разными производителями аппаратных и программных элементов компьютеров). Именно к такому формату приводятся "человеческие" числа для возможности работы с ними компьютеров, используя в качестве порядка - степени двойки (поскольку работа ведется в двоичной системе счисления).

В плане целых чисел Python отличается от многих языков программирования - для работы с целыми числами он реализует длинную арифметику, что позволяет не выставлять ограничения по числу бит для их хранения. Но для работы с типом float используется тип double из C (статья с разбором) - это указано в описании числовых типов в Python - который, в свою очередь, основан на IEEE 754. И это уже накладывает свои ограничения - если мы берем свежие 64-битные системы, то необходимо смотреть на формат binary64.

Почему мы смотрим на float, а не int? Потому что при вычислении среднего значения используется операция деления, что приводит числа к типу float.

В рамках данного формата для представления числа выделяется 64 бита, из которых:

  • 1 бит - знак, то есть это положительное или отрицательное число
  • 11 бит - порядок, то есть степень двойки, на которую умножается мантисса
  • 52 бит - мантисса, то есть число после преобразований к формату числа с плавающей точкой

Например, число 100 будет представлено как:

0100000001011001000000000000000000000000000000000000000000000000

s = 0
e = 10000000101
m = 1001000000000000000000000000000000000000000000000000

Причем изначально мантисса вычислена как 1100100 - первая единица всегда опускается (вычисления приходят к 2**0 == 1), а с конца она дополняется нулями до 52 бит.

Но может быть и обратная ситуация - мантисса вычислена слишком большой, больше 52 знаков. Соответственно, её придется округлить, чтобы уместиться в установленный формат.

Здесь обратим внимание на число, найденное алгоритмом @stanislav-volodarskiy - 15199648742375421 или сколько 80-ок потребуется в списке для того, чтобы Python выдал равенство среднего арифметического и 80.

При нём в итоге будет выполняться следующая операция:

# Числитель
a = 80 * 15199648742375421 + 132
# Знаменатель
b = 15199648742375421 + 3
# Результат
c = 1215971899390033812 / 15199648742375424

К примеру такой калькулятор вычисляет данное выражение как 79.9999999999999929. Но Python выведет 80.0 и True при сравнении с 80.

Попробуем представить указанные числа в формате binary64, а для наглядности выполним первичные преобразования сами, до интересующего нас момента:

def convert_to_binary(num):
    remainders = []
    while num > 0:
        remainders.append(num % 2)
        num //= 2
    remainders.reverse()
    return remainders
    
def get_index(value, values_list):
    if value not in values_list:
        return -1
    return values_list.index(value)
    
def integer_to_ieee754(num):
    # определяем знак числа
    sign = [0] if num >= 0 else [1]
    print(sign)
    
    # получаем список остатков деления на 2 - без 1 единицы от 2**0
    remainders = convert_to_binary(num)[1:]
    print(remainders)
    
    # первичный порядок - по степени двойки для перемножения
    raw_exponent = len(remainders)
    print(raw_exponent)
    
    # убираем незначащие нули с конца
    last_one_index = get_index(1, remainders[::-1])
    raw_mantissa = remainders
    if last_one_index > 0:
        raw_mantissa = remainders[:-last_one_index]
    print(raw_mantissa)
    print(len(raw_mantissa))

На основе данного кода мы получим начальные данные для преставления числа в формате.

Для начала посмотрим на знаменатель 15199648742375424:

integer_to_ieee754(b)
>>> sign = [0]
>>> remainders = [1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
>>> raw_exponent = 53
>>> raw_mantissa = [1, 0, 1, 1]
>>> len(raw_mantissa) = 4

С ним никаких проблем не предвидится - мантисса сводится к 4 символам, что спокойно поместится в 52 бита.

Посмотрим на числитель 1215971899390033812:

integer_to_ieee754(b)
>>> sign = [0]
>>> remainders = [0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0]
>>> raw_exponent = 60
>>> raw_mantissa = [0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1]
>>> len(raw_mantissa) = 58

Вот здесь возникают вопросы - число значений мантиссы 58, что больше 52. Для помещения в формат будет применено округление мантиссы.

Вариантов округления несколько и не уверен, что смогу корректно их описать. Поэтому обратимся к ресурсам, предлагающим конвертировать число к формату IEEE-754, как вариант - данный конвертер:

1215971899390033812 преобразуется в формат с плавающей точкой IEEE-754:
0 10000111011 0000111000000000000000000000000000000000000000000000

До округления:
1.0000110111111111111111111111111111111111111111111111 100101 ("лишние" биты)
                                                     ^данный бит определяется как наименее значащий

После округления:
1.0000111000000000000000000000000000000000000000000000

В результате округления мантисса, а следовательно и само число меняется. Если конвертировать обратно из IEEE-754, то получим большее число: 1215971899390033920

То есть, из-за округления мантиссы для возможности "поместиться" в 52 бита произошло изменение числа и фактически Python вычислял выражение 1215971899390033920 / 15199648742375424

Оно уже определяется как == 80, что остановит цикл вычислений изначального кода.


Кроме того, исходя из пределов формата binary64, начинать искать подходящие значения можно от числителя, превышающего 9007199254740992:

  • числа до и включая разряд 2*52 - помещаются, т.к. предел мантиссы 52 бита
  • 2**53 == 9007199254740992 - тоже помещается, т.к. 1 число в разряде со всеми нулями (представляется как: 0 10000110100 0000000000000000000000000000000000000000000000000000)
  • 2**53 + 1 и далее - зона с рисками искажений из-за округлений

До числа 9007199254740992 (включительно) формата хватает и ошибок не возникает, Python-приложение до этого момента будет работать согласно вышеописанной математической логике. Но далее вступает в силу ограничение стандарта и появляется шанс на получение желаемого результата из-за округлений.


Случайности не случайны

Если же посмотреть ещё пристальней, то искажение числителя 1215971899390033812 равно 108 (1215971899390033920 - 1215971899390033812 = 108)

А сколько нам не хватало в исходном списке [40, 43, 49] до того, чтобы его среднее арифметическое было равно 80?

# При 3-х элементах списка для среднего 80 требуется сумма:
80 * 3 = 240

# Имеющаяся сумма:
40 + 43 + 49 = 132

# Недостаток до среднего в 80
240 - 132 = 108

То есть, желаемая ошибка наступила в момент, когда искажение числителя в положительную сторону стало равно размеру недостатка (при учете, что знаменатель был представлен без искажений).

Таким образом, возможный размер ошибки преобразования числителя должен быть >= размера недостатка. Поскольку преобразование и округление идет в бинарном формате, то и уровни искажения будут степенями двойки. В нашем случае ближе всего к 108 будет 128 или 2**7 (меньше не берем, т.к. интересует объем, в который "поместится" наше число). Иначе говоря, меняя значения 7 бит с конца - мы имеет возможность варьироваться в пределах 128 значений, и, соответственно, потерять их при округлении.

Однако, нам потребуется значение возможной ошибки в 2 раза больше - 256. Всё из-за способов округления мантиссы - упрощенно говоря, число округляется в сторону ближайшего числа, которое без ошибок приводится к формату IEEE-754. Продемонстрируем на нашем значении:

Число 1215971899390033812 для IEEE-754 находится между:
1215971899390033664 - 0 10000111011 0000110111111111111111111111111111111111111111111111
1215971899390033920 - 0 10000111011 0000111000000000000000000000000000000000000000000000

Числа от 1215971899390033665 до 1215971899390033791 будут равны 1215971899390033664
# округляется вниз

Числа от 1215971899390033792 до 1215971899390033919 будут равны 1215971899390033920
# округляется вверх

В результате диапазон между 1215971899390033664 и 1215971899390033920 делится между ними пополам: 256 / 2 = 128 - и уже в рамках половины диапазона можно допустить ошибку в требуемое значение бит.

Из этого следует, что если осуществлять поиск по числителю, то мы можем сократить порог поиска с 2**53 (начало общей зоны с рисками искажений для мантиссы в 52 бита) до той степени, при которой в зоне искажения будет требуемый объем бит:

import math

a = [40, 43, 49]
x = 80
m_limit = 52

gap = x * len(a) - sum(a)
>>> gap = 108

def nearest_power_of_two(n):
    if n <= 0:
        return 1
    return math.ceil(math.log2(n))

numerator_power = nearest_power_of_two(gap) + 1 + m_limit
# ближайшая степень (вверх) для 108
# + 1 (удвоение из-за деления диапазона)
# + предельная точная степень мантиссы
>>> Степень с которой начнутся ошибки числителя требуемого размера = 60

Но, как было указано, при таком способе поиска ещё и должен совпасть такой знаменатель, который не исказится при представлении в формат IEEE-754. В такой ситуации среднее у нас выйдет ровно 80.

Знаменатель тоже может искажаться, но ему требуется меньшее окно искажения: math.ceil(132 / 80) = 2 (число дополнительных 80-ок, которые должны появиться из-за негативного округления знаменателя). Соответственно, и произойти это может в более ранних порядках.

Опять же, на нашем примере такое могло произойти примерно так (при допущении корректного представления числителя):

1215971899390033812 / 15199648742375422 == 80.0
                                      ^ меньше на 2

# а фактически 80.0000000000000034

Какое применение всему этому? ¯_(ツ)_/¯

→ Ссылка