;; hill.lisp - A very bad implementation of a hill cipher along with some
;; unnessesary, insecure PRNG key stuff.
;;
;; This file is free software under the terms of the GPL3 license. You can find
;; a copy of the license at: https://www.gnu.org/licenses/gpl-3.0.en.html

;; PRNG subroutines
(defparameter *system-rng-stream*
  (open "/dev/random" :element-type '(unsigned-byte 8)))

(defun number-bits (num)
  "Return the number of bits in NUM."
  (ceiling (log num 2)))

(defun get-random (&key (size 8) (rng-stream *system-rng-stream*))
  "Read SIZE random bytes from the systems high-quality random source."
  (let ((read-buf (make-array size)))
    (read-sequence read-buf rng-stream)
    (loop with out-val = 0
          for part across read-buf do
            (setq out-val (+ (ash out-val 8) part))
          finally (return out-val))))

(defun xorshift*64 (seed)
  "Generate the next number in the xorshift* sequence starting or continuing
from SEED. That is, SEED should either be an initial, non-zero seed or the last
number generated in the sequence."
  (setq seed (logxor seed (ash seed -12))
        seed (logxor seed (ash seed 25))
        seed (logxor seed (ash seed -27)))
  ;; #x2545F4914F6CDD1DULL
  (logand (* seed 2685821657736338717)
          #xFFFFFFFFFFFFFFFF))

(defun rand-clamp (num &key (bits 64) (min 0 has-min-p)
                         (max (1- (expt 2 bits)) has-max-p))
  "Clamp the randomly generated number NUM to between [MIN, MAX]. The random
numbers is assumed to be between 0 and 2^BITS - 1."
  (if (or has-min-p has-max-p)
      (+ (* (- max min) (/ num (1- (expt 2 bits)))) min)
      num))

(defclass xorshift*64-stream ()
  ((state :initarg :seed
          :initform 1
          :documentation "The last generated number. If no number has yet been
generated, this is the seed."))
  (:documentation "A stream of random numbers generated and xorshift*64
algorithm."))

(defun make-xorshift* (&key (seed (get-random)) (bits 64))
  (case bits
    (64 (make-instance 'xorshift*64-stream :seed seed))
    (t (error "Invalid prandom number size: ~d" bits))))

(defmethod next-int ((stream xorshift*64-stream) &key (max #xFFFFFFFFFFFFFFFF)
                                                   (min 0))
  "Get the next random integer between MIN and MAX from STREAM."
  (with-slots (state) stream
    (setf state (xorshift*64 state))
    (floor (rand-clamp state :bits 64 :max max :min min))))

(defparameter *default-rng* (make-xorshift*)
  "The default RNG source to use when generating matrices.")

;; Matrix subroutines
(defun mat-minor (transpose-mat i j)
  "Find the minor of MAT for I and J."
  (destructuring-bind (height width)
      (array-dimensions transpose-mat)
    (let ((minor (make-array (list (1- height)
                                   (1- width)))))
      (dotimes (row height minor)
        (dotimes (col width)
          (unless (or (= row i)
                      (= col j))
            (let ((out-row (if (> row i)
                               (1- row)
                               row))
                  (out-col (if (> col j)
                               (1- col)
                               col)))
              (setf (aref minor out-row out-col) (aref transpose-mat row col)))))))))

(defun cofactor-sgn (i j)
  "Return the sign of the cofactor at I and J."
  (expt -1 (+ i j)))

(defun cofactor (mat i j)
  "Find the cofactor for I and J in MAT."
  (* (cofactor-sgn i j) (det (mat-minor mat i j))))

(defun first-column-cofactors (mat)
  "Find the cofactors for the first column of MAT."
  (let* ((height (array-dimension mat 0))
         (out-arr (make-array height)))
    (dotimes (i height out-arr)
      (setf (aref out-arr i) (cofactor mat i 0)))))

(defun det2x2 (mat)
  "Find the determinate of a 2x2 matrix MAT."
  (let ((a (aref mat 0 0))
        (b (aref mat 0 1))
        (c (aref mat 1 0))
        (d (aref mat 1 1)))
    (- (* a d) (* b c))))

(defun det (mat &key first-column-cofactors)
  "Find the determinant of MAT. If the cofactors for the first column have
already been calculated, they can be supplied in FIRST-COLUMN-COFACTORS."
  (destructuring-bind (height width)
      (array-dimensions mat)
    (if (and (= height 2)
             (= width 2))
        (det2x2 mat)
        (loop for i below height
              when first-column-cofactors
                summing (* (aref mat i 0) (aref first-column-cofactors i))
              else
                summing (* (aref mat i 0) (cofactor mat i 0))))))

(defun map-mat (function mat)
  "Return a new matrix of the same dimensions as MAT, with each element being
the result of applying FUNCTION to each index of MAT. FUNCTION is passed two
values, (i j), where i and j are the row and column of the element."
  (destructuring-bind (height width)
      (array-dimensions mat)
    (let ((out-mat (make-array (list height width))))
      (dotimes (row height out-mat)
        (dotimes (col width)
          (setf (aref out-mat row col)
                (funcall function row col)))))))

(defun extended-gcd (n1 n2)
  "Return the greatest common denominator of N1 and N2. As a second value,
return a cons of (x . y) such char x * N1 + y * N2 = GCD."
  (if (zerop n1)
      (values n2 (cons 0 1))
      (multiple-value-bind (div rem) (floor n2 n1)
        (destructuring-bind (gcd (x . y))
            (multiple-value-list (extended-gcd rem n1))
          (values gcd (cons (- y (* div x)) x))))))

(defun modular-multiplicative-inverse (n m)
  "Find the modular multiplicative inverse of N with respect to M. That is, some
integer i such that (N * i) mod M = 1."
  (destructuring-bind (gcd (x . y))
      (multiple-value-list (extended-gcd n m))
    (declare (ignorable y))
    (unless (= 1 gcd)
      (error "N and M must be coprime. N: ~d, M: ~d" n m))
    (+ x m)))

(defun modular-inverse-det (mat modulo &key first-column-cofactors)
  "Calculate the inverse determinant of MAT under MODULO. If MODULO is nil, this
is just 1/det."
  (let ((det (det mat :first-column-cofactors first-column-cofactors)))
    (if modulo
        (modular-multiplicative-inverse det modulo)
        (/ det))))

(defun invert-2x2-mat (mat &key modulo)
  "Invert the 2x2 matrix MAT."
  (let ((a (aref mat 0 0))
        (b (aref mat 0 1))
        (c (aref mat 1 0))
        (d (aref mat 1 1))
        (ood (modular-inverse-det mat modulo)))
    (make-array '(2 2)
                :initial-contents (list (list (* ood d) (* ood (- b)))
                                        (list (* ood (- c)) (* ood a))))))

(defun invert-mat (mat &key modulo)
  "Invert MAT. This naively assumes that MAT is not singular."
  (destructuring-bind (height width)
      (array-dimensions mat)
    (if (and (= width 2)
             (= height 2))
        (invert-2x2-mat mat :modulo modulo)
        (let* ((first-column-cofactors (first-column-cofactors mat))
               (one-over-det (modular-inverse-det mat modulo
                                                  :first-column-cofactors
                                                  first-column-cofactors)))
          (map-mat (lambda (i j)
                     ;; this calculates 1/det * adjugate[i][j]
                     (* one-over-det
                        (if (and first-column-cofactors (zerop i))
                            (aref first-column-cofactors j)
                            (cofactor mat j i))))
                   mat)))))

(defun transpose-mat (mat)
  "Transpose MAT."
  (map-mat (lambda (i j)
             (aref mat j i))
           mat))

(defun dot-row-col (mat1 mat2 row col)
  "Take the dot (scalar) product of ROW of MAT1 and COL of MAT2."
  (let ((sum 0))
    (dotimes (n (array-dimension mat2 0) sum)
      (setq sum (+ sum (* (aref mat1 row n)
                          (aref mat2 n col)))))))

(defun *mm (mat1 mat2 &key modulo)
  "Multiply MAT1 by MAT2. This naively assumes that these matrices can be
multiplied without issue. The result will be taken under MODULO, if it is
non-nil."
  (let* ((width (array-dimension mat1 0))
         (height (array-dimension mat2 1))
         (out-mat (make-array (list height width))))
    (dotimes (i height out-mat)
      (dotimes (j width)
        (let ((dot (dot-row-col mat1 mat2 i j)))
          (setf (aref out-mat i j)
                (if modulo
                    (mod dot modulo)
                    dot)))))))

(defun *mv (mat vec &key modulo)
  "Multiply MAT by VEC. This naively assumes that these can be multiplied
without issue. The result will be taken under MODULO, if it is non-nil."
  (let* ((width (array-dimension mat 1))
         (height (array-dimension mat 0))
         (out-vec (make-array height)))
    (dotimes (row height out-vec)
      (setf (aref out-vec row)
            (loop for col below width
                  summing (* (aref vec col) (aref mat row col)) into sum
                  finally (return (if modulo
                                      (mod sum modulo)
                                      sum)))))))

(defun generate-matrix (width height &key (rng *default-rng*) max)
  "Generate a random matrix of size WIDTH x HEIGHT. Use RNG to get the random
data for this matrix."
  (map-mat (lambda (i j)
             (declare (ignorable i j))
             (if max
                 (next-int rng :max max)
                 (next-int rng)))
           (make-array (list height width))))

(defun generate-invertible-matrix (width height &key (rng *default-rng*) max)
  "Generate a random invertible matrix of size WIDTH x HEIGHT. Use RNG to get
the random data for this matrix. This generates a strictly diagonally dominant
matrix."
  (let ((rand-mat (generate-matrix width height :rng rng :max max)))
    (dotimes (row height rand-mat)
      (let ((sum 0))
        (dotimes (col width)
          (setq sum (+ sum (abs (aref rand-mat row col)))))
        (setf (aref rand-mat row row) sum)))))

;; Base64 subroutines
(defparameter +base64-alphabet+
  "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/"
  "Lookup table for base64 encode.")

(defparameter +inverse-base64-alphabet+
  (loop with alph = (make-array 123 :initial-element 0)
        for i below (length +base64-alphabet+)
        for char = (aref +base64-alphabet+ i)
        do (setf (aref alph (char-code char)) i)
        finally (return alph))
  "Inverse lookup table for base64 decode.")

(defun base64-encode-block (b1 &optional b2 b3)
  "Base64 encode B1, B2, and B3, each of which should be a single byte. The
result of this function is a list of the 4 encoded characters."
  ;; allow B2 and B3 to be nil, so we can't just check if they where passed, we
  ;; need to know if they where nil when the function was called
  (let ((has-b2-p b2)
        (has-b3-p b3))
    (unless b2
      (setq b2 0))
    (unless b3
      (setq b3 0))
    (let ((e1 (ash (logand b1 #b11111100) -2))
          (e2 (logior (ash (logand b1 #b11) 4)
                      (ash (logand b2 #b11110000) -4)))
          (e3 (logior (ash (logand b2 #b1111) 2)
                      (ash (logand b3 #b11000000) -6)))
          (e4 (logand b3 #b111111)))
      (cond
        ((and has-b2-p has-b3-p)
         (list (aref +base64-alphabet+ e1)
               (aref +base64-alphabet+ e2)
               (aref +base64-alphabet+ e3)
               (aref +base64-alphabet+ e4)))
        (has-b2-p
         (list (aref +base64-alphabet+ e1)
               (aref +base64-alphabet+ e2)
               (aref +base64-alphabet+ e3)
               #\=))
        (t
         (list (aref +base64-alphabet+ e1)
               (aref +base64-alphabet+ e2)
               #\= #\=))))))

(defun base64-encode (data)
  "Encode DATA, a vector of numbers, as base64."
  (loop with output = ()
        for i upfrom 0 by 3 below (- (length data) 2)
        do (setf output (nconc (nreverse (base64-encode-block
                                          (aref data i)
                                          (aref data (+ i 1))
                                          (aref data (+ i 2))))
                               output))
        finally
           (return
             (coerce (if (< i (length data))
                         (let ((b1 (aref data i))
                               (b2 (when (< (1+ i) (length data))
                                     (aref data (1+ i)))))
                           (nreverse (nconc (nreverse
                                             (base64-encode-block b1 b2))
                                            output)))
                         (nreverse output))
                     'string))))

(defun base64-decode-block (e1 e2 e3 e4)
  "Base64 decode E1, E2, E3, and E4, each of which is a single character. The
result of this function is a list of 3 bytes."
  (let ((i1 (aref +inverse-base64-alphabet+ (char-code e1)))
        (i2 (aref +inverse-base64-alphabet+ (char-code e2)))
        (i3 (aref +inverse-base64-alphabet+ (char-code e3)))
        (i4 (aref +inverse-base64-alphabet+ (char-code e4))))
    (list (logior (ash i1 2)
                  (ash (logand i2 #b110000) -4))
          (unless (eq e3 #\=)
            (logior (ash (logand i2 #b1111) 4)
                    (ash (logand i3 #b111100) -2)))
          (unless (eq e4 #\=)
            (logior (ash (logand i3 #b11) 6)
                    i4)))))

(defun base64-decode (data)
  "Decode DATA, a string, as base64."
  (coerce (loop for i upfrom 0 by 4 below (length data)
                append (delete-if 'null
                                  (base64-decode-block (aref data i)
                                                       (aref data (+ i 1))
                                                       (aref data (+ i 2))
                                                       (aref data (+ i 3)))))
          'vector))

;; Hill subroutines
(defun nth-hill-block (plain n &key (size 4))
  "Get the Nth hill block in PLAIN. If there are not enough characters in PLAIN,
pad the block with zeros. Each block is SIZE characters."
  (let* ((start (* size n))
         (end (+ start size)))
    (when (< (length plain) end)
      (setq end (+ start (- (length plain) start))))
    (cond
      ((> start (length plain))
       (make-array size :initial-element 0))
      ((>= end (length plain))
       (loop with out-vec = (make-array size :initial-element 0 :fill-pointer 0)
             for elem across (subseq plain start (length plain)) do
               (vector-push elem out-vec)
             finally
                (setf (fill-pointer out-vec) size)
                (return out-vec)))
      (t
       (subseq plain start end)))))

(defun hill-block-count (plain &key (size 4))
  "Return the number of hill blocks in PLAIN."
  (ceiling (length plain) size))

(defmacro do-hill-blocks ((var plain &optional retval) (&key (size 4))
                          &body body)
  "Evaluate BODY once for each hill cipher block in PLAIN, with VAR bound to
that block."
  (let ((i-var (gensym))
        (plain-var (gensym)))
    `(let ((,plain-var ,plain))
       (dotimes (,i-var (hill-block-count ,plain-var :size ,size) ,retval)
         (let ((,var (nth-hill-block ,plain-var ,i-var :size ,size)))
           ,@body)))))

(defun make-hill-matrix (key &key (size 4))
  "Get the hill cipher matrix for KEY and block size SIZE."
  (let ((rng (make-xorshift* :seed key)))
    (generate-invertible-matrix size size :rng rng :max 128)))

(defun hill-process (data mat size dict-size)
  "Perform hill encryption or decryption on DATA. Each block is SIZE bytes. MAT
is taken to be the key matrix (or its inverse). DICT-SIZE is the modulus."
  (let ((processed (make-array 0)))
    (do-hill-blocks (bl data (map 'vector (lambda (elem)
                                            (mod elem dict-size))
                                  processed))
        (:size size)
      (setq processed (concatenate 'vector processed
                                   (*mv mat bl))))))

(defun hill-encrypt (plain key &key (size 4) (dict-size 251))
  "Encrypt PLAIN using a hill cipher, using KEY as the secret and SIZE as the
block size. DICT-SIZE is the largest character code in PLAIN."
  (hill-process (map 'vector 'char-code plain)
                (make-hill-matrix key :size size)
                size dict-size))

(defun make-inverse-hill-matrix (key &key (size 4) (dict-size 251))
  "Create an inverse hill matrix for KEY. SIZE and DICT-SIZE are as they are to
`hill-decrypt'."
  (invert-mat (make-hill-matrix key :size size) :modulo dict-size))

(defun hill-remove-padding (plain-vec)
  "Remove trailing null bytes from PLAIN-VEC."
  (loop for i downfrom (1- (length plain-vec))
        while (and (>= i 0) (eq (aref plain-vec i) #\Nul))
        finally (return (subseq plain-vec 0 (1+ i)))))

(defun hill-decrypt (cipher key &key (size 4) (dict-size 251))
  "Perform the inverse operation of `hill-encrypt', restoring the plaintext
of CIPHER that was encrypted using KEY with a block size of SIZE. DICT-SIZE must
be the same value originally used to encrypt the text."
  (coerce (hill-remove-padding
           (map 'vector 'code-char
                (hill-process cipher
                              (make-inverse-hill-matrix key :size size
                                                            :dict-size dict-size)
                              size dict-size)))
          'string))

;; Local Variables:
;; jinx-local-words: "adjugate cofactor cofactors det plaintext prandom unnessesary xorshift"
;; End: