Главная страница

Алгоритм увеличения точности нейронных сетей и его приложения


Скачать 3.63 Mb.
НазваниеАлгоритм увеличения точности нейронных сетей и его приложения
Дата24.06.2022
Размер3.63 Mb.
Формат файлаpdf
Имя файлаdiss_avrutskiy.pdf
ТипДиссертация
#613262
страница7 из 13
1   2   3   4   5   6   7   8   9   10   ...   13
𝜕𝑉
𝜕𝑥
𝑖
,
𝜕
2
𝑉
𝜕𝑥
𝑖
𝜕𝑥
𝑗
,
𝜕
3
𝑉
𝜕𝑥
𝑖
𝜕𝑥
𝑗
𝜕𝑥
𝑘
, 𝑖, 𝑗, 𝑘 ∈ [1, 𝑁 ] .

67
Минимизируемая величина будет иметь вид
𝑒
все
3
= 𝑉
2
+
∑︁
𝑖
(︂ 𝜕𝑉
𝜕𝑥
𝑖
)︂
2
+
∑︁
комб.𝑖,𝑗,𝑘
(︂
𝜕
2
𝑉
𝜕𝑥
𝑖
𝜕𝑥
𝑗
)︂
2
+
∑︁
комб.𝑖,𝑗,𝑘
(︂
𝜕
3
𝑉
𝜕𝑥
𝑖
𝜕𝑥
𝑗
𝜕𝑥
𝑘
)︂
2
,
где в силу симметричности производных, суммирование идет только по различным комби- нациям индексов. Так как нас интересует переобучение, процедуру исключения проводить не требуется.
3.2.2. Производные по осям. Простейшим способом уменьшения количества допол- нительных производных с
𝑂
(︀𝑁
3
)︀
до
𝑂 (𝑁 )
является отбрасывание всех членов кроме про- изводных по осям
𝜕𝑉
𝜕𝑥
𝑖
,
𝜕
2
𝑉
𝜕𝑥
2
𝑖
,
𝜕
2
𝑉
𝜕𝑥
3
𝑖
, 𝑖 ∈ [1, 𝑁 ] ,
𝑒
оси
3
= 𝑉
2
+
∑︁
𝑖
(︂ 𝜕𝑉
𝜕𝑥
𝑖
)︂
2
+
∑︁
𝑖
(︂ 𝜕
2
𝑉
𝜕𝑥
2
𝑖
)︂
2
+
∑︁
𝑖
(︂ 𝜕
3
𝑉
𝜕𝑥
3
𝑖
)︂
2 3.2.3. Случайные направления. Дифференцирование невязки по всем осям фактиче- ски является дифференцированием по
𝑁
ортогональным направлениям, которые одинаковы для всех точек сетки. В силу локальности алгоритма вычисления производных, эти направ- ления можно выбирать независимо в каждой точке. Таким образом, можно рассмотреть ошибку, содержащую производные по некоторому количеству
𝐿
направлений
𝑗
𝑙
𝜕𝑉
𝜕⃗𝑗
𝑙
,
𝜕
2
𝑉
𝜕⃗
𝑗
𝑙
2
,
𝜕
3
𝑉
𝜕⃗𝑗
3
𝑙
, 𝑙 ∈ [1, 𝐿] ,
причём
𝑗
𝑙
= 𝑗
𝑙
(𝑎) ,
то есть направление зависит от индекса паттерна. Так как число неза- висимых первых производных не превосходит
𝑁
, положим
1 ≤ 𝐿 ≤ 𝑁
. Локальная ошибка может быть записана как
𝑒
нап
3
= 𝑉
2
+
𝐿
∑︁
𝑙=1
(︂ 𝜕𝑉
𝜕⃗𝑗
𝑙
)︂
2
+
𝐿
∑︁
𝑙=1
(︃
𝜕
2
𝑉
𝜕⃗𝑗
2
𝑙
)︃
2
+
𝐿
∑︁
𝑙=1
(︃
𝜕
3
𝑉
𝜕⃗𝑗
3
𝑙
)︃
2
Заметим, что при таком подходе, особенно в случае малых размерностей, может образовы- ваться больше производных, чем при осевом дифференцировании. Например, если
𝑉
зави- сит от
𝑣
и
𝜕𝑣/𝜕𝑥
, то её производная
𝜕𝑉 /𝜕𝑥
будет зависеть от трех независимых величин
𝑣
,
𝜕𝑣/𝜕𝑥
,
𝜕
2
𝑣/𝜕𝑥
2
, тогда как
𝜕𝑉 /𝜕⃗𝑗
будет зависеть от четырёх:
𝑣
,
𝜕𝑣/𝜕𝑥
,
𝜕𝑣/𝜕⃗𝑗
,
𝜕
2
𝑣/𝜕𝑥𝜕⃗𝑗.
Для использования случайных направлений требуется модификация начальных условий ал- горитма прямого распространения. Она изложена в Пункте 4.1.1.
3.2.4. Тесты. Определим, насколько сильно будут отличаться результаты решения кра- евой задачи {3.2.1, 3.2.2, 3.2.3} для
𝑁 =
2, 3 и 4 при использовании трех ошибок, предло- женных в предыдущих пунктах. Для каждого случая выполним 2000 итераций минимизации

68
𝑁
2 3
4
𝑛
𝜖
𝑛
𝜖
𝑛
𝜖
𝑒
все
3 21
(1.3 ± 0.2) 10
−4 63
(1.5 ± 0.3) 10
−4 126
(2.3 ± 0.2) 10
−4
𝑒
оси
3 19
(5.5 ± 1.5) 10
−4 40
(5.5 ± 0.5) 10
−4 69
(6.6 ± 1) 10
−4
𝑒
нап 3
𝐿 = 4




117
(2.2 ± 0.4) 10
−4
𝐿 = 3


70
(1.6 ± 0.2) 10
−4 90
(3.1 ± 0.3) 10
−4
𝐿 = 2 35
(1.3 ± 0.3) 10
−4 49
(2.6 ± 0.5) 10
−4 63
(3.2 ± 0.7) 10
−4
𝐿 = 1 20
(4.2 ± 1) 10
−4 28
(4.6 ± 1) 10
−4 36
(4.3 ± 1.2) 10
−4
Таблица 1. Среднеквадратичная невязка при минимизации разных наборов дополнительных производных,
𝑛
– общее количество вычисляемых производ- ных.
RProp, и повторим процедуру решения 10 раз. Среднеквадратичные невязки и их максималь- ные/минимальные отклонения представлены в Таблице 1. Разумеется, для анализа переобу- чения нужно сравнить
𝜖
тест с
𝜖
трен
, однако, абсолютно во всех случаях эти две величины совпадали до первого знака после запятой, и поэтому не были включены в таблицу. Кро- ме того, совпадали и максимальные значения
|𝑉 |
на этих множествах. Это означает, что глобальный экстремум
𝑉
всегда оказывался расположен очень близко к одной из трениро- вочных точек. Очевидно, что условием на близкий экстремум является малость производных
𝑉
, но не совсем ясно, почему тренировочные точки «притягивают» не только минимумы, но и максимумы. Можно заметить, что аналогичная ситуация наблюдалась в Параграфе 2.4.
Итак, даже при вычеркивании всех смешанных производных, вместо появления пере- обучения, мы наблюдаем лишь равномерное увеличение
𝑉
на всём множестве решения. При минимизации
𝜕𝑉 /𝜕⃗𝑗
только для одного случайного направления в 4-мерном пространстве,
коэффициент переобучения лишь несущественно превосходит единицу,
𝑟 = 1.02
. Отличие между наборами производных лишь в том, что чем больше членов минимизируется, тем сред- няя невязка меньше. Заметим, что из-за отсутствия «перекрытия» случайных производных и членов, встречающихся в уравнении, использование двух случайных направлений в двумер- ном случае оказывается в 1.5 раза медленнее, чем обучение всех возможных производных.
Однако, при увеличении размерности до 4, два случайных направления обеспечивают вдвое меньшую нагрузку без ощутимых потерь в качестве обучения. Исходя из представленных результатов, было решено использовать два случайных направления для всех многомерных уравнений. Обозначим их
⃗𝑗
и
⃗𝑘
Вопрос о влиянии смешанных производных можно поставить и по-другому – чему будут равны их значения, если они не включены в ошибку. Так как производные по слу- чайным направлениям иногда могут перекрываться со смешанными производными, сравним случай дифференцирования
𝑉
только по осям со случаем минимизации всех возможных производных. Результаты представлены в Таблице 2.
Как можно видеть, даже если смешанные производные не минимизируются, их зна- чения вполне сравнимы с производными, которые были включены в ошибку – отношение

69
𝑉
𝜕𝑉 /𝜕𝑥
𝑖
𝜕
2
𝑉 /𝜕𝑥
2
𝑖
𝜕
2
𝑉 /𝜕𝑥
𝑖
𝜕𝑥
𝑗
𝜕
3
𝑉 /𝜕𝑥
3
𝑖
𝜕
3
𝑉 /𝜕𝑥
𝑖
𝜕𝑥
𝑗
𝜕𝑥
𝑘
𝑒
все
3 0.00022 0.00046 0.0011 0.0018 0.0049 0.0069
𝑒
оси
3 0.00066 0.0019 0.0026 0.0092 0.0075 0.023
Таблица 2. Среднеквадратичные отклонения различных производных невяз- ки при минимизации всех, либо только осевых производных.
между верхней и нижней величиной примерно одинаково для всех столбцов. Эти результа- ты схожи с эффектом отсутствия переобучения, обсуждаемого в Главе 2 – при исключении точек из тренировочного множества ошибка равномерно растёт для всех точек, а не только для исключённых.
Можно заметить, что довольно много смешанных производных точного решения
𝑢
a равны нулю. Это означает, что в окрестности решения смешанные производные
𝑉
также ма- лы. Такое обстоятельство могло улучшить результат решения при дифференцировании толь- ко по осям, так как фактически выбрасывались члены, близкие к нулю. Это предположение было опровергнуто путём добавления к
𝑢
4D
a слагаемого
(︀1 − 𝑥
2 1
− 𝑥
2 2
− 𝑥
2 3
− 𝑥
2 4
)︀ 𝑥
1
𝑥
2
𝑥
3
𝑥
4
,
с которым все смешанные члены становятся порядка единицы, тем не менее использование всех производных не получает дополнительных преимуществ перед осевым дифференциро- ванием.
3.2.5. Перенормировка. Несмотря на то, что точные решения
𝑢
a
(а также
𝑤
a
) име- ют все производные порядка единицы, в процессе обучения старшие производные
𝑣
могут оказаться разного порядка. Избежать замедления обучения в таком случае можно либо с помощью универсальной нормировки, описанной в Главе 4, либо с помощью специальной процедуры, основанной на произвольности длин
⃗𝑗
и
⃗𝑘
. Суть её состоит в следующем. В
каждой точке и для каждого направления определим вспомогательную величину

как максимум из модулей всех производных
𝑣
, в которых есть это направление:
ℳ = max
(︂⃒



𝜕𝑣
𝜕⃗𝑗




,




𝜕
2
𝑣
𝜕𝑥𝜕⃗𝑗




,




𝜕
2
𝑣
𝜕⃗𝑗
2




,




𝜕
3
𝑣
𝜕𝑥𝜕⃗𝑗
2




, . . .
)︂
Если
ℳ > 4
, то вектор
⃗𝑗
делится на
𝑦


где
𝑦
есть порядок дифференцирования мак- симальной производной по
⃗𝑗
. Аналогичная процедура выполняется для
⃗𝑘
. Перенормировка направлений будет происходить каждые
𝑇
пн итераций.
3.2.6. Процедура обучения. В Главе 2 с помощью обучения с исключением решения были получены со средней относительной точностью
∼ 10
−6
, то есть порядка машинной. Для многомерных случаев было решено воспроизвести этот результат и получить максимальное отклонение от точного решения не больше чем
1 · 10
−5
. Сперва были проведены численные эксперименты для двумерного случая и построена процедура обучения с исключением с использованием локальных ошибок
𝑒
0
= 𝑉
2
,

70
Линейное уравнение
Δ𝑢 = 𝑔
Нелинейное уравнение
Δ𝑢 + 𝑢
2
= ℎ
ошибка
𝛿𝑊 (0)
𝑇
𝑇
пн ошибка
𝛿𝑊 (0)
𝑇
𝑇
пн
1
𝐸
4 2 · 10
−4 2000 200 1
𝐸
4 2 · 10
−4 3000 15–100 2
𝐸
3 2 · 10
−5 2000 200 2
𝐸
3 2 · 10
−5 2000 200 3
𝐸
2 2 · 10
−5 2000 200 3
𝐸
2 2 · 10
−5 2000 200
Таблица 3. Параметры обучения для решения краевых задач.
𝑒
1
= 𝑒
0
+
(︂ 𝜕𝑉
𝜕⃗𝑗
)︂
2
+
(︂ 𝜕𝑉
𝜕⃗𝑘
)︂
2
,
𝑒
2
= 𝑒
1
+
(︂ 𝜕
2
𝑉
𝜕⃗𝑗
2
)︂
2
+
(︂ 𝜕
2
𝑉
𝜕⃗𝑘
2
)︂
2
,
𝑒
3
= 𝑒
2
+
(︂ 𝜕
3
𝑉
𝜕⃗𝑗
3
)︂
2
+
(︂ 𝜕
3
𝑉
𝜕⃗𝑘
3
)︂
2
,
𝑒
4
= 𝑒
3
+
(︂ 𝜕
4
𝑉
𝜕⃗𝑗
4
)︂
2
+
(︂ 𝜕
4
𝑉
𝜕⃗𝑘
4
)︂
2
Напомним, что
⃗𝑗
и
⃗𝑘
это два случайных направления, меняющихся от точки к точке. Гло- бальная ошибка последовательно понижается с четвёртого порядка до второго (здесь
𝑀

число точек сетки):
𝐸
4
=
𝑀
∑︁
𝑎=1
𝑒
4
(𝑎) , 𝐸
3
=
𝑀
∑︁
𝑎=1
𝑒
2
(𝑎) , 𝐸
2
=
𝑀
∑︁
𝑎=1
𝑒
2
(𝑎) ,
так как при анализе точности было установлено, что последние две фазы минимизации
𝐸
1
и
𝐸
0
слабо влияют на конечный результат. Все параметры решения указаны в Таблице 3
(нелинейное уравнение рассмотрено в Параграфе 3.4). Далее оказалось, что все методики обучения, полученные для двумерного случая приводят к практически идентичным резуль- татам и для размерностей 3, 4 и 5. Перейдём непосредственно к результатам численного моделирования.
3.3. Линейные уравнения
Были решены краевые задачи для линейного уравнения Пуассона внутри 2, 3, 4 и
5–мерных единичных шаров:
Δ𝑢 = 𝑔,
Γ :
𝑁
∑︁
𝑖=1
𝑥
𝑖
≤ 1,
𝑢|
𝜕Γ
= 0.

71
x
1
x
2
Рис. 3.3.1. Сетка для решения двумерного уравнения Пуассона. Тот же шаг использовался и для остальных размерностей.
Для каждой размерности было подобрано точное решение, все производные которого имеют порядок единицы, а источник
𝑔
был вычислен как его лапласиан:
𝑢
2D
a
=
10 17
(︀1 − 𝑥
2 1
− 𝑥
2 2
)︀ (︀𝑥
1
+ sin 𝑥
2
+ 𝑥
2 1
+ 𝑥
2
cos 𝑥
1
)︀ ,
(3.3.1)
𝑢
3D
a
=
3 5
(︀1 − 𝑥
2 1
− 𝑥
2 2
− 𝑥
2 3
)︀ (︀𝑥
1
+ sin 𝑥
2
+ 𝑥
2 3
+ 𝑥
2
cos 𝑥
1
)︀ ,
(3.3.2)
𝑢
4D
a
=
7 9
(︀1 − 𝑥
2 1
− 𝑥
2 2
− 𝑥
2 3
− 𝑥
2 4
)︀ (︀𝑥
1
+ sin 𝑥
2
+ 𝑥
2 3
+ 𝑥
4
cos 𝑥
4
)︀ ,
(3.3.3)
𝑢
5D
a
=
7 9
(︀1 − 𝑥
2 1
− 𝑥
2 2
− 𝑥
2 3
− 𝑥
2 4
− 𝑥
2 5
)︀ (︀𝑥
1
+ sin 𝑥
2
+ 𝑥
2 3
+ 𝑥
4
cos 𝑥
5
)︀ .
(3.3.4)
Для решения с помощью нейронной сети производилась замена
𝑢 = 𝑤 ·
(︃
1 −
𝑁
∑︁
𝑖=1
𝑥
2
𝑖
)︃
,
(3.3.5)
и уравнение приобретало вид
(︃
1 −
𝑁
∑︁
𝑖=1
𝑥
2
𝑖
)︃
Δ𝑤 − 4
𝑁
∑︁
𝑖=1
𝑥
𝑖
𝜕𝑤
𝜕𝑥
𝑖
− 6𝑤 − 𝑔 = 0.
(3.3.6)
Сетки состояли из из двух частей: множества точек на поверхности с угловым расстоянием
𝜙 ∼ 𝜋/6
и декартовой сетки внутри области с шагом
𝜆 = 1/3
, случайно повёрнутой и сдвинутой относительно начала координат. Пример для
𝑁 = 2
представлен на Рисунке
3.3.1.

72
𝑁
𝑢
a h
𝜆
𝜙
𝜗
median
𝜗
max
𝜖
𝑀
t, сек
Линейное уравнение
Δ𝑢 = 𝑔
2 3.3.1 96 1/3
𝜋/6 6.8 · 10
−7 3.5 · 10
−6 4 · 10
−5 40 114 3 3.3.2 96 1/3
𝜋/6 1.4 · 10
−6 6.5 · 10
−6 1.3 · 10
−4 159 226 4 3.3.3 148 1/3
𝜋/6 9.7 · 10
−7 6.5 · 10
−6 1.7 · 10
−4 512 560 5 3.3.4 160 1/3
𝜋/6 1.2 · 10
−6 8.7 · 10
−6 2.1 · 10
−4 1536 1280
Нелинейное уравнение
Δ𝑢 + 𝑢
2
= ℎ
2 3.3.1 96 1/4
𝜋/8 2.4 · 10
−6 8 · 10
−6 1 · 10
−4 64 180 3 3.3.2 96 1/4
𝜋/8 1.9 · 10
−6 8.2 · 10
−6 1.8 · 10
−4 343 380 4 3.3.3 148 1/4
𝜋/8 1.2 · 10
−6 8.6 · 10
−6 2.5 · 10
−4 1536 1160 5 3.3.4 160 1/4
𝜋/8 1.5 · 10
−6 8.8 · 10
−6 3.3 · 10
−4 6528 5000
Таблица 4. Результаты решения линейных и нелинейных краевых задач.
Обучение было выполнено с помощью программного комплекса «cpde-2», описанного в Главе 4. Промежуточные операции чтения/записи не использовались. Вычисления выпол- нены на облачной машине Google Cloud с GPU NVIDIA Tesla P100 в режиме 32 битной точности. Все решения для размерностей
𝑁 > 3
были получены с первой попытки. Нейросе- тевая аппроксимация для
𝑤
далее подставлялась в выражение (3.3.5), которое сравнивалось с точным решением
𝑢
a на достаточно плотном множестве равномерно распределённых слу- чайных точек (например, для двумерного случая использовалось
4 · 10 3
а для пятимерного

5 · 10 6
точек). Были посчитано медианное
𝜗
median и максимальное
𝜗
max отклонение от
𝑢
a
Результаты представлены в Таблице 4. Нейронные сети имели архитектуру вида dim
*
, h, h, h, h, h, h, 1
*
,
значение h приведено в Таблице 4. Для того, чтобы
𝜗
max оставалось строго меньше
10
−5
,
число нейронов в скрытых слоях было увеличено для четырёх и пятимерного случая, хо- тя даже при h=96 величина
𝜗
max не превышала
2 · 10
−5
. Несмотря на то, что суммарное количество арифметических операций увеличивается с каждым дополнительным измерени- ем приблизительно в 4 раза, время решения лишь удваивается. Это происходит из-за роста эффективности вычислений на GPU по мере увеличения размера перемножаемых матриц
1
При повышении размерности с 5 до 6 рост уже не будет наблюдаться.
3.4. Нелинейные уравнения
Чтобы показать, что линейность уравнения не является существенным фактором, бы- ли рассмотрены четыре краевые задачи для нелинейного уравнения Пуассона. Так как про- цесс решения схож с линейным случаем, в данном параграфе мы остановимся только на их
1
этот эффект рассмотрен в Пункте 4.4.2

73
отличиях. Итак, имеется нелинейная задача:
Δ𝑢 + 𝑢
2
= ℎ,
Γ :
𝑁
∑︁
𝑖=1
𝑥
2
𝑖
≤ 1,
𝑢|
𝜕Γ
= 0.
выражения для
𝑢
a совпадают с линейными,

определялся аналогично:
ℎ = Δ𝑢
a
+ 𝑢
2
a
. После замены
𝑢 = 𝑤 ·
(︃
1 −
𝑁
∑︁
𝑖
𝑥
2
𝑖
)︃
уравнение приобретает вид
(︃
1 −
𝑁
∑︁
𝑖
𝑥
2
𝑖
)︃
Δ𝑤 − 4
𝑁
∑︁
𝑖=1
𝑥
𝑖
𝜕𝑤
𝜕𝑥
𝑖
− 6𝑤 + 𝑤
2
(︃
1 −
𝑁
∑︁
𝑖=1
𝑥
2
𝑖
)︃
2
− ℎ = 0.
Сетки строятся по такому же принципу, как и в линейном случае. Несмотря на то, что точ- ные решения для нелинейных и линейных уравнений совпадают, шаг
𝜆 = 1/3
не давал желаемую точность
𝜗
max
< 10
−5
. Действительно, теперь невязка уравнения содержит вдвое больше слагаемых и меняется более сложным образом. Тесты в двумерном случае, реше- ние для которого длится порядка трех минут, показали, что для достижения
𝜗
max
< 10
−5
достаточно положить
𝜆 = 1/4
,
𝜙 = 𝜋/8
. Также потребовалось увеличить число итераций минимизации
𝐸
4
до 3000 шагов и дополнительно разбить эту фазу обучения на 5 интервалов,
содержащих 160, 340, 500, 1000 и 1000 шагов для которых
𝑇
пн
= 15
, 50, 50, 50 и 100 соот- ветственно. Инкременты весов
𝛿𝑊
были сброшены до
𝛿𝑊 (0)
в конце каждого интервала кроме последнего. Параметры решения собраны в Таблице 3, результаты в Таблице 4. Нели- нейная задача демонстрирует ещё большую масштабируемость – параметры обучения и шаг сетки, подобранные для двумерного случая, позволяют получить практически идентичный результат даже для
𝑁 = 5
, полуторачасовое решение которого сильно бы затрудняло подбор параметров обучения.

74 3.5. Сравнение с конечными разностями
В настоящем параграфе будет рассмотрено только линейное уравнение Пуассона. На- шей задачей является оценка времени решения {3.2.1, 3.2.2, 3.2.3} с помощью
2𝑁 + 1
точеч- ного центрального шаблона второго порядка. Так, в двумерном случае мы имеем
Δ𝑢 (𝑥
1
, 𝑥
2
) =
1
𝜆
2
(𝑢 (𝑥
1
+ 𝜆, 𝑥
2
) + 𝑢 (𝑥
1
− 𝜆, 𝑥
2
) − 2𝑢 (𝑥
1
, 𝑥
2
)) +
+
1
𝜆
2
(𝑢 (𝑥
1
, 𝑥
2
+ 𝜆) + 𝑢 (𝑥
1
, 𝑥
2
− 𝜆) − 2𝑢 (𝑥
1
, 𝑥
2
)) +
+
𝜆
2 12
(︂ 𝜕
4
𝑢
𝜕𝑥
4 1
+
𝜕
4
𝑢
𝜕𝑥
4 2
)︂
+ 𝑂
(︀𝜆
4
)︀ .
В пятимерном случае он приводит к следующему выражению для ошибки аппроксимации
𝜏
𝜏 =
𝜆
2 12
(︂ 𝜕
4
𝑢
𝜕𝑥
4 1
+
𝜕
4
𝑢
𝜕𝑥
4 2
+
𝜕
4
𝑢
𝜕𝑥
4 3
+
𝜕
4
𝑢
𝜕𝑥
4 4
+
𝜕
4
𝑢
𝜕𝑥
4 5
)︂
+ 𝑂
(︀𝜆
4
)︀ .
(3.5.1)
При подстановке
𝑢
5D
a получаем
𝜏 =
7 108
𝜆
2
(︁
−24 − 8𝑥
4
𝑥
5
sin 𝑥
5
+ 8𝑥
2
cos 𝑥
2

− (sin 𝑥
2
+ 𝑥
4
sin 𝑥
5
)
(︀𝑥
2 1
+ 𝑥
2 2
1   2   3   4   5   6   7   8   9   10   ...   13


написать администратору сайта