Projects : gbw-signer : gbw-signer_usrbin
| 1 | ;;;; bignum.scm: A user-level unsigned bignum arithmetic library |
| 2 | ;;; J. Welsh, August 2017 |
| 3 | ;;; Trimmed for gbw, March 2020 |
| 4 | |
| 5 | ;; A bignum is a list of words, least significant first. It must not have trailing zeros. Thus each number has a unique representation, and zero is the empty list. |
| 6 | |
| 7 | (lambda () |
| 8 | |
| 9 | ;;; Constants |
| 10 | (define base-nibbles (quotient *fixnum-width* 8)) |
| 11 | (define base-bits (delay (fx*/wrap base-nibbles 4))) |
| 12 | (define base/2 (delay (expt 2 (- base-bits 1)))) |
| 13 | (define base (delay (* 2 base/2))) ;; must be <= sqrt of largest fixnum |
| 14 | (define base-1 (delay (- base 1))) |
| 15 | (define neg-base-bits (delay (- base-bits))) |
| 16 | |
| 17 | (define - fx-/wrap) |
| 18 | (define + fx+/wrap) |
| 19 | (define * fx*/wrap) |
| 20 | (define (zero? x) (fx= x 0)) |
| 21 | (define (even? x) (fx= (fxand x 1) 0)) |
| 22 | (define = fx=) |
| 23 | (define < fx<) |
| 24 | (define <= fx<=) |
| 25 | (define hex "0123456789abcdef") |
| 26 | (define char0 (char->integer #\0)) |
| 27 | (define char10-A (fx-/wrap 10 (char->integer #\A))) |
| 28 | (define bn0 '()) |
| 29 | (define bn1 '(1)) |
| 30 | |
| 31 | ;;; Helpers |
| 32 | |
| 33 | (define (fix->hex n) ;; note 0 -> empty string |
| 34 | (do ((n n (fxshift n -4)) |
| 35 | (acc '() (cons (string-ref hex (fxand n 15)) acc))) |
| 36 | ((zero? n) (list->string acc)))) |
| 37 | |
| 38 | (define (hexdigit->fix c) |
| 39 | (if (char-numeric? c) (- (char->integer c) char0) |
| 40 | (let ((i (+ (char->integer (char-upcase c)) char10-A))) |
| 41 | (if (and (<= 10 i) (< i 16)) i |
| 42 | (error "bad hex digit:" c))))) |
| 43 | |
| 44 | (define (decdigit->fix c) |
| 45 | (if (char-numeric? c) (- (char->integer c) char0) |
| 46 | (error "bad decimal digit:" c))) |
| 47 | |
| 48 | (define (left-pad s len char) |
| 49 | (string-append (make-string (- len (string-length s)) char) s)) |
| 50 | |
| 51 | (define (bn-pad-word->hex w) |
| 52 | (left-pad (fix->hex w) base-nibbles #\0)) |
| 53 | |
| 54 | (define (word->bn w) |
| 55 | (if (zero? w) '() (list w))) |
| 56 | |
| 57 | ;; Construct bignum from big-endian, vector-like sequence of nibbles. Ugly, but linear time. |
| 58 | (define (nibbles->bn nibble-ref len) |
| 59 | (define (loop-words start acc) |
| 60 | (if (= start len) acc |
| 61 | (let* ((next (+ start base-nibbles)) |
| 62 | (word (get-word start (- next 1)))) |
| 63 | (loop-words next (if (and (null? acc) (zero? word)) acc |
| 64 | (cons word acc)))))) |
| 65 | (define (get-word start stop) |
| 66 | (define (loop start acc) |
| 67 | (if (< stop start) acc |
| 68 | (loop (+ start 1) (+ (* 16 acc) (nibble-ref start))))) |
| 69 | (loop (+ start 1) (nibble-ref start))) |
| 70 | (if (zero? len) '() |
| 71 | (let* ((msw-end (remainder (- len 1) base-nibbles)) |
| 72 | (msw (get-word 0 msw-end))) |
| 73 | (loop-words (+ msw-end 1) (word->bn msw))))) |
| 74 | |
| 75 | ;; Rather than "shift left/right", which unnecessarily invoke endianness, I'm using "shift" for multiplications and "unshift" for divisions. |
| 76 | (define (shift-words a k) |
| 77 | (if (null? a) a (shift-words-nz a k))) |
| 78 | |
| 79 | (define (shift-words-nz a k) |
| 80 | (if (zero? k) a |
| 81 | (shift-words-nz (cons 0 a) (- k 1)))) |
| 82 | |
| 83 | (define (unshift-words a k) |
| 84 | (if (or (zero? k) (null? a)) a |
| 85 | (unshift-words (cdr a) (- k 1)))) |
| 86 | |
| 87 | ;;; Type conversion |
| 88 | |
| 89 | (define (bn->hex n) |
| 90 | (let ((n (reverse n))) |
| 91 | (if (null? n) "0" |
| 92 | (apply string-append (fix->hex (car n)) ;; optimize? |
| 93 | (map bn-pad-word->hex (cdr n)))))) |
| 94 | |
| 95 | (define (hex->bn s) |
| 96 | (nibbles->bn (lambda (k) (hexdigit->fix (string-ref s k))) |
| 97 | (string-length s))) |
| 98 | |
| 99 | (define (bytes->bn v) |
| 100 | (nibbles->bn (lambda (k) |
| 101 | (let ((byte (vector-ref v (fxshift k -1)))) |
| 102 | (if (even? k) |
| 103 | (fxshift byte -4) ;; big-endian |
| 104 | (fxand byte 15)))) |
| 105 | (* 2 (vector-length v)))) |
| 106 | |
| 107 | ;; ~Cubic algorithm! |
| 108 | (define (bn->dec n) |
| 109 | (let loop ((n n) (acc '())) |
| 110 | (if (null? n) (list->string acc) |
| 111 | (bn-divrem n '(10) |
| 112 | (lambda (q r) |
| 113 | (loop q (cons (string-ref hex (bn->fix r)) acc))))))) |
| 114 | |
| 115 | ;; Quadratic algorithm! |
| 116 | (define (dec->bn s) |
| 117 | (do ((i 0 (+ i 1)) |
| 118 | (acc bn0 (bn+fix (bn*fix acc 10) |
| 119 | (decdigit->fix (string-ref s i))))) |
| 120 | ((= i (string-length s)) acc))) |
| 121 | |
| 122 | ;; Can overflow (obviously) |
| 123 | (define (bn->fix n) |
| 124 | (do ((n (reverse n) (cdr n)) |
| 125 | (acc 0 (+ (* acc base) (car n)))) |
| 126 | ((null? n) acc))) |
| 127 | |
| 128 | (define (fix->bn n) |
| 129 | (do ((n n (fxshift n neg-base-bits)) |
| 130 | (acc '() (cons (fxand n base-1) acc))) |
| 131 | ((zero? n) (reverse acc)))) |
| 132 | |
| 133 | ;;; Predicates |
| 134 | |
| 135 | (define bn-zero? null?) |
| 136 | |
| 137 | (define (bn-even? a) (or (null? a) (even? (car a)))) |
| 138 | (define (bn-odd? a) (not (bn-even? a))) |
| 139 | |
| 140 | (define (cmp a b) |
| 141 | (cond ((< a b) -1) |
| 142 | ((< b a) 1) |
| 143 | (else 0))) |
| 144 | |
| 145 | (define (bn-cmp a b) |
| 146 | (cond ((null? a) (if (null? b) 0 -1)) |
| 147 | ((null? b) 1) |
| 148 | (else (let ((c (bn-cmp (cdr a) (cdr b)))) |
| 149 | (if (zero? c) (cmp (car a) (car b)) |
| 150 | c))))) |
| 151 | |
| 152 | (define bn= equal?) |
| 153 | (define (bn< a b) (< (bn-cmp a b) 0)) |
| 154 | (define (bn> a b) (< 0 (bn-cmp a b))) |
| 155 | (define (bn<= a b) (<= (bn-cmp a b) 0)) |
| 156 | (define (bn>= a b) (<= 0 (bn-cmp a b))) |
| 157 | |
| 158 | ;;; Addition |
| 159 | |
| 160 | (define (bn+1 a) |
| 161 | (if (null? a) bn1 |
| 162 | (let ((head (car a))) |
| 163 | (if (= head base-1) |
| 164 | (cons 0 (bn+1 (cdr a))) |
| 165 | (cons (+ head 1) (cdr a)))))) |
| 166 | |
| 167 | (define (bn+ a b) |
| 168 | (cond ((null? a) b) |
| 169 | ((null? b) a) |
| 170 | (else (let ((sum (+ (car a) (car b)))) |
| 171 | (if (< sum base) |
| 172 | (cons sum (bn+ (cdr a) (cdr b))) |
| 173 | (cons (- sum base) (bn+carry (cdr a) (cdr b)))))))) |
| 174 | |
| 175 | (define (bn+carry a b) |
| 176 | (cond ((null? a) (bn+1 b)) |
| 177 | ((null? b) (bn+1 a)) |
| 178 | (else (let ((sum (+ (car a) (car b) 1))) |
| 179 | (if (< sum base) |
| 180 | (cons sum (bn+ (cdr a) (cdr b))) |
| 181 | (cons (- sum base) (bn+carry (cdr a) (cdr b)))))))) |
| 182 | |
| 183 | ;; CAUTION: assumes 0 <= b < base |
| 184 | (define (bn+fix a b) |
| 185 | (cond ((zero? b) a) |
| 186 | ((null? a) (list b)) |
| 187 | (else (let ((sum (+ (car a) b))) |
| 188 | (if (< sum base) |
| 189 | (cons sum (cdr a)) |
| 190 | (cons (- sum base) (bn+1 (cdr a)))))))) |
| 191 | |
| 192 | ;;; Subtraction |
| 193 | |
| 194 | (define (bn-1 a) |
| 195 | (if (null? a) (error "bn-1: subtract from zero")) |
| 196 | (let ((head (car a)) (tail (cdr a))) |
| 197 | (cond ((zero? head) (cons base-1 (bn-1 tail))) |
| 198 | ((and (= head 1) (null? tail)) '()) |
| 199 | (else (cons (- head 1) tail))))) |
| 200 | |
| 201 | (define (bn- a b) |
| 202 | (cond ((null? a) (if (null? b) b (error "bn-: subtract from zero"))) |
| 203 | ((null? b) a) |
| 204 | (else (let ((diff (- (car a) (car b)))) |
| 205 | (if (< diff 0) |
| 206 | (cons (+ diff base) (bn-sub-borrow (cdr a) (cdr b))) |
| 207 | (let ((tail (bn- (cdr a) (cdr b)))) |
| 208 | (if (and (= diff 0) (null? tail)) '() |
| 209 | (cons diff tail)))))))) |
| 210 | |
| 211 | (define (bn-sub-borrow a b) |
| 212 | (cond ((null? a) (error "bn-: subtract from zero")) |
| 213 | ((null? b) (bn-1 a)) |
| 214 | (else (let ((diff (- (car a) (car b) 1))) |
| 215 | (if (< diff 0) |
| 216 | (let ((tail (bn-sub-borrow (cdr a) (cdr b))) |
| 217 | (diff (+ diff base))) |
| 218 | (if (and (= diff 0) (null? tail)) '() |
| 219 | (cons diff tail))) |
| 220 | (let ((tail (bn- (cdr a) (cdr b)))) |
| 221 | (if (and (= diff 0) (null? tail)) '() |
| 222 | (cons diff tail)))))))) |
| 223 | |
| 224 | ;;; Multiplication |
| 225 | |
| 226 | (define (bn*2 a) |
| 227 | (if (null? a) '() |
| 228 | (let ((product (* (car a) 2))) |
| 229 | (if (< product base) |
| 230 | (cons product (bn*2 (cdr a))) |
| 231 | (cons (- product base) (bn*2+carry (cdr a))))))) |
| 232 | |
| 233 | (define (bn*2+carry a) |
| 234 | (if (null? a) bn1 |
| 235 | (let ((product (+ (* (car a) 2) 1))) |
| 236 | (if (< product base) |
| 237 | (cons product (bn*2 (cdr a))) |
| 238 | (cons (- product base) (bn*2+carry (cdr a))))))) |
| 239 | |
| 240 | ;; CAUTION: assumes 0 <= scale < base |
| 241 | (define (bn*fix a scale) |
| 242 | (if (or (null? a) (zero? scale)) '() |
| 243 | (let ((product (* (car a) scale))) |
| 244 | (if (< product base) |
| 245 | (cons product (bn*fix (cdr a) scale)) |
| 246 | (cons (fxand product base-1) |
| 247 | (bn*fix+carry (cdr a) scale |
| 248 | (fxshift product neg-base-bits))))))) |
| 249 | |
| 250 | (define (bn*fix+carry a scale carry) |
| 251 | (if (or (null? a) (zero? scale)) (list carry) |
| 252 | (let ((product (+ (* (car a) scale) carry))) |
| 253 | (if (< product base) |
| 254 | (cons product (bn*fix (cdr a) scale)) |
| 255 | (cons (fxand product base-1) |
| 256 | (bn*fix+carry (cdr a) scale |
| 257 | (fxshift product neg-base-bits))))))) |
| 258 | |
| 259 | (define (bn-shift a bits) |
| 260 | (cond ((< bits 0) (error "bn-shift: negative bits")) |
| 261 | ((null? a) a) |
| 262 | (else (bn*fix (shift-words-nz a (quotient bits base-bits)) |
| 263 | (expt 2 (remainder bits base-bits)))))) |
| 264 | |
| 265 | (define (bn* a b) |
| 266 | (define (a* b) |
| 267 | (if (null? b) b |
| 268 | (bn+ (bn*fix a (car b)) |
| 269 | (shift-words (a* (cdr b)) 1)))) |
| 270 | (if (null? a) a (a* b))) |
| 271 | |
| 272 | ;; Still quadratic, but ~30% faster than generic multiplication |
| 273 | (define (bn^2 a) |
| 274 | (if (null? a) '() |
| 275 | (let* ((hd (car a)) |
| 276 | (tl (cdr a)) |
| 277 | (hd^2 (* hd hd)) |
| 278 | (hd^2 (if (< hd^2 base) (word->bn hd^2) |
| 279 | (list (fxand hd^2 base-1) |
| 280 | (fxshift hd^2 neg-base-bits))))) |
| 281 | (if (null? tl) hd^2 |
| 282 | (bn+ hd^2 |
| 283 | (cons 0 (bn+ (cons 0 (bn^2 tl)) |
| 284 | (bn*fix (bn*2 tl) hd)))))))) |
| 285 | |
| 286 | (define (strip-leading-zeros l) |
| 287 | (cond ((null? l) l) |
| 288 | ((zero? (car l)) (strip-leading-zeros (cdr l))) |
| 289 | (else l))) |
| 290 | |
| 291 | (define (bn-split a k cont) |
| 292 | (do ((head '() (cons (car tail) head)) |
| 293 | (tail a (cdr tail)) |
| 294 | (k k (- k 1))) |
| 295 | ((or (null? tail) (zero? k)) |
| 296 | (cont (reverse (strip-leading-zeros head)) tail)))) |
| 297 | |
| 298 | ;;; Division |
| 299 | |
| 300 | (define (bn/2 a) |
| 301 | (if (null? a) a |
| 302 | (cdr (bn*fix a base/2)))) |
| 303 | |
| 304 | (define (bn-unshift a bits) |
| 305 | (if (< bits 0) (error "bn-unshift: negative bits")) |
| 306 | (let* ((full-words (quotient bits base-bits)) |
| 307 | (extra-bits (remainder bits base-bits)) |
| 308 | (a (unshift-words a full-words))) |
| 309 | (if (or (null? a) (zero? extra-bits)) a |
| 310 | (cdr (bn*fix a (expt 2 (- base-bits extra-bits))))))) |
| 311 | |
| 312 | (define (num-bit-shifts start target) ;; optimize? |
| 313 | (define (loop start n) |
| 314 | (if (<= target start) n |
| 315 | (loop (* start 2) (+ n 1)))) |
| 316 | (loop start 0)) |
| 317 | |
| 318 | (define (last l) |
| 319 | (if (null? (cdr l)) (car l) |
| 320 | (last (cdr l)))) |
| 321 | |
| 322 | (define (bn-divrem a b return) |
| 323 | (if (null? b) (error "division by zero")) |
| 324 | ;; Normalize: most sig. bit of most sig. word of divisor must be 1 |
| 325 | (let ((msw (last b))) |
| 326 | (if (<= base/2 msw) (div-normalized a b return) |
| 327 | (let ((s (expt 2 (num-bit-shifts msw base/2)))) |
| 328 | (div-normalized |
| 329 | (bn*fix a s) (bn*fix b s) ;; optimize |
| 330 | (lambda (q r) |
| 331 | (return q (if (null? r) r |
| 332 | (cdr (bn*fix r (quotient base s))))))))))) |
| 333 | |
| 334 | (define (divrem-q q r) q) |
| 335 | (define (divrem-r q r) r) |
| 336 | (define (bn-quotient a b) (bn-divrem a b divrem-q)) ;; optimize? |
| 337 | (define (bn-remainder a b) (bn-divrem a b divrem-r)) |
| 338 | |
| 339 | (define (slice-2 a k) |
| 340 | (if (null? a) 0 |
| 341 | (let ((tail (cdr a))) |
| 342 | (if (zero? k) |
| 343 | (if (null? tail) (car a) |
| 344 | (+ (car a) (* (car tail) base))) |
| 345 | (slice-2 tail (- k 1)))))) |
| 346 | |
| 347 | (define (most-sig-word a) ;; assumes not null |
| 348 | (define (loop a tail) |
| 349 | (if (null? tail) (car a) |
| 350 | (loop tail (cdr tail)))) |
| 351 | (loop a (cdr a))) |
| 352 | |
| 353 | (define (div-normalized A B return) |
| 354 | (let* ((n (length B)) |
| 355 | (m (- (length A) n))) |
| 356 | (if (< m 0) (return '() A) |
| 357 | (let ((n-1 (- n 1)) |
| 358 | (B-msw (most-sig-word B))) |
| 359 | (define (get-qi i Q A) |
| 360 | (define (try-qi qi) |
| 361 | (let ((prod (shift-words (bn*fix B qi) i))) |
| 362 | (if (bn< A prod) |
| 363 | (try-qi (- qi 1)) |
| 364 | (get-qi (- i 1) (cons qi Q) (bn- A prod))))) |
| 365 | (if (< i 0) (return Q A) |
| 366 | (try-qi (min base-1 (quotient (slice-2 A (+ n-1 i)) |
| 367 | B-msw))))) |
| 368 | (let ((B-shift (shift-words B m))) |
| 369 | (if (bn< A B-shift) |
| 370 | (get-qi (- m 1) '() A) |
| 371 | (get-qi (- m 1) '(1) (bn- A B-shift)))))))) |
| 372 | |
| 373 | ;; Multiplicative inverse of a mod n: |
| 374 | ;; (bn-remainder (bn* a (bn-mod-inverse a n)) n) -> bn1 |
| 375 | ;; Assumes reduced input (0 < a < n) |
| 376 | (define (bn-mod-inverse a n) |
| 377 | ;; Extended Euclidean algorithm: find x where ax + by = gcd(a, b) |
| 378 | ;; If the gcd is 1, it follows that ax == 1 mod b |
| 379 | ;; See [HAC] Algorithm 2.142 / 2.107 |
| 380 | ;; Simplified / adjusted for unsigned bignums |
| 381 | ;; Invariants: |
| 382 | ;; [1] 0 <= r < b (loop terminates when b would reach 0) |
| 383 | ;; [2] 0 < b < a (by [1], as b is last r and a is last b) |
| 384 | ;; [3] q > 0 (by [2]) |
| 385 | ;; [4] If x > 0, then last x <= 0 and next x < 0 |
| 386 | ;; If x < 0, then last x > 0 and next x > 0 |
| 387 | ;; (by [3], as next x is last x - qx) |
| 388 | ;; Full proofs (mostly by induction) left as an exercise to the reader. |
| 389 | (define (loop a b x neg last-x) |
| 390 | (bn-divrem a b |
| 391 | (lambda (q r) |
| 392 | (if (bn-zero? r) |
| 393 | (if (bn= b bn1) (if neg (bn- n x) x) |
| 394 | (error "not invertible (modulus not prime?)")) |
| 395 | (loop b r (bn+ last-x (bn* q x)) (not neg) x))))) |
| 396 | (loop n a bn1 #f bn0)) |
| 397 | |
| 398 | ;;; Misc |
| 399 | |
| 400 | ;; Number of significant bits |
| 401 | ;; = least integer b such that 2^b > a |
| 402 | ;; = ceil(log_2(a+1)) |
| 403 | (define (bn-bits a) |
| 404 | (if (null? a) 0 |
| 405 | (+ (num-bit-shifts 1 (+ (last a) 1)) |
| 406 | (* (- (length a) 1) base-bits)))) |
| 407 | |
| 408 | (define (read-bytes n port) |
| 409 | (let ((v (make-vector n))) |
| 410 | (do ((k 0 (+ k 1))) ((= k n) v) |
| 411 | (vector-set! v k (char->integer (read-char port)))))) |
| 412 | |
| 413 | (define (ceil-quotient a b) |
| 414 | (quotient (+ a b -1) b)) |
| 415 | |
| 416 | ; Unbiased random integer generator in the interval [0, n) |
| 417 | (define (rand-bn n) |
| 418 | (if (bn-zero? n) (error "rand-bn: zero range")) |
| 419 | ;; Collecting one more byte than strictly necessary avoids cases where a large part of the range is invalid (e.g. n=130) |
| 420 | (let* ((nbytes (+ (ceil-quotient (bn-bits (bn-1 n)) 8) 1)) |
| 421 | (rand-range (bn-shift bn1 (* nbytes 8))) |
| 422 | (unbiased-range (bn- rand-range (bn-remainder rand-range n)))) |
| 423 | (lambda (rng-port) |
| 424 | (define (retry) |
| 425 | (let ((r (bytes->bn (read-bytes nbytes rng-port)))) |
| 426 | (if (bn< r unbiased-range) |
| 427 | (bn-remainder r n) |
| 428 | (retry)))) |
| 429 | (retry)))) |
| 430 | |
| 431 | ;; Deferred initializations |
| 432 | (set! base-bits (force base-bits)) |
| 433 | (set! base/2 (force base/2)) |
| 434 | (set! base (force base)) |
| 435 | (set! base-1 (force base-1)) |
| 436 | (set! neg-base-bits (force neg-base-bits)) |
| 437 | |
| 438 | (export |
| 439 | bn0 bn1 |
| 440 | hexdigit->fix decdigit->fix ;; not strictly bignum ops, but handy |
| 441 | bn->hex hex->bn bytes->bn |
| 442 | bn->dec dec->bn |
| 443 | bn->fix fix->bn |
| 444 | bn-zero? bn-even? bn-odd? bn= bn< bn> bn<= bn>= |
| 445 | bn+1 bn+ bn-1 bn- |
| 446 | bn*2 bn-shift bn*fix bn* bn^2 |
| 447 | bn/2 bn-unshift bn-divrem bn-quotient bn-remainder |
| 448 | bn-mod-inverse |
| 449 | bn-bits |
| 450 | rand-bn)) |