Алгоритм разложения простого числа на сумму двух квадратов - Алгебра - Страница 2

@x3mEn · 07.02.2013, 01:10 **[ТС]**

Студворк — интернет-сервис помощи студентам

Ок, ввожу следующий фактор.
Вот описание этого супер-пупер быстрого метода.
Кто сможет пересказать его человеческим языком?

Primes congruent to 1 mod 4
We begin by showing that if p is a prime and p = 1 (mod 4), then there is a number a such that a² = -1 (mod p). We say that -1 is a quadratic residue mod p. One way to see this is the following. If a, b, and c are three integers and p does not divide a, then ab = ac (mod p) implies that b = c (mod p). The reason is that the conditions imply that p divides ab-ac and p does not divide a so that p divides b-c, which is the conclusion. Now for any a such that 0 < a < p the numbers a, 2a, 3a, ..., (p-1)a are all non-congruent mod p and therefore must be congruent to 1, 2, 3, ..., p-1 in some order. By taking the product of all the integers in the two sets, we see that a^p-1(p-1)! = (p-1)! (mod p). Since p divides no factor
of (p-1)!, it can be canceled to conclude that a^p-1 = 1 (mod p) (Fermat's little theorem).

Next we note that for any a not divisible by p, if we let b=a^(p-1)/2, then b² = 1 (mod p). This means that p divides (b²-1) = (b-1)(b+1) which means that a^(p-1)/2 = b = ±1 (mod p).

Can a^(p-1)/2=1 for all a? No, it cannot. I will not give all the details here, but the underlying reason is that arithmetic mod a prime is really like ordinary arithmetic in nearly all respects. In particular, the usual argument that a polynomial of degree d has at most d distinct roots is valid. If f(x) is a polynomial, then f(a) = 0 (mod p) if and only if x-a divides f(x), with division taking place mod p. This is not true for non-primes, by the way. For example, x²-1 = (x-1)(x+1) has four distinct roots mod 8 because 2 can divide one of the factors and 4 can divide the other. Getting back to the main argument, the polynomial x^(p-1)/2-1 can have at most (p-1)/2 roots as can the other factor of x^p-1-1 = (x^(p-1)/2-1) (x^(p-1)/2+1). Since the polynomial does have p-1 roots (namely, 1, 2, ..., p-1), it follows that exactly (p-1)/2 of the factors must satisfy the equation x^(p-1)/2 = -1 (mod p). Up to now, this argument is valid for all odd primes. The next step is limited to primes p, for which 4 divides p-1.

If a^(p-1)/2 = -1 (mod p) and we let b=a^(p-1)/4, then b² = -1 (mod p), as desired.

So the first step of the process is to find a square root of -1 mod p. We have shown, in fact, that exactly half of the numbers 1, ..., p-1 have the desired property. I do not know a deterministic process for finding such a number, but they appear to be distributed randomly so you can just try numbers at random until you find one and
each trial has a 50% chance of success. It is not hard to show that the smallest solution must be a prime so that you can just try primes 2, 3, 5, 7, 11,... until you succeed. Even if you are trying 1000 digit numbers it is unlikely that you will not find a solution before the 2000th prime, so for computational purposes, it is probably worth keeping a list of the first 2000 or so primes.

Actually, you can do better than that if a higher power of 2 divides p-1. If 2^k divides p-1 and we let m=(p-1)/2^k, then as soon as a^m $https://www.cyberforum.ru/cgi-bin/latex.cgi?\neq$ 1 (mod p), it follows that one among the numbers a, a², a⁴ will be a square root of -1 mod p. The probability of a number chosen at random satisfying a^m = 1 (mod p) is only 1 in 2^k-1.

Now a²+1² is a multiple of p, say a²+1²=kp. Since a <= p-1, we see that a²+1 < p² so that k < p. What we will now describe is a procedure that continues to produce solutions of a²+b²=kp, but with smaller values of k until k=1 and we are through. Assuming we have such a pair a and b, let a' and b' be the absolutely least residues a mod k and b mod k, respectively. This means that a' and b' are integers in the range (-k/2,k/2]. It follows that a'²+b'² <= (k/2)²+(k/2)² = k²/2. Also, a'²+b'² = a²+b² = 0 (mod k) since a = a' (mod k) and b = b' (mod k).

Thus a'²+b'²=mk with m < k. Multiplying, we find that (a²+b²)(a'²+b'²) = mk²p. Now we have the identity,

(a²+b²)(a'²+b'²) = (aa'+bb')²+(ab'-ba')² (*)

Since aa'+bb' = a²+b² = 0 (mod k) and ab'-ba' = ab-ba = 0 (mod k) it follows that both terms on the right hand side of (*) are divisible by k and so

((aa'+bb')/k)² + ((ab'-ba')/k)² = mp

and m < k. It can be shown that the number of reduction steps is bounded above by log₂ p.

iifat · 07.02.2013, 01:41

Ну, примерно так.
Сначала находим $https://www.cyberforum.ru/cgi-bin/latex.cgi?a^2\equiv-1(mod p)$ . Поскольку $https://www.cyberforum.ru/cgi-bin/latex.cgi?a^{\frac{p-1}2\equiv\pm1(mod p)$ и притом в равных количествах, можно пробовать возводить в степень 2, 3, и т.п., достаточно скоро должно найтись число. В результате имеем $https://www.cyberforum.ru/cgi-bin/latex.cgi?a^2+1=kp$ .
После этого начинаем процесс спуска. $https://www.cyberforum.ru/cgi-bin/latex.cgi?a^2+b^2=kp$ . Берём $https://www.cyberforum.ru/cgi-bin/latex.cgi?a'=a\,mod\,k,\,b'=b\,mod\,k$ и переходим к числам $https://www.cyberforum.ru/cgi-bin/latex.cgi?\frac{aa'+bb'}k,\frac{ab'-ba'}k$ -- и так пока не дойдём до p.

Добавлено через 5 минут
Ах да, искать надо среди простых чисел. Таблицы простых до 2000 должно хватить.
Также, если p-1 делится на более высокую степень двойки, начальный шаг можно ускорить.

@x3mEn · 07.02.2013, 02:22 **[ТС]**

Сообщение от iifat

Ну, примерно так.
Сначала находим $https://www.cyberforum.ru/cgi-bin/latex.cgi?a^2\equiv-1(mod p)$ .

Вот начнём с этого. Как найти $https://www.cyberforum.ru/cgi-bin/latex.cgi?a^2$ ? Ну воть хотя бы для нашего знакомого p = 1844674407370954349.
Это какая математика нужна, чтоб найти такое число? 128-битная, что ли?

Добавлено через 35 минут
Если имеется в виду, что нужно вычислить
2^922337203685477174 (mod 1844674407370954349), сравнить с -1
3^922337203685477174 (mod 1844674407370954349), сравнить с -1
5^922337203685477174 (mod 1844674407370954349), сравнить с -1
и т.д. пока сравнение не будет равенством, то пусть меня покрасят, у меня нет ни малейшего представления, как это можно посчитать и на каких дисках можно хранить эти петабайты цифр.

iifat · 07.02.2013, 04:45

Ну, напоминаю: всё ж по модулю p. 128-битной уж точно хватит; возможно, хватит и 64, хотя сходу что-то не придумывается. И при чём тут петабайты цифр? Эта с виду громадная степень -- 64- возведения в квадрат плюс не более чем столько умножений.
Куда-то я потерял своего Кнута, ёлки. Надо поискать.

Добавлено через 5 минут
Ну, например: $https://www.cyberforum.ru/cgi-bin/latex.cgi?(ax+b)(cx+d)=acx^2+(ad+bc)x+bd$ . Если $https://www.cyberforum.ru/cgi-bin/latex.cgi?x=2^{32}$ , то можно считать остатки по частям. Думаю, вполне можно остаться в рамках 64 бит.

@x3mEn · 07.02.2013, 08:52 **[ТС]**

Сообщение от iifat

Куда-то я потерял своего Кнута, ёлки. Надо поискать.

Не надо меня бить, я хороший

Если я правильно понимаю, возведение в степень нужно сделать по модулю. В результате мы получим +1 или -1.
Ок, предположим взяли за основу 2, получилось -1.
2^922337203685477174 + 1 = k * 1844674407370954349
Как найти k? Тоже по модулю p?

Признаюсь, я вот не силён в этом разделе математики, в делении с остатками.
Киньте, что ли какой-нибудь ссылкой на эту тему.

iifat · 07.02.2013, 09:16

Ну, я всё больше по книгам, насчёт ссылок вряд ли помогу. Теория чисел, арифметика по модулю -- думаю, найти будет несложно.
Насчёт k -- не по модулю. Гарантируется только, что оно будет меньше p. Пожалуй, не сильно это радует, конечно.

Добавлено через 1 минуту
Хм. А то и того не гарантируется... И что тут придумать, не очень понимаю...

@x3mEn · 07.12.2013, 19:54 **[ТС]**

Для того, чтобы сравнить 2^922337203685477174 (mod 1844674407370954349) с -1
(в общем случае b^n mod k с -1)
нужно посчитать, в какую степень нужно возвести k, чтобы максимально приблизиться к b^n.
Нужно вычислить m = [n * ln(b) / ln(k)]+1
В данном случае m = [922337203685477174*ln(2)/ln(1844674407370954349)]+1 = 15200502829552871
Теперь нужно представить m в виде суммы степеней 2.
Для этого достаточно представить число m в двоичной системе:
15200502829552871(10) = 1101100000000011000110110110111001000000 01010011100111(2)
Получается, что 15200502829552870 = 2^53 + 2^52 + 2^50 + 2^49 + 2^39 + 2^38+ 2^34 + 2^33 + 2^31 + 2^30 + 2^28 + 2^27 + 2^25 + 2^24 + 2^23 + 2^20 + 2^12 + 2^10 + 2^7 + 2^6 + 2^5 + 2^2 + 2^1 + 2^0
Возвести k в степень 15200502829552870 эквивалентно посчитать произведение k в степенях
2^53, 2^52, 2^50, 2^49, 2^39, 2^38+ 2^34, 2^33, 2^31, 2^30, 2^28, 2^27, 2^25, 2^24, 2^23, 2^20, 2^12, 2^10, 2^7, 2^6, 2^5, 2^2, 2^1, 2^0
Поскольку k^(2^(n+1)) = (k^2^n)^2, подсчёт всех k в указанных выше степенях потребует всего 53 операции умножения.

Высчитывать все значения с точностью до последнего знака не потребуется.
Результат mod k - это число от 0 до k-1
Заведомо известно, что k < 2^64, а значит достаточно вычислить 2^64 последних бита числа k^m,
высчитать 2^64 последних бита числа b^n, добавить 1 и сравнить между собой два числа.

При произведении двух чисел < 2^64 получится число, не более 2^128.
Из полученного результата нужно взять 64 последних бита, а старшие разряды отбросить.
Таким образом за 53 шага мы найдём все 64-битные остатки для k во всех степенях 2^n, где n=0..53
0: 1844674407370954349
1: 11805916207174773353
2: 17089064145905393425
3: 1928239932235652897
4: 4845331104372818497
5: 6592285539504428161
6: 736195842894620929
7: 5896668691730649601
8: 13498404477061800961
9: 18347673574999803905
10: 1222870319323058177
11: 9228562819106807809
12: 8402366234573291521
13: 17999937209791447041
14: 2259030824146501633
15: 8138469837573849089
16: 4954072467344982017
17: 16435093416004026369
18: 12356717414724927489
19: 12411308189031596033
20: 17119283982236647425
21: 3425238381167116289
22: 12724366945076379649
23: 12050806473702244353
24: 7403391429021532161
25: 3354129005639892993
26: 16237874822795755521
27: 15253984670526734337
28: 16961141661923016705
29: 16628460754743328769
30: 975119380494942209
31: 1950238760989884417
32: 3900477521979768833
33: 7800955043959537665
34: 15601910087919075329
35: 12757076102128599041
36: 7067408130547646465
37: 14134816261095292929
38: 9822888448481034241
39: 1199032823252516865
40: 2398065646505033729
41: 4796131293010067457
42: 9592262586020134913
43: 737781098330718209
44: 1475562196661436417
45: 2951124393322872833
46: 5902248786645745665
47: 11804497573291491329
48: 5162251072873431041
49: 10324502145746862081
50: 2202260217784172545
51: 4404520435568345089
52: 8809040871136690177
53: 17618081742273380353
Жирным выделены остатки, которые понадобятся для расчёта финального результата.

Теперь по такому же принципу последовательно высчитаем 64-битные остатки произведения k=1844674407370954349 в степенях
2^53, 2^52, 2^50, 2^49, 2^39, 2^38+ 2^34, 2^33, 2^31, 2^30, 2^28, 2^27, 2^25, 2^24, 2^23, 2^20, 2^12, 2^10, 2^7, 2^6, 2^5, 2^2, 2^1, 2^0:
k^2^0: 1844674407370954349
*(k^2^1): 4353431600858879157
*(k^2^2): 649142647038184197
*(k^2^5): 8384006120983589253
*(k^2^6): 15156778244432133765
*(k^2^7): 4908112718770280581
*(k^2^10): 10089594507726284933
*(k^2^12): 5230735132893456517
*(k^2^20): 7482045653221451909
*(k^2^23): 14062182855346563205
*(k^2^24): 3154094951022012549
*(k^2^25): 14191678424040679557
*(k^2^27): 4478960772764215429
*(k^2^28): 3212039167769127045
*(k^2^30): 450195757002467461
*(k^2^31): 13373253009178699909
*(k^2^33): 9725249796754974853
*(k^2^34): 2429243371907524741
*(k^2^38): 14820349090315184261
*(k^2^39): 2709072379711400069
*(k^2^49): 15420482327964625029
*(k^2^50): 3949814077051971717
*(k^2^52): 13407373294530013317
*(k^2^53): 13875747655776544901

13875747655776544901(10) = 1001101000001101001111110010111001001100 011010101000011011010(2)
Дополним до 64 бит ведущими нулями
0001001101000001101001111110010111001001 100011010101000011011010
Это и есть последние 64 бита числа 1844674407370954349^15200502829552871

Проделав похожую операцию над числом b^n+1, получим 64-битный остаток.
В случае b=2 особо и высчитывать не надо, остаток будет равен
0000000000000000000000000000000000000000 000000000000000000000001

В данном случае 64-битные остатки не равны, значит берём следующее значение b из множества простых (3) и проделываем операцию повторно.
Итого на проверку b^n mod k = -1 нужно потратить не более 256 операций умножения (если не брать во внимание вычисление m = [n * ln(b) / ln(k)]+1):
- не более 64 операций на получение 64-битных остатков чисел k^2^n
- не более 64 операций на получение 64-битного остатка произведений нужных k^2^n
- те же не более 128 операций для числа b^n

@barbalion · 15.05.2019, 16:51

вот реализация описанного выше алгоритма на python, если кому-то пригодится:

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
def primeAsSumOfSquares(p):
  if p % 4 != 1:
    return None
  p2 = (p-1)//2
  for a in allPrimes():
    if pow(a, p2, p) == p-1:
      break
  a, b = pow(a, p2//2, p), 1
  while True:
    k = (a*a + b*b) // p
    if k == 1:
      return (b, a)
    a1, b1 = a % k, b % k
    a, b = (a * a1 + b * b1) // k, (a * b1 - b * a1) // k

Реализация allPrimes() остается на домашнее задание

@barbalion · 15.05.2019, 23:19

патч в 13 строчке, а то бесконечные циклы иногда встречаются (правка сообщений мне пока недоступна):

Code
1
2
-    a1, b1 = a % k, b % k
+    a1, b1 = (a + k//2) % k - k//2, (b + k//2) % k - k//2

@x3mEn · 27.10.2020, 11:06 **[ТС]**

Реализация на С:

Кликните здесь для просмотра всего текста

C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
#include <math.h>
#include <inttypes.h>
 
#define xchgu64(a,b) \
do {uint64_t c = *a; *a = *b; *b = c;} while (0)
 
// 6542 primes less then 2^16 = 65536
#define SMALL_PRIMES_CNT 6542
uint32_t SmallPrimes[SMALL_PRIMES_CNT];
 
/* (a*b) % m */
//  m must be less than 2^63
static __inline__ uint64_t mod_mul(uint64_t a, uint64_t b, const uint64_t m)
{
   return (uint64_t)(((__uint128_t)a * b) % m);
}
 
/* (a^b) % m */
static __inline__ uint64_t mod_pow(uint64_t a, uint64_t b, const uint64_t m)
{
    uint64_t r = 1;
    while (b > 0) {
        if(b & 1)
            r = mod_mul(r, a, m);
        b = b >> 1;
        a = mod_mul(a, a, m);
    }
    return r;
}
 
static __inline__ uint64_t root4(const uint64_t n)
{
    uint64_t a, b, k = n/4;
    for (uint16_t i=0; ; i++) {
        uint16_t j = SmallPrimes[i];
        a = mod_pow(j, k, n);
        b = mod_mul(a, a, n);
        if (b == n-1)
            return a;
    }
}
 
static __inline__ void decompose_4kp1(uint64_t n, uint64_t *x, uint64_t *y)
{
    int64_t c, r, s;
    s = rintl(sqrtl(n));
    c = root4(n);
    r = n % c;
    while (c > s) {
      n = c;
      c = r;
      r = n % c;
    };
    (*x) = r;
    (*y) = c;
    if (*x > *y)
        xchgu64(x, y);
}

Заполнение массива SmallPrimes в качестве домашнего задания. )

Новые блоги и статьи Все статьи Все блоги /
Беседа с ИИ о программистах, недопускающих к созданию и правке кода генеративные ИИ и причины этого zorxor 21.09.2026 Раньше я радовался или получал некоторые эмоции, пусть небольшие, но всё же, от самого процесса написания кода, рекомпиляции и запуска, видя постепенное развитие программы и прочее. А теперь лень. . .	Мобильное приложение ColorStep pavlinmavlin 17.09.2026 Реализовал приложение Красный, Зеленый, Синий в Unity3d + c#. Название изменил на ColorStep. Приложение прошло модерацию и теперь доступно для скачивания. Делал его сам, шаг за шагом — и вот,. . .	Запрет дублирования строк в табличной части Maks 13.09.2026 Реализация из решения ниже выполнена на нетиповом справочнике "Нормы ТО" с табличной часть "Виды ТО", разработанного в КА2, со следующими реквизитами: - ВидТО (СправочникСсылка. ВидыТО); - ВидГСМ. . .	Скрипты Tampermonkey для CyberForum, ChatGPT, Claude и пр. Jin X 06.09.2026 Скрипты Tampermonkey для CyberForum, ChatGPT, Claude и пр. Работая с форумом и нейросетями в браузере часто хочется что-то подкорректировать или добавить какого-то функционала. Ниже прикреплён. . .
Программа опроса у.з. расходомера SLS-720F Argus19 02.09.2026 Программа опроса у. з. расходомера SLS-720F Программа опрашивает один раз в минуту три ультразвуковых расходомера SLS-720F через интерфейс RS-485 по протоколу Modbus RTU. Опрашиваются регистры. . .	Hyper-V: Компьютер должен поддерживать доверенный платформенный модуль 2.0. Maks 31.08.2026 При установке Windows 11 на виртуальную машину Hyper-V 2-го поколения вылезла такая ошибка: Решение: в параметрах виртуальной машины, в разделе "Безопасность" (Security) активировать флаг. . .	Архитектура биовида Стива в Майнкрафте: Зачем бонобо кубический каннибализм anaschu 30.08.2026 Кубический Вагинокапитализм в Minecraft: Математический инвариант ОДУ и рок Стивов-бонобо Главная задача разработанной «Модели Всего» — наглядно продемонстрировать наличие системной «судьбы». . .	Оттачиваю умение писать js программы. russiannick 30.08.2026 Проектом выходного дня стало написание Книги шифров Виженера. Итогом стала версия 200, синий туман. Синий туман назван так, потому что замораживает текст под собой. Нажатие синих кнопок управляют. . .

@x3mEn 2 / 2 / 0 Регистрация: 24.01.2013 Сообщений: 15
	07.02.2013, 01:10 [ТС]
	Ок, ввожу следующий фактор. Вот описание этого супер-пупер быстрого метода. Кто сможет пересказать его человеческим языком? Primes congruent to 1 mod 4 We begin by showing that if p is a prime and p = 1 (mod 4), then there is a number a such that a² = -1 (mod p). We say that -1 is a quadratic residue mod p. One way to see this is the following. If a, b, and c are three integers and p does not divide a, then ab = ac (mod p) implies that b = c (mod p). The reason is that the conditions imply that p divides ab-ac and p does not divide a so that p divides b-c, which is the conclusion. Now for any a such that 0 < a < p the numbers a, 2a, 3a, ..., (p-1)a are all non-congruent mod p and therefore must be congruent to 1, 2, 3, ..., p-1 in some order. By taking the product of all the integers in the two sets, we see that a^p-1(p-1)! = (p-1)! (mod p). Since p divides no factor of (p-1)!, it can be canceled to conclude that a^p-1 = 1 (mod p) (Fermat's little theorem). Next we note that for any a not divisible by p, if we let b=a^(p-1)/2, then b² = 1 (mod p). This means that p divides (b²-1) = (b-1)(b+1) which means that a^(p-1)/2 = b = ±1 (mod p). Can a^(p-1)/2=1 for all a? No, it cannot. I will not give all the details here, but the underlying reason is that arithmetic mod a prime is really like ordinary arithmetic in nearly all respects. In particular, the usual argument that a polynomial of degree d has at most d distinct roots is valid. If f(x) is a polynomial, then f(a) = 0 (mod p) if and only if x-a divides f(x), with division taking place mod p. This is not true for non-primes, by the way. For example, x²-1 = (x-1)(x+1) has four distinct roots mod 8 because 2 can divide one of the factors and 4 can divide the other. Getting back to the main argument, the polynomial x^(p-1)/2-1 can have at most (p-1)/2 roots as can the other factor of x^p-1-1 = (x^(p-1)/2-1) (x^(p-1)/2+1). Since the polynomial does have p-1 roots (namely, 1, 2, ..., p-1), it follows that exactly (p-1)/2 of the factors must satisfy the equation x^(p-1)/2 = -1 (mod p). Up to now, this argument is valid for all odd primes. The next step is limited to primes p, for which 4 divides p-1. If a^(p-1)/2 = -1 (mod p) and we let b=a^(p-1)/4, then b² = -1 (mod p), as desired. So the first step of the process is to find a square root of -1 mod p. We have shown, in fact, that exactly half of the numbers 1, ..., p-1 have the desired property. I do not know a deterministic process for finding such a number, but they appear to be distributed randomly so you can just try numbers at random until you find one and each trial has a 50% chance of success. It is not hard to show that the smallest solution must be a prime so that you can just try primes 2, 3, 5, 7, 11,... until you succeed. Even if you are trying 1000 digit numbers it is unlikely that you will not find a solution before the 2000th prime, so for computational purposes, it is probably worth keeping a list of the first 2000 or so primes. Actually, you can do better than that if a higher power of 2 divides p-1. If 2^k divides p-1 and we let m=(p-1)/2^k, then as soon as a^m $https://www.cyberforum.ru/cgi-bin/latex.cgi?\neq$ 1 (mod p), it follows that one among the numbers a, a², a⁴ will be a square root of -1 mod p. The probability of a number chosen at random satisfying a^m = 1 (mod p) is only 1 in 2^k-1. Now a²+1² is a multiple of p, say a²+1²=kp. Since a <= p-1, we see that a²+1 < p² so that k < p. What we will now describe is a procedure that continues to produce solutions of a²+b²=kp, but with smaller values of k until k=1 and we are through. Assuming we have such a pair a and b, let a' and b' be the absolutely least residues a mod k and b mod k, respectively. This means that a' and b' are integers in the range (-k/2,k/2]. It follows that a'²+b'² <= (k/2)²+(k/2)² = k²/2. Also, a'²+b'² = a²+b² = 0 (mod k) since a = a' (mod k) and b = b' (mod k). Thus a'²+b'²=mk with m < k. Multiplying, we find that (a²+b²)(a'²+b'²) = mk²p. Now we have the identity, *(a²+b²)(a'²+b'²) = (aa'+bb')²+(ab'-ba')² () Since aa'+bb' = a²+b² = 0 (mod k) and ab'-ba' = ab-ba = 0 (mod k) it follows that both terms on the right hand side of ()* are divisible by k and so ((aa'+bb')/k)² + ((ab'-ba')/k)² = mp and m < k. It can be shown that the number of reduction steps is bounded above by log₂ p. 1

iifat 2923 / 1953 / 215 Регистрация: 05.06.2011 Сообщений: 5,791
	07.02.2013, 01:41
	Ну, примерно так. Сначала находим $https://www.cyberforum.ru/cgi-bin/latex.cgi?a^2\equiv-1(mod p)$ . Поскольку $https://www.cyberforum.ru/cgi-bin/latex.cgi?a^{\frac{p-1}2\equiv\pm1(mod p)$ и притом в равных количествах, можно пробовать возводить в степень 2, 3, и т.п., достаточно скоро должно найтись число. В результате имеем $https://www.cyberforum.ru/cgi-bin/latex.cgi?a^2+1=kp$ . После этого начинаем процесс спуска. $https://www.cyberforum.ru/cgi-bin/latex.cgi?a^2+b^2=kp$ . Берём $https://www.cyberforum.ru/cgi-bin/latex.cgi?a'=a\,mod\,k,\,b'=b\,mod\,k$ и переходим к числам $https://www.cyberforum.ru/cgi-bin/latex.cgi?\frac{aa'+bb'}k,\frac{ab'-ba'}k$ -- и так пока не дойдём до p. Добавлено через 5 минут Ах да, искать надо среди простых чисел. Таблицы простых до 2000 должно хватить. Также, если p-1 делится на более высокую степень двойки, начальный шаг можно ускорить. 0

iifat 2923 / 1953 / 215 Регистрация: 05.06.2011 Сообщений: 5,791
	07.02.2013, 04:45
	Ну, напоминаю: всё ж по модулю p. 128-битной уж точно хватит; возможно, хватит и 64, хотя сходу что-то не придумывается. И при чём тут петабайты цифр? Эта с виду громадная степень -- 64- возведения в квадрат плюс не более чем столько умножений. Куда-то я потерял своего Кнута, ёлки. Надо поискать. Добавлено через 5 минут Ну, например: $https://www.cyberforum.ru/cgi-bin/latex.cgi?(ax+b)(cx+d)=acx^2+(ad+bc)x+bd$ . Если $https://www.cyberforum.ru/cgi-bin/latex.cgi?x=2^{32}$ , то можно считать остатки по частям. Думаю, вполне можно остаться в рамках 64 бит. 0

iifat 2923 / 1953 / 215 Регистрация: 05.06.2011 Сообщений: 5,791
	07.02.2013, 09:16
	Ну, я всё больше по книгам, насчёт ссылок вряд ли помогу. Теория чисел, арифметика по модулю -- думаю, найти будет несложно. Насчёт k -- не по модулю. Гарантируется только, что оно будет меньше p. Пожалуй, не сильно это радует, конечно. Добавлено через 1 минуту Хм. А то и того не гарантируется... И что тут придумать, не очень понимаю... 0

@x3mEn 2 / 2 / 0 Регистрация: 24.01.2013 Сообщений: 15
	07.12.2013, 19:54 [ТС]
	Для того, чтобы сравнить 2^922337203685477174 (mod 1844674407370954349) с -1 (в общем случае b^n mod k с -1) нужно посчитать, в какую степень нужно возвести k, чтобы максимально приблизиться к b^n. Нужно вычислить m = [n * ln(b) / ln(k)]+1 В данном случае m = [922337203685477174ln(2)/ln(1844674407370954349)]+1 = 15200502829552871 Теперь нужно представить m в виде суммы степеней 2. Для этого достаточно представить число m в двоичной системе: 15200502829552871(10) = 1101100000000011000110110110111001000000 01010011100111(2) Получается, что 15200502829552870 = 2^53 + 2^52 + 2^50 + 2^49 + 2^39 + 2^38+ 2^34 + 2^33 + 2^31 + 2^30 + 2^28 + 2^27 + 2^25 + 2^24 + 2^23 + 2^20 + 2^12 + 2^10 + 2^7 + 2^6 + 2^5 + 2^2 + 2^1 + 2^0 Возвести k в степень 15200502829552870 эквивалентно посчитать произведение k в степенях 2^53, 2^52, 2^50, 2^49, 2^39, 2^38+ 2^34, 2^33, 2^31, 2^30, 2^28, 2^27, 2^25, 2^24, 2^23, 2^20, 2^12, 2^10, 2^7, 2^6, 2^5, 2^2, 2^1, 2^0 Поскольку k^(2^(n+1)) = (k^2^n)^2, подсчёт всех k в указанных выше степенях потребует всего 53 операции умножения. Высчитывать все значения с точностью до последнего знака не потребуется. Результат mod k - это число от 0 до k-1 Заведомо известно, что k < 2^64, а значит достаточно вычислить 2^64 последних бита числа k^m, высчитать 2^64 последних бита числа b^n, добавить 1 и сравнить между собой два числа. При произведении двух чисел < 2^64 получится число, не более 2^128. Из полученного результата нужно взять 64 последних бита, а старшие разряды отбросить. Таким образом за 53 шага мы найдём все 64-битные остатки для k во всех степенях 2^n, где n=0..53 0: 1844674407370954349* 1: 11805916207174773353 2: 17089064145905393425 3: 1928239932235652897 4: 4845331104372818497 5: 6592285539504428161 6: 736195842894620929 7: 5896668691730649601 8: 13498404477061800961 9: 18347673574999803905 10: 1222870319323058177 11: 9228562819106807809 12: 8402366234573291521 13: 17999937209791447041 14: 2259030824146501633 15: 8138469837573849089 16: 4954072467344982017 17: 16435093416004026369 18: 12356717414724927489 19: 12411308189031596033 20: 17119283982236647425 21: 3425238381167116289 22: 12724366945076379649 23: 12050806473702244353 24: 7403391429021532161 25: 3354129005639892993 26: 16237874822795755521 27: 15253984670526734337 28: 16961141661923016705 29: 16628460754743328769 30: 975119380494942209 31: 1950238760989884417 32: 3900477521979768833 33: 7800955043959537665 34: 15601910087919075329 35: 12757076102128599041 36: 7067408130547646465 37: 14134816261095292929 38: 9822888448481034241 39: 1199032823252516865 40: 2398065646505033729 41: 4796131293010067457 42: 9592262586020134913 43: 737781098330718209 44: 1475562196661436417 45: 2951124393322872833 46: 5902248786645745665 47: 11804497573291491329 48: 5162251072873431041 49: 10324502145746862081 50: 2202260217784172545 51: 4404520435568345089 52: 8809040871136690177 53: 17618081742273380353 Жирным выделены остатки, которые понадобятся для расчёта финального результата. Теперь по такому же принципу последовательно высчитаем 64-битные остатки произведения k=1844674407370954349 в степенях 2^53, 2^52, 2^50, 2^49, 2^39, 2^38+ 2^34, 2^33, 2^31, 2^30, 2^28, 2^27, 2^25, 2^24, 2^23, 2^20, 2^12, 2^10, 2^7, 2^6, 2^5, 2^2, 2^1, 2^0: k^2^0: 1844674407370954349 (k^2^1): 4353431600858879157 (k^2^2): 649142647038184197 (k^2^5): 8384006120983589253 (k^2^6): 15156778244432133765 (k^2^7): 4908112718770280581 (k^2^10): 10089594507726284933 (k^2^12): 5230735132893456517 (k^2^20): 7482045653221451909 (k^2^23): 14062182855346563205 (k^2^24): 3154094951022012549 (k^2^25): 14191678424040679557 (k^2^27): 4478960772764215429 (k^2^28): 3212039167769127045 (k^2^30): 450195757002467461 (k^2^31): 13373253009178699909 (k^2^33): 9725249796754974853 (k^2^34): 2429243371907524741 (k^2^38): 14820349090315184261 (k^2^39): 2709072379711400069 (k^2^49): 15420482327964625029 (k^2^50): 3949814077051971717 (k^2^52): 13407373294530013317 (k^2^53): 13875747655776544901 13875747655776544901(10) = 1001101000001101001111110010111001001100 011010101000011011010(2) Дополним до 64 бит ведущими нулями 0001001101000001101001111110010111001001 100011010101000011011010 Это и есть последние 64 бита числа 1844674407370954349^15200502829552871 Проделав похожую операцию над числом b^n+1, получим 64-битный остаток. В случае b=2 особо и высчитывать не надо, остаток будет равен 0000000000000000000000000000000000000000 000000000000000000000001 В данном случае 64-битные остатки не равны, значит берём следующее значение b из множества простых (3) и проделываем операцию повторно. Итого на проверку b^n mod k = -1 нужно потратить не более 256 операций умножения (если не брать во внимание вычисление m = [n ln(b) / ln(k)]+1): - не более 64 операций на получение 64-битных остатков чисел k^2^n - не более 64 операций на получение 64-битного остатка произведений нужных k^2^n - те же не более 128 операций для числа b^n 0