Решение упражнения 1.28 из SICP

15 October, 2007 (20:14) | Решения упражнений

Это упражнение немного заковыристо сформулировано, но сложным не является.

Вот код модифицированной процедуры expmod, которая сигнализирует о нахождении нетривиального квадратного корня из 1, возвращая 0 (как и рекомендовано в условии).

(define (square x) 
  (* x x))
(define (apply-trivial-check k m r) 
  (if (and (not (= k 1)) 
           (not (= k (- m 1))) 
           (= r 1)) 
      0 
      r))
(define (remainder-or-trivial k m) 
  (apply-trivial-check k m (remainder (square k) m)))
(define (modified-expmod base exp m) 
  (cond ((= exp 0) 1) 
        ((even? exp) 
         (remainder-or-trivial (modified-expmod base (/ exp 2) m) m)) 
        (else 
         (remainder (* base (modified-expmod base (- exp 1) m)) 
                    m))))

Теперь запишем всю процедуру, реализующую тест Миллера-Рабина. Она прогоняет для заданного числа n модифицированный тест Ферма на всех числах a<n и подсчитывает количество успехов. Если их больше половины, то число в соответствии с тестом Миллера-Рабина идентифицируется как простое. В противном случае оно составное.

(define (miller-rabin-test n) 
  (define (miller-rabin-iteration a t n) 
    (define (try-it a) 
      (= (modified-expmod a (- n 1) n) 1)) 
    (if (= a n) 
        (> t (/ n 2)) 
        (if (try-it a) 
            (miller-rabin-iteration (+ a 1) (+ t 1) n) 
            (miller-rabin-iteration (+ a 1) t n)))) 
  (miller-rabin-iteration 1 0 n))

Проверим написанную процедуру на нескольких простых и составных числах:

> (miller-rabin-test 5) 
#t 
> (miller-rabin-test 15) 
#f 
> (miller-rabin-test 97) 
#t 
> (miller-rabin-test 121) 
#f 
> (miller-rabin-test 1003) 
#f 
> (miller-rabin-test 1009) 
#t 
> (miller-rabin-test 100003) 
#t 
> (miller-rabin-test 100005) 
#f

Теперь проверим первые 6 чисел Кармайкла:

> (miller-rabin-test 561) 
#f 
> (miller-rabin-test 1105) 
#f 
> (miller-rabin-test 1729) 
#f 
> (miller-rabin-test 2465) 
#f 
> (miller-rabin-test 2821) 
#f 
> (miller-rabin-test 6601) 
#f

Как видим, определение простоты проходит корректно во всех случаях.

Comments

Comment from thror
Date: February 14, 2008, 11:45 am

И все же, несмотря на корректность работы процедуры на всех входных данных…

Общая схема работы теста Миллера-Рабина для числа n такова:

[1] Выбирается s случайных значений величины a из {1,…,n-1}
[2] Если одно из них свидетельствует о том, что n—составное, то n признается составным. Такой вывод всегда верный.
[3] Если в ходе s попыток не было таких свидетельств, то предполагается, что нет причин считать число n составным, так что оно считается простым.

Справедливы следующие теоремы (Кормен, Лейзерсон, Ривест, Штайн. Алгоритмы—построение и анализ, 2-е издание. Теоремы 31.38, 31.39)

Если n—нечетное составное число, то количество свидетельств того, что n—составное, не меньше (n-1)/2.

При любом нечетном n > 2 и положительном целом s вероятность того, что тест Миллера-Рабина выдаст неправильный результат не превышает 1/2**s.

Comment from thror
Date: February 14, 2008, 11:48 am

В дополнение два варианта, реализующих этот тест.

;; SICP

(define (square x)
  (* x x))

(define (sqrmod x m)
  (let ((y (remainder (square x) m)))
    (if (and (= y 1) (not (= x 1)) (not (= x (- m 1))))
        0
        y)))

(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp) (sqrmod (expmod base (/ exp 2) m) m))  
        (else (remainder (* base (expmod base (- exp 1) m)) m))))

(define (witness? a n)
  (not (= (expmod a (- n 1) n) 1)))

(define (miller-rabin-test n)
  (not (witness? (+ 1 (random (- n 1))) n)))

(define (fast-prime? n times)
  (cond ((= times 0) true)
        ((miller-rabin-test n) (fast-prime? n (- times 1)))
        (else false)))

Comment from thror
Date: February 14, 2008, 11:49 am

;; CLRS

(define (square x)
  (* x x))

(define (sqrmod x m)
  (remainder (square x) m))

(define (ntz n)
  (cond ((or (= n 0) (odd? n)) 0)
        (else (+ 1 (ntz (/ n 2))))))

(define (ntz n)
  (define (iter m r)
    (cond ((or (= m 0) (odd? m)) r)
          (else (iter (/ m 2) (+ r 1)))))

  (iter n 0))

(define (rtz n)
  (cond ((or (= n 0) (odd? n)) n)
        (else (rtz (/ n 2)))))

(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp) (sqrmod (expmod base (/ exp 2) m) m))  
        (else (remainder (* base (expmod base (- exp 1) m)) m))))

(define (witness? a n)
  (define (iter x t)
    (let ((y (remainder (square x) n)))
      (cond ((= t 0) (not (= x 1)))
            ((and (= y 1) (not (= x 1)) (not (= x (- n 1)))) true)
            (else (iter y (- t 1))))))
  
  (iter (expmod a (rtz (- n 1)) n) (ntz (- n 1))))

(define (miller-rabin-test n)
  (not (witness? (+ 1 (random (- n 1))) n)))

(define (fast-prime? n times)
  (cond ((= times 0) true)
        ((miller-rabin-test n) (fast-prime? n (- times 1)))
        (else false)))

Comment from thror
Date: February 14, 2008, 11:54 am

Один из них, как ясно из комментариев, “навеян” SICP, другой—CLRS =)

Оба приведены к общей схеме, и соответствующие функции названы одинаково, для большей ясности.

–thror

Comment from Sergey Khenkin
Date: February 21, 2008, 11:52 pm

Спасибо за комментарии. К сожалению, долго не возвращался к этому упражнению. Вернувшись, перечитал соответствующий раздел SICP, а также главу из замечательной книги Кормена, Лейзерсона и Ривеста “Алгоритмы: построение и анализ”.

В итоге хочу сказать, что наши подходы во многом сходны, хотя есть и одно очень существенное различие. Мне больше нравится разбивка на процедуры, сделанная не мной, а throrом (в первом варианте по мотивам SICP), поэтому я буду объяснять в ее терминах.

Действительно, в книге Кормена и др. тест Миллера-Рабина дается как вероятностный тест. Запущенный s раз на случайных значениях a<n он ошибается с вероятностью не более 1/2s в силу того, что если n не простое, то модифицированная процедура expmod, описанная выше, обнаружит нетривиальный квадратный корень из 1 по модулю n по крайней мере для половины a<n. Именно этот подход положен в основу проверки на простоту, предложенной throrом. И, как мне сейчас кажется, именно так на практике и используется тест Миллера-Рабина.

Но есть и другой ход мыслей. В условии упражнения сказано, что в нем идет речь о тесте “который невозможно обмануть”. Из этого можно заключить, что речь идет не о вероятностном тесте, а о тесте, который достоверно дает ответ на вопрос, является ли число простым. Мы знаем, что что если n не простое, то модифицированная процедура expmod, описанная выше, обнаружит нетривиальный квадратный корень из 1 по модулю n по крайней мере для половины a (define (fast-prime? n) (define (count-witnesses from) (if (= from n) 0 (if (witness? a n) (+ 1 (count-witnesses (+ from 1))) (count-witnesses (+ from 1))))) (< (count-witnesses 1) (/ (- n 1) 2)))

Я полагаю, что этот вариант решения также имеет право на жизнь и является, по сути, правильным.
Хотя с практической точки зрения, скорее всего, он менее полезен и применим.

Comment from Sergey Khenkin
Date: February 22, 2008, 9:03 am

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

(define (fast-prime? n)
  (define (any-witness? from)
    (cond ((= from n) false)
          ((witness? from n) true)
          (else (any-witness? (+ from 1)))))
  (not (any-witness? 1)))

Однако это все же далеко от практического применения.

Ведь определение простоты числа этим методом требует порядка Θ(n log n) шагов, что для больших чисел n неприемлемо. В то же время вероятностный тест, предложенный thror-ом, для s попыток требует только Θ(s log n) шагов и дает корректный результат с вероятностью не меньше, чем 1-2-s.

Comment from gorilych
Date: September 27, 2008, 7:45 am

кхм. Я нигде не вижу сравнения результата modified-expmod с нулём. Таким образом, вся ваша модификация коту под хвост – решение неправильное!

Общий алгоритм решения таков:

1) Берём n/2+1 чисел от 1 до n.
2) Последовательно вычисляем для каждого числа a из этого списка (modified-expmod a (- n 1) n)=res
Если res=0 – это значит, что мы нашли “нетривиальный квадратный корень из 1 по модулю n” – n – составное (однозначно)
Если res!=1 – это значит, что необходимое условие простоты n не удовлетворено, то есть опять n – составное (однозначно)
Если res=1 – продолжаем тестировать оставшиеся числа.
3) Если ни одно из чисел не выявило “нетривиальный квадратный корень из 1 по модулю n” – n – простое.

Comment from gorilych
Date: September 27, 2008, 8:16 am

а, ну да, оно и не должно сравниваться с нулём. однако у вас в miller-rabin-test подсчитывается, сколько раз результат modified-expmod был не равен единице. Достаточно одного раза – это уже гарантирует, что число n – составное.

Comment from трор
Date: November 2, 2008, 4:05 pm

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

Итак, задача в том, чтоб потратив к*лог(н) действий выяснить-составное число, или простое?

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

Витнесс, или на латинском-свидетель, в моих процедурах возвращает правду, когда число не простое. В этом вы легко убедитесь и сами. Для составного числа модифицированный экспмод вернет либо ноль, либо какое нибудь не равное единице число. Но не единицу.

А далее, следуем простой схеме-если к свидетелей не утверждают, что число составное, то считаем его простым. Как в жизни.

Откуда же взялось это загадочное н/2, которое способно даже грамотных (судя по 1.45) людей ввести в заблуждение? Зачем нужна эта громадная цифра? А для чисел в четыре тясячи бит длиной она действительно потрясает воображение.

Ответ кроется в простой идее, которую не трудно доказать-для нечетного составного числа н не менее (н-1)/2 свидетелей того, что оно составное.

Это обеспечивает экспоненциальное убывание вероятности ошибки при росте к.

Для сергея-на практике именно так он и реализован. Вы можете, к примеру, обратиться к старым реализациям библиотек шифрования из 4.4бсд-не помню сейчас точно, или План9-а там очень хорошая идея как модифицируют алгоритм на практике.

Вообщем-на самом деле то не слишком важно, правильное решение или нет-важно то, на какие идеи оно нас наталкивает, и о чем заставляет задуматься.

Comment from трор
Date: November 4, 2008, 10:05 am

вы можете взглянуть на код полной криптосистемы РСА (включающей свозможность генерировать длинный ключ случайно-но прежде чем утверждать что ключ получается не слишком случайный прочтите пожалуйста кнута 3.2.2 и упражнение 18 к ней, а также 3.3.2-критерии серий, и не говорите, что я не в курсе=), позволяет шифровать и расшифровывать сообщения, создавать цифровую подпись, создавать подписанные сообщения, и проверять подписи, и предоставлющую функции взлома приватного ключа по публичному, основанные как на классическом методе разложения, так и на эвристическом ро-методе полларда. На моей странице http://www.uic.nnov.ru/~kuaa29
–трор

Write a comment